Sediment Quality Mapping Project Final Report.pdf

South Florida Water Management District Final Report

4600004016-WO03 PO No 950

December 2021

Lake Okeechobee Sediment Quality Mapping Project

Contract 4600004016-WO03 PO No 950

Produced for

3301 Gun Club Road West Palm Beach, FL 33406

Produced by

Estuarine Biogeochemistry Laboratory

Whitney Laboratory for Marine Bioscience University of Florida

9505 Ocean Shore Blvd St. Augustine, FL 32080

December 2021

December 2021 i contributing authors

Todd Z. Osborne Tracey B. Schafer Paul Julian II Patrick W. Inglett

Paul Jones

Draft Final Report

EXECUTIVE SUMMARY

Lake Okeechobee is the largest lake in the southeastern United States and a major channel of water transport from central Florida to the Atlantic and Gulf coasts. Runoff from agricultural and urban areas in central Florida are transported down the Kissimmee River and into the lake along with a plethora of sediments, nutrients, and chemicals from upstream. Since the majority of new sediment introduction to the lake is composed of fine particulate mud, some of the new inputs can also be swiftly transported to other portions of the lake. Lake Okeechobee is extremely shallow with an average depth of less than 3 meters that creates a unique environment for wind transport and mixing of sediment. In order to assess the changes in the sediment profile that occur over time due to these variables, Lake Okeechobee sediment has been mapped four times over the past 40 years. Previous surveys have occurred in 1988, 1998, 2006, and 2020 (which is the focus of this report).

To assess sediment change over time, mud depth was measured in situ and surface sediment samples were collected and brought back to the laboratory for further analysis.

Sediment samples were extracted under an anaerobic environment and analyzed for a variety of chemical and physical characteristics. Sediment biogeochemical variables include total phosphorus (TP), nitrogen (TN) and carbon (TC), select metals, and physical parameters.

Additionally, porewater analysis was conducted to assess nutrient and metal concentrations. It was found that overall TP and TN have not varied significantly over the past 40 years, but spatial distribution of these analytes across the lake has changed over time. Mud depth has also varied between sampling periods with a noticeable increase since 2006. This mud is linked to the highest concentrations of metals and nutrients such as phosphorus. Sand, peat, and rock sites (which usually have a thin layer of mud sediments) have much lower concentrations of these analytes. Although it is not possible to measure these characteristics at a fine enough scale to capture all the system wide variability, there are clear spatial changes in biogeochemical patterns over time across the lake. These changes are likely due to intense storm activity that initiates resuspension of fine sediments and aids in redistribution around the lake. The deeper areas of the lake consistently function as settling areas for mud sediments and associated nutrients (N and P).

TABLE OF CONTENTS

EXECUTIVE SUMMARY

TABLE OF CONTENTS

LIST OF FIGURES

LIST OF TABLES

1.0 INTRODUCTION

2.0 METHODS

2.1 Field Methods

2.2 Laboratory Methods

2.3 Statistical Methods

3.0 RESULTS

3.1 Mud Sediment Distribution

3.2 Sediment Biogeochemical Characteristics

3.3 Biogeochemical Characteristics by Sediment Type

3.4 Porewater Bigoechemistry

3.5 Spatial Analysis of Sediment N and P

4.0 DISCUSSION

4.1 Mud Sediment Distribution

4.2 Sediment Biogeochemical Characteristics

4.3 Biogeochemcial Characteristics by Sediment Type

4.4 Porewater Bigoechemistry

4.5 Spatial Analysis of Sediment N and P

5.0 CONCLUSIONS

6.0 REFERENCES

LIST OF FIGURES

Figure 1 – Sediment sampling sites in Lake Okeechobee Figure 2 -Spatial distribution of sediment zones/types based on weighted nearest-neighbor classification Figure 3 – Schematic of piston coring device used to obtain lake sediments Figure 4 – Bar graphs of average mud depth (cm) Figure 5 – Changes in mud thickness (cm) across each sampling period with 10 cm contour lines identified Figure 6 – Change maps of mud thickness (cm) between each survey (1988, 1998, 2006 and 2020) Figure 7 - Mud volume (km3) in Lake Okeechobee based on spatially interpolated data Figure 8 – Bar graphs of sediment average total phosphorus Figure 9 – Bar graphs of sediment average total nitrogen Figure 10 – Bar graphs of sediment average total calcium Figure 11 – Bar graphs of sediment average total iron Figure 12 – Bar graphs of sediment average total magnesium Figure 13 – Bar graphs of sediment average organic matter Figure 14 – Mean Total P content by major substrate type Figure 15 – Mean Total N content by major substrate type Figure 16 – Mean Total Ca content by major substrate type Figure 17 – Mean Total Mg content by major substrate type Figure 18 – Mean Total Al content by major substrate type Figure 19 – Porewater concentrations (g/L) of total phosphorus and soluble reactive phosphorus Figure 20 – Spatial interpolation time-series of sediment total phosphorus across Lake Okeechobee Figure 21 – Spatial interpolation time-series of sediment total nitrogen (across Lake Okeechobee.

Figure 22 – Spatial interpolation time-series of sediment N:P ratio across Lake Okeechobee Figure 23 – Change maps of total phosphorus (mg/kg) between each survey (1988, 1998, 2006, and 2020) Figure 24 – Change maps of total nitrogen (g/kg) between each survey (1988, 1998, 2006, and 2020) Figure 25 - Sediment totals across Lake Okeechobee during the 2020/2021 sampling period for the top 10 cm Figure 26 – Inflow and outflow annual discharge Figure 27 -Time series of spatially average total phosphorus concentration Figure 28 – Time series of net phosphorus sedimentation coefficient for Lake Okeechobee Figure 29- Estimates of P mass in Lake Okeechobee surface sediments for all sediment types . 40 Figure 30 – Spatial interpolation of surface sediment P mass across Lake Okeechobee sediment for each timestep

LIST OF TABLES

Table 1 -Porewater extraction analytes and methods Table 2 -Sediment analytes and methods Table 3 -Sediment residue analytes and methods Table 5 – Summary statistics of porewater characteristics Table 6 - Summary of geostatistical parameters

1.0 INTRODUCTION

Lake Okeechobee, located in the center of peninsular Florida, is a large shallow lake (Figure 1).

It covers an area of approximately 1730 km2 at an average depth of only 2.7 m, making it highly susceptible to wind-driven water movement. Over the latter half of the 20th century, Lake Okeechobee was significantly impacted by cultural eutrophication, primarily through elevated phosphorus (P) loading associated with runoff from intensive agriculture in its watershed (Havens et al., 1996; Havens & James, 2005; Schottler & Engstrom, 2006). One direct consequence of eutrophication has been the substantial accumulation of P enriched sediment throughout the lake basin. In 1988, a comprehensive survey reported that the upper 10 cm of mud sediments within the lake contained an estimated 28,600 metric tons of P (Reddy et al., 1995). Two additional surveys in 1998 and 2006 made similar estimates of sediment P retention within Lake Okeechobee (Fisher et al., 2001) (BEM 2006).

Within large, shallow lakes such as Lake Okeechobee, sediment resuspension driven by wind-induced water movement may play a significant role in nutrient cycling as particulate bound forms of P can be released into the overlying water column (Søndergaard et al., 2003). Given that over 44% of Lake Okeechobee has been found to be covered with P enriched sediments, resuspension of these sediments may be one of the most important mechanisms affecting P availability within this lake ecosystem (James et al., 1995a, b). Additionally, the easily resuspended thixotrophic sediment has been shown to produce up to 50% more nutrient mass as compared to consolidated sediments (Missimer et al., 2021). These changes make creating restoration projects for Lake Okeechobee increasingly difficult and highlight the importance of regular benthic mapping (Canfield et al., 2020).

Studies conducted on the lake since the 1998 survey indicate potential changes to the sediment profile. Between September 2004 and October 2005, three major hurricanes (Frances, Jeanne, and Wilma) passed directly over the lake creating large seiches, waves, and currents.

Observations of the lake indicate that the mud sediments were remixed to depths between 10 and 20 cm (Engstrom & Pollman, 2005) and that these mud sediments were resuspended and transported to other areas of the lake. Since the 2005 study, several other hurricane and storm events have impacted the lake. Additionally, potential changes to lake inputs since 2006 have also been recorded, as shown by a 2014 model of P imports and exports to Lake Okeechobee that estimated that P imports to the lake had decreased by 25% and P storage had increased by 29% (He et al., 2014). To monitor and examine changes to Lake Okeechobee sediments, this project addresses the questions:

1) What is the extent and change in mud sediment throughout Lake Okeechobee?

2) What are the current pools of some forms of P and nitrogen (N) in the sediments?

3) What is the extent and change of sediment nutrients throughout Lake Okeechobee?

This project provides essential information regarding the location and amount of P and N pools within sediments that may potentially impact the health and future management of Lake Okeechobee and is authorized through the Northern Everglades and Estuaries Protection Plan (NEPP) program.

2.0 METHODS

2.1 Field Methods

Lake Okeechobee sediment samples were retrieved from sites corresponding to sampling locations of three previous lake sediment surveys conducted in 1988, 1998, and 2006 (Figure 1).

Locations of sites were provided by the South Florida Water Management District (SFWMD) in Latitude/Longitude format (degree, minute, seconds). These were transformed to decimal degrees by University of Florida (UF) staff for navigation to sites in 2006 and the same sites were visited again for this report. Field sampling was conducted in August-December of 2020 using a 26’ power boat equipped with GPS. Due to prolific marsh vegetation, a subset of samples (n = 25) could not be accessed during August-December with a skiff. These sites were visited in January 2021 via airboat. A portion of marsh sites were excluded all together from the 2020- 2021 sampling as they were difficult to access due to woody vegetation growth in 2005 and thus were determined to be inaccessible (see fringing marsh areas Figure 2). This decision was made prior to contracting for the 2020 sampling effort.

Sediment cores were retrieved using a piston coring device with removable 7.5 cm ID polycarbonate core tubes (Fisher et al, 1992; Figure 3). Cores were retrieved intact with headspace water and sealed immediately for transport back to the Estuarine Biogeochemistry Laboratory (EBL) at the UF’s Whitney Laboratory for Marine Biosciences. Mud depths were determined on each core at time of sampling with the exception of some sites in the center of the lake where mud depth was beyond the depth of the measuring rod due to excessive water depth.

Ekman grabs (6” square) were taken where cores could not be obtained due to the shallow nature of mud over rock or consolidated sand/rock sediments. These samples were sealed in resealable plastic bags and stored on ice along with the core samples. Secchi disk measurements and water depth were determined at each site visited.

Field QA/QC procedures included personnel training, dual GPS position checks, and duplicate coring efforts at select sites. A senior staff member was present for all field sampling efforts.

Field data was recorded on waterproof, numbered page field notebooks. Ten percent of field data were checked by senior staff for transcription error.

Figure 1 – Sediment sampling sites in Lake Okeechobee

Figure 2 -Spatial distribution of sediment zones/types based on weighted nearest-neighbor classification.

Figure 3 – Schematic of piston coring device used to obtain lake sediments.

2.2 Laboratory Methods

Upon return to the laboratory, cores were extruded in an anaerobic hood (N2 atmosphere) to ensure minimal oxygen exposure. Extruded samples were then sent to the Wetland Biogeochemistry Laboratory (WBL) at UF in Gainesville and centrifuged at 10,000 rpm for 10 minutes to extract porewater. Post centrifugation, porewater was removed from centrifuge tubes and filtered at 0.45 um with Millipore® regenerated cellulose filters. Non-acidified samples were analyzed for SRP within 48 hrs. Non-acidified samples were also analyzed for pH and conductivity, while acidified samples were analyzed for Total P, Total N, Ammonia N, Total Al, Ca, Fe, and Mg using standard methods (see Table 1).

Table 1 -Porewater extraction analytes and methods

Analyte Method PQL QA/QC Type pH EPA 150.2 (U.S. EPA 1982b) 0.01 SU D,R,St,B Conductivity EPA 120.1 (U.S. EPA 1982a) 0.1 uS/cm D,R,St,B SRP EPA 365.1 (U.S. EPA 1993c) 0.005 mg/l D,R,St,B,MS, Total Dissolved P EPA 365.1 (U.S. EPA 1993c) 0.01 mg/l D,R,St,B,MS Total Dissolved N EPA 351.2 (U.S. EPA 1993b) 0.1 mg/l D,R,St,B,MS NH4-N EPA 350.1 (U.S. EPA 1993a) 0.06 mg/l D,R,St,B,MS Total Al EPA 200.7 (U.S. EPA 1994) 0.05 mg/l D,R,St,B,MS Total Ca EPA 200.7 (U.S. EPA 1994) 0.05 mg/l D,R,St,B,MS Total Fe EPA 200.7 (U.S. EPA 1994) 0.05 mg/l D,R,St,B,MS Total Mg EPA 200.7 (U.S. EPA 1994) 0.05 mg/l D,R,St,B,MS

(D denotes duplicates, R repeats, St standards, B blanks, MS matrix spikes)

Sediment samples were first weighed to determine bulk density; subsamples were dried for 3 days at 70°C to determine moisture content, then ground, milled, extracted, and analyzed for Total C, N, and P, Al, Ca, Fe, Mg, inorganic P, and loss on ignition using standard methods (Table 2).

Table 2 -Sediment analytes and methods

Analyte Method PQL QA/QC Type Moisture Content Eq.1 0.01 D,R,B Bulk Density USDA 3.3 (USDA 2004) 0.01 g D,R, Total C EPA 440.0 (Zimmerman 1997) 0.035 g/kg D,R,St,B Total N EPA 440.0 (Zimmerman 1997) 0.035 g/kg D,R,St,B Total P EPA 365.1 (U.S. EPA 1993c) 0.02 mg/kg D,R,St,B,MS Total Inorganic P EPA 365.1 (U.S. EPA 1993c) 0.02 mg/kg D,R,St,B,MS Total Al EPA 200.7 (U.S. EPA 1994) 0.05 mg/kg D,R,St,B,MS Total Ca EPA 200.7 (U.S. EPA 1994) 0.05 mg/kg D,R,St,B,MS Total Fe EPA 200.7 (U.S. EPA 1994) 0.05 mg/kg D,R,St,B,MS Total Mg EPA 200.7 (U.S. EPA 1994) 0.05 mg/kg D,R,St,B,MS

Sediment residue (post-centrifugation sediment) moisture content was calculated and then subsampled to determine KCl extractable P and Ammonia, and NaHCO3 extractable SRP and TP employing standard methods (Table 3).

Table 3 -Sediment residue analytes and methods

Analyte Method PQL QA/QC Type Moisture Content Eq.1 0.01% D,R,B KCl Pi EPA 365.3 (U.S. EPA 1978) 0.01 mg/kg D,R,St,B,MS

KCl NH4-N EPA 350.1 (U.S. EPA 1993a) 0.01 g/kg D,R,St,B,MS NaHCO3- SRP EPA 365.1 (U.S. EPA 1993c) 0.02 mg/kg D,R,St,B,MS

NaHCO3-TP EPA 365.1 (U.S. EPA 1993c) 0.02 mg/kg D,R,St,B,MS

Moisture content was determined by measuring a wet subsample of each soil, then oven drying the soil for 48 hours and taking measurement of the dry soil. Percent moisture content was then calculated (Eq. 1).

Eq. 1: Moisture Content (%) = (Wt Wet soil – Wt oven dried soil)/ Wt oven dried Soil * 100

All analytical procedures conducted in the WBL employ National Institute of Standards and Testing (NIST) traceable standards, analytical and field duplicates, matrix spikes, method blanks and duplicates, and analyte blanks. Personnel must complete required training before they can perform analyses without supervision and all data is QC checked by the Laboratory QC Officer.

Further detail of analytical QA/QC can be found in the project analyte database provided by

WBL/UF.

2.3 Statistical Methods

Kruskal-Wallis was used to examine changes in sediment properties (mud depth, major nutrient concentrations) across survey years. For this analysis, we ignore the confounding effect of sediment type because of the strong concordance between sampling events. Statistical analysis and graphing was implemented in R. Significant differences among years were evaluated using Wilcoxon Rank Sum post-hoc tests in the ggplot 2 R-package (Kassambara, 2020). Statistical significance was set at p = 0.05.

Second, all analyte data was shown to be non-normal upon examination of qq plots, and therefore a Kruskal-Wallis was used to examine differences between nutrient concentrations of different substrate types. Rock, marl, and shell substrates were grouped together and labelled as rock for ease of analysis. Statistical analysis and graphing was implemented in R. Significant differences among years were evaluated using Wilcoxon Rank Sum post-hoc tests. Statistical significance was set at p = 0.05.

To explore spatial patterns of substrate type, mud depth, and biogeochemical attributes across Lake Okeechobee, common geostatistical methods were employed. Sample variograms were evaluated with continuous variables using the variogram function, and variogram models were fit to the sample variogram in the gstat R-package (Benedikt et al., 2016). Each sample variogram was constructed using a width or lag value of 2,700 m and a maximum distance (cutoff) of 40,000 m. Several models were fit to the sample variogram and the models with the best sill ratio and goodness of fit (R2) value were used for ordinary kriging. In cases where model fit was poor, data were log-transformed. The final kriged surface was back-transformed and corrected consistent with Journel (1980) and Yamamoto (2007). The spatial structure for each variogram was evaluated by comparing the partial sill to the total sill. Changes in the strength of spatial structure over time are suggestive of changes in forcing factors affecting the lake.

Due to variability in the data and poor fit of the variograms, mud depth, bulk density, and mineral content spatial distribution were evaluated using thin-plate spline regression (Tps function in the field R-package; Nychka et al., 2017). To visualize sediment zone/type, weighted k-nearest neighbor classification (kknn function in the kknn R-package; Schliep and Hechenbichler, 2016) was used to spatially interpolate sediment zone/type data with a k value (number of nearest neighbors considered) of 10. All interpolations were conducted on a 400m x

400m grid resolution consistent with prior efforts. All geostatistical analyses were performed using R (Ver 4.0.4, R Foundation for Statistical Computing, Vienna, Austria). Spatial models for mud depth, sediment TP and TN, and porewater TP and SRP from historic sampling efforts were recreated using the methods outlined above. Difference maps (2020 minus 1988, 2020 minus 1998, etc.) were constructed for all sampling periods. These maps were created to visualize the magnitude, extent, and pattern of change in key parameters across the lake.

3.0 RESULTS

3.1 Mud Sediment Distribution

Comparison of average mud sediment depths across the entire dataset of the most recent survey indicate that current lake mud depths are statistically different from those observed in the 1998 and 2006 surveys (Figure 4). While observing the differences between mean depths suggests change, a more spatially explicit observation of mud depth distribution provides more insight to the actual changes in the lake sediments (Figure 5) Spatial distribution of mud depth across the lake varied between sampling periods with slight changes in the distribution of mud across the lake (Figures 5 and 6). While mean mud depth varied between years (Figure 4), the zones of the deepest mud remain similar between periods.

Figure 4 – Bar graphs of average mud depth (cm) with error bars representing standard deviation. Kruskal- Wallis test indicated significant differences between years, and a Wilcoxon rank sum indicated that 2020 is significantly different than 1988, 1998, and 2006 (p = 0.0014).

Figure 5 – Changes in mud thickness (cm) across each sampling period with 10 cm contour lines identified.

Difference maps of mud thickness provide a visual representation of change between different sampling efforts (Figure 6). The major changes observed spatially in 2020 are located in the center of the lake and at both the north and south ends. Additionally, total mud volume within the lake has continued to increase from 0.22 km3 in 1988 to 0.31 km3 in 2020 (Figure 7).

Figure 6 – Change maps of mud thickness (cm) between each survey (1988, 1998, 2006 and 2020). The change in mud depth/thickness was calculated as each column subtracted by each row (i.e., 1998 – 1988, 2006 – 1988, 2020 – 1988, etc.).

Figure 7 - Mud volume (km3) in Lake Okeechobee based on spatially interpolated data.

3.2 Sediment Biogeochemical Characteristics

Table 4 provides an overview of sediment biogeochemical characteristics and correlations between selected traits. As with the mud analysis, we examined aggregated evidence for change to control for site level (random) effects. Figures 8 – 12 show the estimated means for sediment totals of P, N Ca, Fe and Mg. For each contrast, the overall model significance (p-value that any year was different from any other) is shown. In the figure legends, the significance of pairwise differences is given.

Table 4 Summary statistics and cross-correlations of measured analytes

M ea n in ax

St d.

D ev

L O

I (

T ot al P

(m g/ kg

H C l-P (m g/ kg

T ot al N

(g /k g)

T ot al C

(g /k g)

T ot al C a (m g/ kg

T ot al M g (m g/ kg

T ot al F e (m g/ kg

LOI (%) 23.5 0.0 93.5 25.1 -

Tot. P (mg/kg) 566.1 12.1 1700.4 488.3 0.23 -

HCl-P (mg/kg) 354.7 6.1 1569.1 372.1 0.05 0.89 -

Tot. N (g/kg) 6.77 0.13 36 8 0.88 0.35 0.15 -

Tot. C (g/kg) 108.3 0 534.0 121.0 0.90 0.33 0.17 0.96 -

Tot. Ca (mg/kg) 41786 106.1 154045 37788 0.09 0.63 0.67 0.10 0.20 -

Tot. Mg (mg/kg) 7572.1 4.3 52417.1 9317.3 0.003 0.49 0.48 0.04 0.15 0.73 -

Tot. Fe (mg/kg) 2682.4 65.2 11480.6 2354.4 0.22 0.90 0.89 0.32 0.31 0.65 0.55 -

Tot. Al (mg/kg) 520.8 13.6 1588.3 423.7 0.30 0.86 0.76 0.35 0.35 0.63 0.55 0.89

Figure 8 – Bar graphs of sediment average total phosphorus (mg/kg) with error bars representing standard deviation. Kruskal-Wallis test indicated no significant differences between years.

Figure 9 – Bar graphs of sediment average total nitrogen (g/kg) with error bars representing standard deviation. Kruskal-Wallis test indicated no significant differences between years.

Figure 10 – Bar graphs of sediment average total calcium (mg/kg) with error bars representing standard deviation. Kruskal-Wallis test indicated significant differences between years (p = 0.04). Wilcoxon rank sum test indicated significant differences between 2006 and 2020 (p = 0.022).

Figure 11 – Bar graphs of sediment average total iron (mg/kg) with error bars representing standard deviation. Kruskal-Wallis test indicated significant differences between years (p = 7.75e-09). Wilcoxon rank sum test indicated significant differences between1988 and 1998 (p= 3.0e-07), 1988 and 2006 (p = 0.024), 1998 and 2006 (p = 9.1e-05), and 1998 and 2020 (p = 1.2e-06).

Figure 12 – Bar graphs of sediment average total magnesium (mg/kg) with error bars representing standard deviation. Kruskal-Wallis test indicated significant differences between years (p = 3.35e-05). Wilcoxon rank sum test indicated significant differences between1988 and 1998 (p= 0.00732), 1988 and 2006 (p = 0.00043), 1998 and 2020 (p = 0.00732), and 2006 and 2020 (p = 0.00043).

3.3 Biogeochemical Characteristics by Sediment Type

Substrate distributions are clearly changing with time in the lake. As such, we might expect changes in the biogeochemistry of the sediment water interface as well. To quantify the magnitude of those changes would require measurements of flux that were not made during this work. However, we can draw some inferences from the analysis of the surface sediments; these include the biogeochemical properties, such as TP, TN, and metal cation concentrations.

(Figures 13 – 18).

Figure 13 – Bar graphs of sediment average organic matter (% LOI) with error bars representing standard deviation. Kruskal-Wallis test indicated significant differences between substrate types (p = 9.806e-10).

Wilcoxon rank sum test indicated significant differences between sand and mud (p= 1.1e-10), sand and peat (p = 5.7e-06), and sand and rock (p = 0.0046).

Figure 14 – Mean Total P content by major substrate type with error bars representing 1 standard deviation from the mean. Kruskal-Wallis test indicated significant differences between substrate types (p < 2.2e-16).

Wilcoxon rank sum test indicated significant differences between mud and peat (p = 4.6e-12), mud and rock (p = 0.0049), mud and sand (p = 2.6e-13), peat and rock (p = 0.0472), and rock and sand (p = 0.0049).

Figure 15 – Mean Total N content (g/kg) by major substrate type with error bars representing 1 standard deviation from the mean. Kruskal-Wallis test indicated significant differences between substrate types (p < 8.27e-10). Wilcoxon rank sum test indicated significant differences between mud and sand (p = 4.1e-11), peat and sand (p = 2.2e-05), and rock and sand (p = 0.0056).

Figure 16 – Mean Total Ca content (mg/kg) by major substrate type with error bars representing 1 standard deviation from the mean. Kruskal-Wallis test indicated significant differences between substrate types (p < 2.2e-16). Wilcoxon rank sum test indicated significant differences between mud and peat (p = 1.6e-15), mud and sand (p = 2.4e-10), peat and rock (p = 0.0056), and rock and sand (p = 0.0345).

Figure 17 – Mean Total Mg content (mg/kg) by major substrate type with error bars representing 1 standard deviation from the mean. Kruskal-Wallis test indicated significant differences between substrate types (p < 2.2e-16). Wilcoxon rank sum test indicated significant differences between mud and peat (p = 4.3e-16), mud and sand (p = 3.7e-12), peat and rock (p = 0.0038), and rock and sand (p = 0.0113).

Figure 18 – Mean Total Al content (mg/kg) by major substrate type with error bars representing 1 standard deviation from the mean. Kruskal-Wallis test indicated significant differences between substrate types (p = 3.268e-16). Wilcoxon rank sum test indicated significant differences between all substrate types, including mud and peat (p = 1.4e-10), mud and rock (p = 0.0424), peat and rock (p = 0.0424), sand and mud (p = 1.7e-12), peat and sand (p = 0.0424), and rock and sand (p = 0.0028).

3.4 Porewater Biogeochemistry

Mean porewater biogeochemical characteristics and parameter correlations are summarized in Table 5. Of note, SRP and TP max values were extremely high at 3,509 and 3,567 mg/L, respectively. These values are believed to be outliers and originated from the same sample; thus they were removed from the geospatial analysis (Figure 19). In the context of previous surveys, porewater TP and SRP are comparatively very low, with SRP falling below 200 mg/L for the entire lake.

Geospatial analysis of porewater SRP (Figure 19) shows an increase from 1988 to 1998, then a significant decline in concentration between 1998 and 2006, which continues into 2020.

However, this trend is not followed in the TP analysis. The 1998 survey suggests an overall decline in porewater TP across the lake from 1988, then a significant increase in 2006 followed by a clear decline by 2020.

Table 5 – Summary statistics of porewater characteristics.

ea n in ax

St d.

D ev pH

C on du ct iv ity

(µ

S)

SR

P

(µ g/

L)

To ta l P (µ

A m m on ia

(m

L

To ta l C a

(m l M g

(m l F e

(m pH 7.49 3.84 8.67 0.57

Conductivity (µS)

951.77 2.25 11127 894.26 0.18 -

SRP (µg/L) 55.22 4.72 3509 280.05 -0.02 -0.01 -

TP(µg/L) 106.37 14.24 3567 280.05 -0.02 -0.03 0.96 -

NH3 (mg/L) 2.96 0 30.35 5.28 0.13 0.02 -0.06 -0.01 - Tot. Ca (mg/L) 138.77 11.34 521.45 97.78 -0.20 0.14 -0.12 -0.2 0.05 -

Tot. Mg (mg/kg) 32.1 1.97 146.79 20.86 -0.07 0.15 -0.11 -0.2 -0.1 0.86 -

Tot. Fe (mg/kg) 0.02 0 0.61 0.06 -0.61 -0.13 0.12 0.17 0.04 0.15 -0.04 -

Tot. Al (mg/kg) 0.13 0 9.54 0.76 -0.07 0.15 -0.01 -0.01 0.06 0.32 0.12 0.82

Figure 19 – Porewater concentrations (µg/L) of total phosphorus (top) and soluble reactive phosphorus (bottom) for 1988, 1998, 2006, and 2020 sampling.

3.5 Spatial Analysis of Sediment N and P

Geostatistical analysis of sediment TP show highest concentrations in the central and eastern portions of the lake and a steady decline in concentration in the sediments of the western portion of the lake (Figure 20). Similar interpolations of sediment TN indicate peak levels of TN in 1998 followed by redistribution and TN reduction in 2020 (Figure 21). Modeling of N:P ratios is also instructive to provide a combined interpretation of both TP and TN changes in time (Figure 22). Over the four surveys, N:P ratios have steadily declined across the lake.

Figure 20 – Spatial interpolation time-series of sediment total phosphorus (mg/kg) across Lake Okeechobee.

Figure 21 – Spatial interpolation time-series of sediment total nitrogen (g/kg) across Lake Okeechobee.

Figure 22 – Spatial interpolation time-series of sediment N:P ratio across Lake Okeechobee.

Table 6 summarizes these parameters for N, P and N:P. The relative structure is similarly strong for N and P, while the structure of N:P is generally slightly weaker. In almost all cases, the relative structure metric is suggestive of moderate spatial dependence (sill ratio = 25 – 75%;

Cambardella et al., 1994). Also, in all cases, there is substantial nugget variability, suggestive of sampling densities and strategies that fail to capture the scale of spatial variability.

Table 6 - Summary of geostatistical parameters describing the spatial structure of the sediment nutrient observations. Relative structure (i.e., sill ratio) is a measure of the variance explained by the semi-variogram.

Model Exp = Expomemtial Sph =spherical)

Year Parameter Model Nugget Sill Range Sill Ratio Log- Transformed

1988 TP Exp 102720.72 4920798.81 946843.63 97.9 N TN Sph 0.86 2.40 20855.05 64.3 Y N:P Exp 0.14 1.71 18479.66 91.6 Y

1998 TP Sph 0.38 1.59 19868.69 76.3 N TN Sph 2.21 4.29 19792.08 48.5 Y N:P Exp 0.89 3.48 14952.28 74.4 Y

2006 TP Sph 0.72 7.30 457082.65 90.1 N TN Exp 0.78 32.66 1141518.46 97.6 Y N:P Exp 0.28 1.82 11079.18 84.8 Y

2020 TP Sph 27757.58 308227.09 40077.22 91.0 N TN Sph 2.69 6.24 37012.81 57.0 Y N:P Sph 0.72 7.30 457082.65 90.1 Y

Figure 23 – Change maps of total phosphorus (mg/kg) between each survey (1988, 1998, 2006, and 2020). The change in total phosphorus was calculated as each column subtracted by each row (i.e., 1998 – 1988, 2006 – 1988, 2020 – 1988, etc.).

Figure 24 – Change maps of total nitrogen (g/kg) between each survey (1988, 1998, 2006, and 2020). The change in total nitrogen was calculated as each column subtracted by each row (i.e., 1998 – 1988, 2006 – 1988, 2020 – 1988, etc.).

Figure 25 - Sediment total phosphorus (TP), HCl extracted inorganic phosphorus (Pi, HCl), total nitrogen (TN), HCl extracted calcium (CaHCl), iron (FeHCl), magnesium (MgHCl) and aluminum (AlHCl) across Lake Okeechobee during the 2020/2021 sampling period for the top 10 cm.

4.0 DISCUSSION

4.1 Mud Sediment Distribution

We compared mud depths across years controlling for random site level effects and observed significant evidence for increased depths in 2020 (Figure 4). Sampling protocols differed between different years, which could manifest in coring depth dissimilarities. Samples for the 1988 and 1998 study were obtained by divers collecting sediment cores, while in 2006 and 2020 they were obtained using a piston corer. However, our general assertion is that the piston corer is likely to be as effective at collecting complete profiles as diver-obtained cores, and that the observed increase in mud depths is not likely an artifact of sampling methodology; particularly since the largest change in mud depths occurred between the last two sample periods (2006 and 2020), which shared the same sampling technique.

In prior surveys, the areas of deepest mud appeared in the center stretching to the north and northeast (i.e., Taylor Creek/Nubbin slough (S-133/S-191) inflow region). However, in this most recent survey, a new zone has appeared extending south (i.e., Pelican Bay/Pahokee; Figure 5).

Changes in mud depth (Figure 6) between years can be informative regarding the potential effects of wind and hurricane induced mixing. Between 1988 and 1998 (Figure 6), the change vector is spatially heterogeneous, with local areas in primarily deeper water showing marked increases juxtaposed with areas in shallower regions showing equally marked decreases in mud depths. Notably, the areas in the center of the lake where mud is currently deepest accumulated large quantities, while outlying areas where mud is largely absent today (western lake) lost appreciable mud depth. Near the inflow of the Kissimmee River (northern lake), mud depths also declined. Overall, the weak spatial structure of the change pattern for 1988-to-1998 is suggestive of the major current gyres that exist in the lake (Jin and Ji, 2004).

In contrast, mud depth changes between 1998 and 2006 occurred over much larger scales. The middle of the lake appeared to have lost appreciable mud, while areas around the central zone appeared to have accumulated mud. Furthermore, mud depth changes between the 2006 to 2020 period (and other periods) show a significant redistribution as well as accumulation of mud.

These temporal dynamics relative to significant disturbance events suggest redeposition processes, perhaps in response to high wind mixing events.

4.2 Sediment Biogeochemical Characteristics

While no significant effect of time was observed for TP and TN (Figures 8 and 9), we did observe a general decline in sediment TP with time. There is no overall trend in TN concentrations with time. For Ca and Mg (Figures 10 and 11), concentrations have remained similar in average and range across the lake across the four sampling periods, with only slight variations. However, Fe exhibited substantially higher concentrations in 1998 than in any other sampling year. The range of Fe values then decreased across the lake in the 2006 sampling period, possibly indicating some temporary change in the late 1990s near the time of the sampling period. Additionally, recent observations were made within the lake that the flocculent structure of the mud sediments was due in part to filamentous magnesium silicates, which would suggest that the flocculent characteristic is abiotic and durable (Harris et al., 2003). Although Mg concentrations have not changed drastically over the years (Figure 12), this is an interesting note that can help in understanding the source and distribution of Mg within the lake. Moreover, during the 2020 sampling effort higher concentrations of nutrients and acid extracted metals were higher in the sediments of the limnetic zone (predominately mud; Figure 25).

4.3 Biogeochemical Characteristics by Sediment Type

Substrate distribution has changed throughout the lake over the past 30 years, and these changes have impacted some biogeochemical analytes, such as TP, TN, and metal cation concentrations (Figures 13 – 18). The mud substrate had the highest TP concentration and amongst the highest concentrations of all of the metal cations measured. TP might be especially high in mud due to the mobility of mud and therefore the regular oxidation of the metal cations from wind and other disturbance events (Moore & Reddy, 1994). This matches previous observations and is likely due to many factors such as metal cation concentration, oxidation status, and organic matter derived from algal biomass (Reddy, 1991; Engstrom & Pollman, 2005). Sand, unsurprisingly, contains the lowest concentrations of all measured analytes due to low surface area and high porosity. Peat is present in the marsh areas of Okeechobee with prolific vegetation. Peat was found to have low concentrations of TP and metal cations, likely due to high rates of vegetative uptake and assimilation. Rock sites (some of which were composed of primarily shell) were shown to have moderate concentrations of metal cations and TP; however, it should be noted that most of the rock sites contained a thin layer of mud that was the substrate analyzed for this report. If mud wasn’t present, it was usually not possible to recover substrate for analysis. TN was variable and low in most substrate types, but it was highest on average in mud sediments (Figure 15). The highest individual values were observed in peat samples, ostensibly due to the organic fraction. Due to the organic nature of the mud sediments, we assert that organic forms of nitrogen persist, providing a nitrogen signature that has not been altered from denitrification or vegetative uptake. The rapid cycling of N and P through repeated algal blooms may also contribute significantly to the organic nature of N in the mud sediments.

4.4 Porewater Biogeochemistry

Among the key management challenges in Lake Okeechobee is the regulation of water column P concentrations; P has been linked to phytoplankton blooms in the lake, and a total maximum daily load (TMDL) for total P has been established to help attenuate this undesirable biological condition. While source reduction is essential to lake recovery, the sediments are a significant source of phosphorus (James et al., 2005; Havens and James, 2005; Fisher et al., 2001) via resuspension and diffusion. To examine the extent to which the sediment is acting as a chemostat for overlying lake water, we examined the spatial and temporal patterns of total and soluble reactive phosphorus in the porewater. The spatial distribution of porewater concentrations (Figure 19) demonstrate two key results, the patterns of both TP and SRP are highly localized and individual points can exert high leverage in the interpolation, resulting in high regional variability. Additionally, in the most recent sampling period, both porewater TP and SRP concentrations across the lake were lower by an order of magnitude than prior sampling events.

Geospatial analysis of porewater SRP (Figure 19) shows a significant decline in concentration between 1998 and 2006, which continues into 2020. This suggests that the increase in SRP between 1988 and 1998 was principally due to changes in P form (i.e., more dissolved organic P in 1988). Dissolved organic P was a significant fraction of the porewater P only in 1988.

Evidently there has been a significant decrease in porewater concentrations of TP since the last survey, likely due to increased disturbance from wind events before and during period of sampling..

4.5 Spatial Analysis of Sediment N and P

We observe a substantially concordant spatial pattern and concentration for TP (Figure 20).

Notably, the main TP hot-spot is the central part of the lake where mud dominates the sediment.

However, other zones of high concentrations are notable in other areas across both space and time; for instance, in 2006 high sediment TP was observed near S-77 (Moore Haven lock) and could be linked to the limited backflow events near this structure during extreme rainfall events.

Meanwhile, in 2020, high sediment TP concentrations (the highest across the lake) were observed near the S-308 structure (Port Mayaca lock) and S-133 (Taylor Creek inflow), driven by frequent backflowing events from the C-44 canal at S-308, and high TP loads from the Taylor Creek/Nubbin Slough watershed, respectively. However, the magnitude of this signal is in stark contrast to the porewater observations, underscoring the need to interpret those patterns with substantial caution.

Patterns of sediment total N are more variable (Figure 21). The mud zone was strongly N enriched in 1998, and moderately enriched in 1988 and 2006. A nitrogen hot-spot in the western marshes around Monkey Box and Moonshine Bay emerged in 1998. Meanwhile, sediment total N concentrations observed in 2020 suggest the strongly N enriched zone observed in 2006 has dissipated slightly, though the enriched region in the southern portion of the Lake remains persistent. In all years, there is no evidence from the sediment data that watershed loading is a significant source of sediment N, despite N-loading to the lake increasing considerably in the past two decades. More directly, the inflow point of the Kissimmee River just west of the northern tip of the lake is consistently among the lowest N concentrations anywhere in the lake, likely due to reduced sedimentation of fine mud materials in areas with inflow currents. N fixation, perhaps due to cyanobacteria, is another potential source (Phlips & Ihnat, 1995).

The spatial distribution of the molar ratios of N to P (Figure 22) suggest a strong north-to-south gradient in potential nutrient limitation; we infer nutrient limitation based on stoichiometric ratios, with N:P ratios greater than 16:1 suggestive of P limitation, and lower values suggestive of N limitation. We note that this inference assumes that 16:1 is indicative of growth limitation in this system (Julian et al., 2019) and one or the other is limiting. In the case of Lake Okeechobee, light limitation is the principal regulator of primary production because of a hyper-abundance of both nutrients. Moreover, light limitation across the lake can vary given changes in water levels throughout the seasons in response to water management activities. Changes in N:P molar ratios are unlikely to affect ecological processes at current elevated nutrient levels. In the southern portion of the lake around Clewiston and Belle Glade, and the western marshes, there is massive N enrichment compared with P (N:P > 100). In the northern lake and the central mud zone, the sediments are more suggestive of high P enrichment (except for the mud zone in 1998).

The observation that, with regard to stoichiometry alone, lake sediments are P limited (N:P > 16:1) is surprising given the legacy and magnitude of P loading. We note that inference about nutrient limitation in the lake water column from observations of sediments is somewhat perilous; sediment concentrations are only weakly associated with the water column concentrations (see Figure 20 for P; correlations for N in sediment and porewater are much weaker).

Nutrient change maps highlight change patterns not immediately evident from the interpolated maps of concentrations. For example, between 1988 and 1998, there was an increase in TP through a central band of the lake, but a general decrease elsewhere; the magnitude of change was moderate in most areas (absolute change value ~ 100 mg/kg) (Figure 23). Between 1998 and 2006, the changes were more localized (i.e., no lake scale gradients) but of equal or higher magnitude. Changes in TP between 2006 to 2020 were both localized (western marsh near S-77 and northern portion near Taylor Creek) and regional, with decreases of sediment TP along the northeast littoral zone and increases in the northern and southern portions. For example, the northern tip of the lake and the sediments near the edge of the marsh fringe (e.g., near Kings Bar) saw large TP increases of more than 500 mg/kg, while the marshes and bays in the western lake and the sediments near Pelican Bay in the southeastern lake saw declines of TP of > 250 mg/kg.

Notably, TP change vectors are not strongly spatially structured, which is suggestive of two conclusions. First, large-scale effects of hurricanes or other lake-scale forcing functions on TP pattern have not been significant. Changes in TP appear more localized, driven by processes operating at scales below our sampling frame. Second, localized patterns of significant change are somewhat suggestive of an insufficient temporal and spatial density of sampling. High spatial leverage in change maps suggests that whatever processes are driving change are occurring at spatial scales smaller than the resolution of this survey. Moreover, the inconsistent temporal patterns suggest that sediments may be changing at time scales shorter than the 10-year monitoring frequency. In either case, there is a need to establish what the appropriate monitoring density in time and space is via short time scale, high density sampling efforts.

A similar change map for N tells a dramatically different story. The changes observed between sampling dates are strongly spatially structured (generally low leverage of individual sampling points in the interpolated surface). Changes between 1988 and 1998 (Figure 24) suggest a moderate but systematic increase in N in the western marshes and the central mud zone emanating from Port Mayaca in the east (> 2,000 mg/kg), accompanied by a general decline in N concentrations in the northern and southern lake zones (> 2,000 mg/kg). In the ensuing 10 years, the trend reverses dramatically, with the western marshes and central mud zone exhibiting an equivalent in magnitude and also largely uniform decline in N, while the northern and southern lake zones saw spatially uniform increases. Meanwhile, between 1998, 2006 and 2020, sediment TN concentrations have generally declined across the lake (Figure 21 and 24).

Additionally, spatial maps for comparisons of all measured sediment characteristics from 2020 sampling were generated (Figure 25). TP and FEHCl appear to have similarities spatially, which is additionally shown by a correlation R2 value of 0.9 (Table 4). Since oxidized Fe will readily bind P and has been previously observed in Lake Okeechobee, it is unsurprising that these would appear well correlated (Moore & Reddy 1994). TN does not appear to show close spatial relationships with the other variables.

Geostatistical interpolation methods allow characterization of various attributes of spatial structure. Paramount among them are descriptors of how much spatial variance (semi-variance) is explained by the model semi-variograms. In particular, the nugget semi-variance, which characterizes the zero-separation variability, can be combined with the partial sill semi-variance, which is that fraction of the variability that can be described using geography, to establish the sill semi-variance, or the variability expected when samples are spatially independent. The ratio of the partial sill to the sill is the fraction of the total variability that is spatially dependent and is a useful metric of structure. Another useful parameter extracted from semi-variance analyses is the range of autocorrelation, which quantifies the distance over which samples are expected to be dependent. High spatial structure and long ranges are both suggestive of exogenous gradients driving spatial self-organization, while short range and weak structure suggest a dominance of local or autogenic processes. That being said, the geostatistical model results presented in Table 6 show a relative robust model results with high sill ratios for P and N:P interpolations. Sill ratios for N were generally acceptable but lower than those of P and N:P despite transformation.

While trends in TP and TN inflow loads to the lake have not significantly changed over the long term, in the past decade and a half loading has gradually increased (Figure 26). The load conditions to Lake Okeechobee have facilitated a strong internal loading as suggested by the differences between the spatially averaged water column P and predicted P concentrations using the Vollenweider model (Figure 27). Moreover, the annual net settling estimate suggests that lake sediments behave as a source of nutrients to the overlying water column.

Figure 26 – Inflow and outflow annual discharge (top), total phosphorus and total nitrogen load (middle), and flow-weighted mean concentrations (bottom) for Lake Okeechobee.

Figure 27 -Time series of spatially average total phosphorus concentration (µg/L) across Lake Okeechobee compared to the simple Vollenweider net settling model for Lake Okeechobee. Data utilize same subset as the South Florida Environmental Report (SFER) Chapter 8B (2021), which removes nearshore influences.

Figure 28 – Time series of net phosphorus sedimentation coefficient for Lake Okeechobee.

Figure 29- Estimates of P mass in Lake Okeechobee surface sediments for all sediment types and mud only over the four monitoring periods.

Figure 30 – Spatial interpolation of surface sediment P mass across Lake Okeechobee sediment for each timestep.

In 1988 it was estimated that the upper 10 cm of mud sediments (i.e. surface sedimented) within the lake contained an estimated 28,600 metric tons of P (Reddy, et al., 1995). To provide a change over time estimate, TP estimate mass was recalculated using available data. The recalculated estimates of mass P storage for the 1988 survey differed. We assume that Reddy et al., (1995) provided a total mass across the lake rather than a mud only, based on the recalculated values. The recalculated P mass for all sediment types was approximately 33% different than the estimate provided by Reddy et al., (1995). This variation could be due to several factors involved with the calculation across a large area. It was assumed that the mean storage (g m-2) for each sediment type was multiplied by the area of each sediments type across the lake. This method was applied to all four sampling periods. Differences in mean TP concentration, bulk density, sediment depth and mud area values all can contribute the disparity. Regardless P mass in surface sediment for mud only and all sediment types have steadily increased from 1988 to 2006 with a 20% and 13% decline in total and mud only P mass (Figure 29). This peak in mass storage in the Lake could be driven by variability in surface sediment storage across the Lake and changes in the extent of mud sediments across the lake over time (Figure 30). Based on the mean P storage for mud sediment (20.89 g m-2) and an area of 802 km2 the 2020 estimated surface sediment P mass in mud sediments is 16,757 metric tons. The cumulative surface sediment P mass across the lake for 2020 is 25,236 metric tons.

5.0 CONCLUSIONS

This project was designed to answer three questions about the status of Lake Okeechobee sediments.