Acoustic Doppler Velocimeter Measurements at Fixed Platforms

In the wetlands the ADV instruments were attached to custom cradles deployed on the upstream site of the fixed platforms at the location of the main monitoring sites referred to as sentinel sites (table 3) near where other hydrologic and biogeochemical measurements were being made (fig. 1). We followed the general operating procedures outlined in Harvey and others (2009) to measure flow velocity at a fixed depth located approximately at the midpoint of the water column. Flow velocity measurements were made continuously for months at a time in which velocities were recorded for one-minute bursts (600 samples) and collected every 15 minutes to save battery power. Quality control on velocity datasets were followed to maintain standards suggested by the instrument manufacturer (SonTek, 2001), as well as specific criteria developed and refined in a prior Everglades study (Riscassi and Schaffranek, 2004). Beam checks were performed on each visit as recommended by the manufacturer. To assess beam check results, the noise count was evaluated. If the noise level of any one beam is outside the range of 60–70 counts, this can be an indicator of an interference with the ADV operation problem by a noisy environment or a larger problem with the electronics. For this case, the system was moved away from any motors or electronics to see if the noise level changed. If no improvement was achieved, the instrument was sent for repair.

To further ensure data quality, ADV velocity data had to be filtered so that samples within a burst passed a 40-percent-minimum correlation filter of signal strength between the three receivers. Within a burst, at least 70 percent of samples need to pass the correlation filter to be retained (Martin and others, 2002). Data were also filtered using an acoustic signal-to-noise ratio (SNR) threshold to eliminate data with SNR values below 5 decibels (dB) and 3.5 dB for ADV and Argonaut ADV, respectively. A sound-speed correction was applied to ADV data using either data from an attached convention temperature device (CTD) or average temperature and salinity recorded by a nearby CTD during the deployment period. The Nortek Vectrino also applied a sound-speed correction using an embedded thermistor and measured salinity.

After filtering, a phase-space threshold “despiking” (that is, removing unrealistic data) and smoothing algorithm (Goring and Nikora, 2002) was applied to ADV, Argonaut ADV, and Vectrino data. Unrealistic data are identified by values that exceed 4 times the standard deviation of a 2-day window from the running median and then removed. Data points identified for smoothing are determined by values that exceed 3 times the standard deviation of a 4-hour window from the running median and then replaced with the running median. These quality-assured data were time averaged to produce hourly and daily values and statistics such as the standard deviation. Flow direction was determined using the Yamartino method (Yamartino, 1984).

We visited the monitoring sites for regular (quarterly) maintenance of the instruments during continuous deployments, during which data were downloaded, batteries were replaced, the compass was calibrated, and diagnostic beam checks were performed on the instrument. On each site visit, the height of probes was adjusted to keep the sampling volume approximately at or slightly below the mid-point of the water column anticipated for the deployment period. Instruments were continuously deployed on average for six months during the study period, between September and March, and we retrieved the instruments from the sites to avoid any potential damage owing to the probe’s exposure to air during the low water-level season. The velocity data collected during times when we visited the monitoring sites were removed to exclude the data from the disturbed system.

Late in the study (July 2021), telemetry systems using a SUTRON SatLink 3 Logger/Transmitter were set up at UB2, DB1, 3B-350E, 3B-1100S, S633GNU-A, and S633GSU-B). which the USGS uses for visualizing and processing downloaded hydrologic time-series data that later are made accessible in the NWIS (USGS, 2016). In the processing of ADV data that was accessed from telemetry using AQUARIUS (Aquatic Informatics, 2024), a threshold was implemented so that flow-speed values below −5 ft/s and above 5 ft/s were suppressed to filter out extremely erroneous data that can result from instrument connection failures. Unfortunately, due to the Covid-19 pandemic, field operations during the last year and a half of the study were negatively affected and the instrument data were never manually downloaded. As a result, flow speed data that were downloaded remotely and processed through AQUARIUS did not retain SNR values and additional quality control information, so processing of data from the last year of the study excluded the correlation and SNR filtering that had been performed on the data collected earlier. Ideally, in the future, ADV data should also be manually downloaded manually, and quality assurance and quality control (QA/QC) should be performed manually using the filtering and despiking steps described earlier in this section to satisfactorily assure the data. Downloading ADV data by SatLink will only be useful as a method for remotely screening data to identify potential problems, until a time when all quality control information can be remotely downloaded with the speed data.

Several methods were used for estimating the flow speed at RS1D-S after the USGS ADV was removed late in the study. For high flows after June 19, 2019, the ADV data from a SFWMD emplacement were used. For high-flow periods with no data, the velocities at RS1D-S were estimated using a linear regression between ADV data and S152 discharge (speed [in cm/s] = 0.0147*(S152 discharge [in ft³/s]) + 0.7893). After June 19, 2019, and when the S-152 culverts were closed, we estimated the resulting low-flow velocities as equal to the average of the measured low-flow velocities between 2010 and 2019 (0.5 cm/s).

Early in the study (2011–15), velocity profiles were measured vertically in the water column at selected sites. For vertical profiling, the instrument cradles were further modified for vertical movement that permitted measurements at multiple depths in the water column. At each depth increment, velocities were measured in 1-minute bursts yielding 600 samples with SNRs monitored continuously to determine if the sample volume was obstructed by vegetation and as an indicator of the vertical location of the top of the floc. Velocity profile data underwent the same corrections and filters as the continuous velocity data. Velocity profiles can be used to refine depth-averaged estimates based on continuously measured velocities measured at approximately the mid-depth. Essentially, point velocities at the mid-depth can be adjusted using velocity profile shape factors described in the procedures documented in Lightbody and Nepf (2006) and in Harvey and others (2009).

FlowTracker Measurements Away from Fixed Platforms

Discrete measurements of flow velocity were made using a SonTek ADV FlowTracker. Those measurements permitted data collection at sites where installation of fixed platforms with ADV cradles was impractical. FlowTracker data collection in the wetland was performed away from any artificial obstructions, such as the airboat. The FlowTracker was positioned in a location with approximately 30 cm of clear space around the probes; vegetation was clipped if necessary. Data were recorded for one minute at the midpoint of water column. Sometimes data were collected at 3 different depths: 5 cm above floc surface, in the middle of the water column, and 5 cm below water surface. Instrument operating procedure performed was followed according to the FlowTracker Handheld ADV Technical Manual (SonTek, 2007). Quality control of the measured data is continuously performed on all variables collected during a measurement. The SNR is the most important quality control data provided by the FlowTracker2, the default threshold of 10 dB is provided in FlowTracker2, User's Manual 1.6 (SonTek, 2019).

Bed Shear Stress Measurements

Bed shear stress measurements were taken pre-release and during the flow release at 6 discrete points along the RS1U ridge-to-slough transect using a Nortek Vectrino II Profiler ADV, manufactured by Nortek. To compute bed shear stress, the Vectrino Profiler was mounted on a vertically sliding rod and deployed to sample flow immediately above the bed. Each instrument reading consisted of instantaneous velocity collected over 18 points within a 2.5-cm vertical profile, at a spacing of 2.0 mm. The instrument was operated at high power to achieve the maximum SNR, with pinging set to “max interval.” Each reading was acquired over a period of 5 minutes. At the end of 5 minutes, the profiler was moved to its next location along the sliding rod (which overlapped with the first by about 0.5 cm) and again deployed. In this manner, four stacked 2.5-cm profiles were obtained near the bed. A fifth 2.5-cm profile was acquired in the middle of the water column at each site.

Raw velocity records were filtered to remove data points that did not meet standards for SNR (at least 5 dB), correlation (at least 40 percent), or consistency in redundant measurements of the z-direction velocity component. Instrument noise spikes were removed using Goring and Nikora’s (2002) threshold despiking algorithm. Following Biron and others (2004), profiles of total stress were calculated as the sum of the total kinetic energy derived stress and laminar stress, as follows:

where is a proportionality constant (0.19), and u’, v’, and w’ are the velocity fluctuations of the streamwise vertical and lateral components where < > denotes a temporal average, μ is dynamic viscosity (in pascal seconds, [Pa·s]), and U is the flow speed (in cm/s) at a height of z above the bed (in cm), where U is given as the resultant of the component-wise velocities, ${\left({〈u〉}^2+{〈v〉}^2+{〈w〉}^2\right)}^{1/2}$. The maximum near-bed total stress was selected as the bed shear stress.

Water-Column Sampling

Water column samples were collected for analysis of total phosphorus concentration (TP) in surface water. During 2013, 2014, and 2015, the water column samples were collected concurrently with SSC samples (see “Suspended Sediment Sampling” section) by peristaltic pumping from the middle of the water column at a rate of 60 mL/min through size 16 Masterflex tubing into a 1-liter acid washed Nalgene sample bottle. To avoid inadvertently capturing periphyton in the sample, the pump tubing had been fitted with a 500-μm “coarse” Nitex screen over the tubing inlet. Samples were preserved by adding sulfuric acid (H₂SO₄) to the sample bottle to reach a pH of 2 and placed on ice before shipment to the SFWMD laboratory for analysis. Samples were acidified with H₂SO₄ to reach a pH of 2 and packed on fresh ice and shipped overnight for analysis at the Wetland Biological Laboratory (WBL), University of Florida, Gainesville, Florida.

Epiphyton and Bed Floc Particulate Sampling

Epiphyton samples were collected that represented fine particulates principally of an organic nature attached to macrophyte stems as well as to leaves of floating aquatic plants (also described as periphyton mats). Epiphyton samples were collected using a wet/dry vacuum operated at approximately the water column’s mid-depth within the canopy of vegetation stems and floating vegetation, as described in Larsen and others (2009a). Epiphyton was dislodged from vegetation by moving the vacuum hose though the canopy and pumping a slurry of epiphyton into a plastic tub. If needed to collect a representative sample the macrophyte stems and floating vegetation were cut and collected in the plastic tub and shaken with some Everglades water to dislodge epiphyton and to create a slurry. The tub was tipped to pour the slurry through a 500-μm Nitex “coarse” filter into a 1-liter Nalgene acid-washed bottle that was allowed to settle and decanted until the sample fit within a 250-mL acid-washed Nalgene bottle.

Bed floc samples were collected by coring the bed with a clear sharpened acrylic tube (1.9-inch inside diameter) to obtain all the top layer of loose, flocculent organic sediment known as “floc” that lies atop more consolidated peat. Floc is important because it represents recently settled material that is biogeochemically active and that often is phosphorus rich. After the core tube was removed it was held upright, allowing time for the floc to settle back onto the more consolidated peat that was also captured in the core tube. Then the clear liquid was gently poured off the top of the core, retaining the slurry of floc above the well consolidated peat. The thickness of floc was recorded by measuring its height with a tape held outside the core tube. Next the slurry of floc was poured off through a 500-μm Nitex “coarse” filter into a 1-liter acid-washed sample bottle that was allowed to settle and decanted until the sample fit within a 250-mL acid-washed Nalgene bottle. Note that after September 27, 2019, bed floc samples were collected by peristaltic pumping of the floc layer without disturbing the more consolidated peat layer through a 500-μm Nitex “coarse” filter into collection bottles. The area of peristaltic pumping (approximately 10×10 cm) was kept to a minimum to collect floc from above the underlying peat in that area without also collecting the peat. The area of collection was expanded as necessary to collect enough sample for analysis. Samples were stored in coolers containing wet ice to maintain a sample temperature of 4 °C. Samples were shipped either the same day or the next day with fresh ice using one-day service to the WBL.

Water-Column Biogeochemical Analysis

The SFWMD analyzed TP concentrations of water samples by first undertaking an acid-persulfate digestion to convert organic and inorganic phosphates into reactive phosphate. Total phosphorus analysis is based on the formation of a blue-colored antimony-phosphomolybdate complex from the reaction of ammonium molybdate, and antimony potassium tartrate with reactive phosphate in an acidic medium followed by reduction with ascorbic acid. The blue color complex was measured by photometric analysis at 880 nanometers (nm) on the Lachat Flow Injection Analyzer (FIA) with a detection limit of 0.002 mg/L and detection range of 0.002 to 0.4 mg/L. The total phosphorus concentration includes phosphorus associated with suspended sediment and dissolved phosphorus in the water column.

Beginning in September 2016, water column samples were analyzed for both TP and total dissolved phosphorus (TDP) by the WBL. Samples were first subjected to an acid-persulfate digestion to hydrolyze all condensed inorganic and oxidize organic phosphates and convert them into orthophosphate. The digestate was then analyzed on a Bran+Luebbe AutoAnalyzer 3 (AA3) instrument following the guidance of standard method (SM) 4500P-F for automated ascorbic acid reduction method for TP or TDP analysis (Lipps and others, 2023). This method is applicable to surface waters, ground waters, and porewater for either TP or TDP. The applicable range for low level TP analysis using AA3 is from 0.002 to 0.15 mg/L, or higher can be measured by using sample dilution. Detection limit for water TP is 0.002 mg TP/L using a standard curve that ranged from 0–0.15 mg/L. For analyses that passed the QA and QC, the particulate phosphorus concentration (TPP) was determined as follows:

Epiphyton and Bed Floc Biogeochemical Analysis

All biogeochemical analyses of epiphyton and bed floc sediments were performed by the WBL. Analysis included: dry weight of solids, loss on ignition (LOI), total nitrogen (TN), total carbon (TC), and labile, microbial, and refractory (total remaining) phosphorus content of epiphyton and floc samples determined by several extractions and a total digestion. For dry weight of solids, samples were dried to a constant weight (detection limit of 0.01 percent), and for LOI, samples were ignited at 550±50 °C in a muffle furnace. Weights were determined using scales and balances verified daily using National Institute of Standards and Technology (NIST) certified weights. TN was found using a Thermo Electron FlashEA 1112. Epiphyton and floc samples underwent a sequence of NaHCO₃ extraction, NaHCO₃ plus chloroform extractions, and finally a perchloric acid digestion using a Technicon Autoanalyzer AA3.

The results from NaHCO₃ and NaHCO₃ plus chloroform extraction were used to determine organic phosphorus fractions consisting of refractory, microbial and labile phosphorus (Ivanoff and others, 1998) as follows:

where L($P_{NaHC{O}_3}$) is the labile phosphorus concentration (in milligrams per kilogram [mg/kg]) removed through extraction with NaHCO₃ (in mg/kg), M is microbial phosphorus concentration (in mg/kg) which is the difference between the phosphorus extracted by NaHCO₃ and NaHCO₃ + Chloroform ($P_{NaHC{O}_3+Chloroform}$), R is refractory phosphorus concentration (in mg/kg), and TP is the total phosphorus concentration (in mg/kg). The subscript “fraction” indicates the fractionation of TP.

Laboratory Analytical Results Assessment

Laboratory results from the WBL were delivered in batches corresponding to sample trips. Each batch of analytical results had been flagged to indicate if reported concentrations were less than the minimum detectable level (MDL), that is, the concentration representing 99 percent confidence that the analyte was detected at a value is greater than zero, and a practical quantitation limit (PQL), that is, the estimated lowest concentration which can be reliably determined within the highest confidence limits of precision and accuracy. The MDL and PQL varied across different sampling runs, and their range is reported in table 6.

Laboratory QA and QC included determination of the accuracy of digestion process and calibration curves using commonly used standards such as WBL, NIST and ERA-Waters Corporation (ERA) certified reference materials for complex nutrients including NIST Peach Leave, WBL Standard Soil, ERA Simple Nutrients for total phosphorus of floc samples. When calibrations were inadequate the laboratory reanalyzed the sample. If the same response was obtained, the laboratory prepared new primary and calibration standards. If that failed, the laboratory checked and re-prepared stock solutions, and instrument service as a final resort. All outcomes were noted in a QC result log.

Further laboratory considerations for flagging included sample holding time, actual sample temperature, and minimum amount of sample for reliable data. The actual flags used by the laboratory included “U,” value less than the MDL (non-detect); “I,” value greater than or equal to MDL but less than PQL; “Q,” actual holding time greater than approved holding time; “T,” actual temperature greater than approved temperature; “F,” QC fail; “D,” delay in receiving sample shipment; “B,” sample amount below minimum required; “N,” not enough sample for analysis. To that list, we added a flag “E,” excessively high concentration considered highly unlikely to be accurate. Those flagged values were identified as more than 2 standard deviations larger than the mean of all samples.

Dissolved Oxygen and Temperature Modeling

Over several years of high-flow experiments (2016–22), water temperature, and dissolved oxygen concentration were monitored on 30-minute intervals with YSI 6920 V2 sondes, and later with EXO3 Multiparameter sondes, for the purpose of modeling aquatic metabolism.

The sensors were deployed using the aluminum pipes supporting boardwalk structures or PVC pipes such that the sensors were located approximately at the midpoint of the water column. The conductance probes were calibrated with 1,000 μS/cm standard. If values were not within ±20 μS/cm, a 1-point calibration was performed. The YSI specific conductance probe has a resolution of 0.001 to 0.1 mS/cm and accuracy ±0.5 percent of reading plus 0.001 mS/cm, and the EXO specific conductance probe has a resolution of 0.0001 to 0.01 mS/cm and accuracy being the greater of either ±1 percent of the reading or 0.002 mS/cm (YSI, Inc., 2022b, d, e).

The optical dissolved oxygen probes were checked for calibration and recalibrated, if necessary, on each site visit. First, the instrument was checked for calibration by comparing the instrument reading in air-saturated water against the expected reading based on the local temperature and “true” barometric pressure. Barometric pressure (BP) readings from the weather service are corrected for elevation and have to be uncorrected using the following equation: “true” uncorrected BP = (corrected BP)−(2.5*[local altitude in feet above sea level/100]). The Everglades are just 9.8 feet above sea level, indicating that the percent saturation in an air-saturated water sample should read approximately 100 percent saturation ±1 percent or 2 percent in an air-saturated water sample. The acceptable uncertainty is selected based on manufacturer specifications for the YSI DO probe, which measured dissolved oxygen (DO) at a resolution of 0.01 mg/L and an accuracy the greater than ±2 percent of readings, or 0.2 mg/L for the range from 0 to 20 mg/L, and for the EXO DO probe, which has an accuracy that is greater than ±0.1 mg/L, or 1 percent of readings for the range from 0 to 20 mg/L (YSI, Inc., 2022c, e84). Calibration checks were performed during each site visit and recalibrations were conducted whenever necessary. It should also be noted that the YSI temperature probe has a resolution of 0.01 °C and an accuracy of ±0.15 °C, and the EXO temperature probe has a resolution of 0.001 °C and an accuracy of ±0.2 °C (YSI, Inc., 2022b, d).

The period and sites where continuous monitoring of dissolved oxygen, temperature, and specific conductivity was conducted is summarized in table 3. Sondes were positioned in the middle of the water column unless noted. Between November 7, 2018, and March 26, 2019, three sondes were deployed at RS1DS to measure conditions in the upper-third, middle-third, and low-third of the water column. During roughly the same period two sondes were deployed at RS1DR to measure in the upper-half and lower-half of the water column.

Late in the study (July 2021), telemetry systems using a SUTRON SatLink 3 Logger/Transmitter were set up at 3B-350E, 3B-1100S, S633GNU-A, and S633GSU-B. In the processing of ADV data that was accessed from telemetry using AQUARIUS (Aquatic Informatics, 2024), thresholds were implemented so that temperature values below 1 °C and above 40 °C were suppressed, specific conductivity values below 20 µS/cm and above 2,000 µS/cm were suppressed, and dissolved oxygen values below −1 mg/L and above 10 mg/L were suppressed to filter out extremely erroneous data that can result from instrument connection failures. Water quality was measured to understand the effects of changes in flow conditions during the high-flow experiments on metabolism at DPM study sites. Eventually, it will help us to predict the effect of the Everglades restoration activities on the aquatic metabolism processes in the Everglades wetland environment.

Temperature Vertical Profile Measurements

Surface-subsurface water interactions are often difficult to detect in wetlands by hydrometric measurements. Attempts have been made in the Everglades using seepage meters (Choi and Harvey, 2000) as well as measuring the hydraulic conductivity by slug tests in piezometers and using it to estimate the vertical flux by combining it with hydraulic head gradient measurements in piezometers using Darcy’s Law (Harvey and others, 2004). Vertical flux across a peat interface also may be detectable by deploying temperature sensors vertically in sediment and modeling heat transport to constrain the vertical flux of water as discussed by Stallman (1965). Traditional piezometers were installed with four temperature sensors in the DPM experimental footprint in the central Everglades at sites RS1D-S, S1, UB2 and Z51-USGS-S. Each piezometer is made of a 0.5-m-long and 1-inch-wide PVC pipe enclosing four TMCx-HD Water/Soil temperature sensors connected to a single HOBO U-12 4-channel data logger (with a resolution of 0.03 °C at 20 °C and an accuracy of ±0.25 °C) (Onset, 2022). The piezometer was inserted into the peat with temperature sensors located at the floc-peat interface (0 cm depth) and three additional depths within the peat sediment. Three depth configurations were used. Configuration 1 had probes located at 0, 5, 15, and 35 cm below the peat surface. Configuration 2 had probes located at 0, 5, 10, and 15 cm below the peat surface. Configuration 3 varied by site and had a probe at the floc-surface water interface, floc-peat interface (0 cm), at half of the peat thickness (69–87 cm), and at the deepest possible point given the length of pipe (133–143 cm). Each sensor collected temperature every 10 minutes. Monitoring periods ranging from 2013 to 2016 of temperature profiles at several sites performed in this way are shown in table 3. Before installation, measurements in a well-mixed, warm bath were used to compare the accuracy of the sensors. For this controlled experiment, the sensors showed an accuracy of approximately 0.05 °C, which is better than the accuracy reported by the manufacturer.

In a later study, HOBO U-12 4-channel data loggers with TMCx water/soil temperature sensors were also used to monitor temperature in more detail at RS1D-S. To install the temperature sensors in the water column, a 1.5-inch PVC pipe was driven into the peat. From January 9, 2018, to February 1, 2018, 8 temperature sensors spaced 7 cm apart, with the bottom-most temperature sensors located 10 cm above the floc surface, were secured to the PVC pipe. The temperature sensor head was positioned to avoid the direct contact with the pipe. On February 1, 2018, one additional temperature sensor was installed 7 cm above the former top-most temperature sensor, as well as one temperature sensor for air temperature. To reduce the effects of solar radiation on temperature measurements, aluminum foil-wrapped thermal shield was installed over the air-temperature sensor. Furthermore, to reduce the effects of thermal conductance along the wires leading to the temperature sensors, the wiring was wrapped with aluminum foil. Temperature profile measurements were taken in this manner from February 1, 2018, to April 10, 2018.

On April 10, 2018, more temperature sensors were installed in the peat and floc. A wooden stake was driven into the peat, with temperature sensors located 120 cm below the peat surface, 75 cm below the peat surface, 25 cm below the peat surface, 5 cm below the peat surface, 1 cm above the peat surface (in floc), 4 cm above the peat surface (in floc), 7 cm above the peat surface (in floc), and 15 cm above the peat surface (in water column). Temperature profile measurements were taken in this manner from April 10, 2018, to March 26, 2019.

Vertical temperature profiles were used to calculate daily Richardson number values at RS1D-S. The Richardson number is a dimensionless number to express the ratio of the buoyancy term to the flow shear term (Rutherford, 1994). The daily maximum water density gradient between adjacent sensors, distance between the sensors at maximum water density gradient, daily average water velocity, and depth to the middle of the water column are needed to calculate the Richardson number.

where is gravitational acceleration (9.81 meters per square second [m/s²]), $\frac{{\rho }_{upper}- {\rho }_{lower}}{z_{ul}}$ is the daily water density gradient (in kilograms per meters raised to the fourth power [kg/m⁴]), which is the difference of water density between the upper and lower adjacent sensors at which the maximum temperature difference occurs divided by the distance between them, ${\rho }_{average}$ is the average of the water densities at the minimum and maximum gradients (in kilograms per cubic meter [kg/m³]), and $\frac{u_{upper}- u_{lower}}{z_{ul}}$ is the water velocity gradient (1/s), which was approximated as the daily average water velocity divided by the mid-depth.

Introduction

Everglades water resources and ecology have been degraded by a loss of patterned microtopography and vegetation communities that are integral to sustaining water resources and high-quality habitats. Returning high flows to the wetlands is central to Everglades restoration, and is expected to improve water storage, flood conveyance, and ecological habitats. Although restoration has the potential to restore lost functions, but with the potential benefits comes an acknowledgement of risks, for example, propagation of contaminants to downstream areas and continued over drying of Everglades wetlands.

There is a recognized need for physically based modeling of Everglades flow dynamics to account for changes in landscape condition. As primary roughness elements that impede flow, vegetation architecture and microtopography are important in determining water velocities and sediment movement as well as water depth in the Everglades. Biophysical modeling is needed to show how changes in Everglades vegetation and microtopography will affect flow dynamics under future conditions of flow restoration with additional challenges of a changing climate.

Here, we implement a simple flow model in the DPM to assess the relative importance of management factors, specifically the discharge and spreading angle of flow releases through culverts, and to compare those with effects of landscape factors that control velocity and depth of sheet flow, as well as the size of the wetland area where bed sediments are entrained in sloughs and deposited on ridges, a condition favorable for maintaining microtopography and open slough habitats. Landscape factors include vegetation community type and architecture that determine hydraulic roughness, as well as the proportions of ridge and slough vegetation communities and their spatial arrangement including the degree of connectivity of the sloughs. We combine the modeling with observations before and after implementation of the Active Marsh Improvement project that increased slough area and connectivity by removing sawgrass that had invaded historical sloughs. We also model scenarios for a well-functioning Everglades with landscape attributes like the historical ridge and slough landscape.

Methods—Radial Flow Dynamics of S-152 High-Flow Releases in DPM

Flow releases between the large, levee-enclosed basins of the Everglades occur actively through gated culverts or passively through gaps in the levee. Active flow releases through gated culverts concentrate the inflows such that they spread radially into the downstream basin (fig. 3). Therefore, a flow model was needed that could track how a radial pattern of flow interacts with downstream vegetation and landscape characteristics.

Governing Equations.—Flow in the DPM is described by a steady-state mass continuity equation that accounts for a radial pattern of water dispersal through the wetland. The continuity equation states that local unit-width discharge, q (in square length over time [L²/T]) declines logarithmically as a function of the radial distance, r in length [L]), total flow through structure, Q (in cubic length over time [L³/T]), and spreading angle, θ (in radians):

and, in terms of the velocity: where U is the vertically integrated flow velocity in length over time [L/T], h is water depth in length [L].

To model the flow in terms of wetland properties that control how the flow spreads, we use the biophysical flow rate expression (BioFRE) of Harvey and Choi (2022a). The biophysical flow rate expression is distinguished from other rate expressions like Manning’s rate law in being parameterized from field measurements rather than calibrated versus field-measured water velocities. Using BioFRE can help anticipate how changing wetland conditions will affect flow and how it is partitioned between a change in flow velocity and a change in flow depth. The BioFRE equation is:

and, in terms of velocity: where n_eff is the upscaled effective roughness coefficient for the wetland landscape (in seconds over meters raised to the 1/3 power [s/m^1/3], S is water slope in meters over meters [m/m]), β is the exponent of water depth, for which a value of 0.667 was used, and α is exponent of water slope, for which a value of 1 was used for consistency with previous simulations in the Everglades.

The BioFRE rate law expression uses the effective upscaled roughness coefficient for the landscape from field measurements. For example, the depth-dependent roughness of ridge and slough vegetation was calculated from a force balance equation (Harvey and others, 2009) based on theory (Nepf, 1999) and measurements of stem frontal area and stem diameter (Larsen and others, 2009b). Furthermore, those vegetation community specific roughnesses are upscaled to represent an effective landscape roughness (n_eff) based on the scaling analysis of Larsen and others (2017a) for which we used estimates of aerial coverage, p, microtopographic height difference between ridge and slough communities (z_p) fractal dimension (f_d), and directional connectivity (DCI), of the slough-ridge landscape pattern (Harvey and others, 2022b). Flow is a power function of water depth and slope for which the exponents, 1 for α and 0.667 for β, were selected to be consistent with previous Everglades modeling (Said and Brown, 2013; Lal and others, 2015).

To use BioFRE expression for radial flow modeling in DPM, the BioFRE equation (28) was substituted into the mass continuity equation (26) to yield a governing equation for wetland flow velocity from a point release into a radial flow field:

where the predicted water depth is:

We used the resulting governing equations in radial coordinates to investigate the relative importance of the wetland variables n_eff and S and the human-managed variables Q and θ in controlling outcomes for U, h, and resulting changes over time, including sediment dynamics, that can determine the effectiveness of high-flow restoration in the Everglades. The total culvert discharge (Q) and the radial spreading angle (θ) are largely under the control of human engineering projects and larger-scale patterns of water availability. In contrast, the hydraulic variables (n_eff, S, and h) depend on vegetation and landscape conditions (for example, ridge proportion and associated microtopography), which can evolve over time because of changes in the managed variables. For example, a change in Q or θ will interact with the effective roughness (n_eff) in ways that adjust water slope (S) and water depth (h). In turn, any change in water depth h will feed back by favoring increases or decreases in vegetation stem density or invasion by new vegetation species (Larsen and Harvey, 2011). Structural changes in vegetation communities will further alter hydraulic variables n_eff, S, and h through biophysical feedbacks that operate on relatively short timescales of just a few years. Because of the interdependencies among variables, we performed many simulations to investigate how various combinations of managed variables and vegetation and landscape conditions can affect future hydrologic and ecological conditions in the Everglades. Although challenging to model, the delicate balance between managed flows and vegetation and landscape conditions is a primary driver of the Everglades ecosystem function and health.

The BioFRE is uniquely capable of simulating changes that could occur in the hydroecology of the Everglades. Unlike the typical hydrologic models used in the Everglades, the effective roughness parameter of BioFRE is calculated from vegetation and landscape measurements rather than being adjusted by trial and error to fit hydrologic observations. As a result, BioFRE can help anticipate future changes brought about by the interdependencies discussed earlier. BioFRE simulations can be verified by comparing simulation against independent observations of hydrology that were not used as inputs to BioFRE. The system-wide hydrologic observations that we used for comparisons are reported in Harvey and Choi (2022a) and Harvey and Choi (2022b), along with the more detailed BioFRE simulations for DPM that we describe in the present document and in Choi and Harvey (2023).

For the present study we used the BioFRE to simulate overland flow in the DPM over the ten-year period of experimental high-flow releases. We used the simulations to assess the relative importance of changing water supply caused by year-to-year precipitation variability and high-flow releases through the S-152 culvert versus changes in wetland vegetation and related landscape factors that control outcomes such as water depth, a master variable that determines Everglades health. We focused on understanding how water velocity was affected, because it has immediate effects on sediment entrainment and downstream delivery of sediment-associated contaminants, such as phosphorus. The product of depth and velocity is the unit-width discharge, which is also an important outcome that controls the rate of downstream transport of sediment and phosphorus.

Three major groups of BioFRE simulations were undertaken for DPM: (1) DPM simulations of observed high flow, (2) DPM scenario simulations that held total discharge constant and varied landscape conditions, and (3) DPM sensitivity analyses that examined the relative importance of factors driving degradation or restoration of wetland functions.

Simulation group 1 tracked observed high-flow conditions in the DPM based on the observations between 2013 and 2020, during which the DPM overland flow conditions varied as a function of changing S-152 culvert discharges as well as changes in wetland landscape resulting from AMI implementation. Simulation group 2 isolated the effects of AMI implementation by holding culvert discharge constant as well as simulating several full restoration scenarios that decreased ridge proportion to historical, pre-drainage levels and altered slough vegetation type. Simulation group 3 quantified sensitivity of the simulation outputs with respect to inputs, as a means to understand the relative importance of the different drivers of Everglades degradation and restoration. Data sources and estimation procedures to create inputs for the three DPM simulation groups are presented below.

Group 1 simulations of observed high-flow releases.—DPM overland flow was simulated for eight steady high flow releases between 2013 and 2021. Simulation group 1 contrasted pre- and post-implementation of the AMI project. Four measurement periods occurred prior to AMI (pre-AMI periods 1, 2, 3, 4) and four occurred after (post-AMI periods 6, 7, 8, 10). Each measurement period also differed in total culvert discharge, ranging from 5.30 to 11.04 m³/s as well as in water surface slopes (table 8).

Eight of the eleven high-flow releases had periods of a month or more when the total culvert discharge was judged to be steady enough for our analysis. To calculate the water slope during each period, a logarithmic function was regressed to achieve a best match to measured surface water elevations for each period as a function of distance from the S-152 culvert. At each distance, slope was computed by taking the derivative of the regressed logarithmic function from table 8. For situations where surface water elevations were not available at an important target site during a measurement period, we used a regression between S-152 discharge and the offset in water elevation between RS1D and that target site, measured at times when both were available, to estimate the missing surface water elevations. The regression parameters used were regression slope coefficients 0.00013, 0.0001198, and 0.00018, respectively, and intercepts 0.08998, 0.02222, and 0.07577, respectively, to predict water surface elevations at Z51_USGS, RS2, and UB2.

Effective upscaled roughness (n_eff) for simulation group 1 was calculated based on observations of vegetation architecture, microtopography, slough water depth, and related landscape parameters, including fractal dimension of ridge-slough edges (f_d) and directional connectivity index (DCI) (Choi and Harvey, 2023; Larsen and others, 2017a). Effective roughness was affected beginning in late 2016 by implementation of the AMI project that opened two historical sloughs in the DPM. The pre-AMI ridge proportion (0.835) was relatively high and was lowered to a post-AMI ridge proportion of 0.730 by the AMI project’s removal of sawgrass from several historical sloughs. To account for the slough openings in calculation of effective upscaled roughness, we adjusted the pre-AMI binarized ridge and slough map by reassigning the opened historical sloughs as sloughs rather than ridges, which adjusted the ridge proportion (p) in the southern DPM flow way. Other needed inputs, such as fractal dimension of ridge-slough edges (f_d) and DCI were re-computed using the correlations between p and other metrics established previously by Harvey and Choi (2022a). Those observations are used to create inputs using equations of Larsen and others (2017a) to recompute effective roughness, as documented in Harvey and Choi (2022a).

Table 9 reports the measurements used to calculate effective upscaled roughness for the eight selected measurement periods of steady high flow in the DPM. Effective roughness was calculated from microtopographic height, ridge proportion, directional connectivity, fractal dimension, and ridge and slough vegetation roughness using the equation outlined in Harvey and Choi (2022a). The ridge-slough microtopographic height difference was held constant for all simulations, but the landscape pattern parameters (p, DCI, f_d,) were adjusted between pre- and post-AMI periods. Ridge proportion (p) reflects how the relatively high pre-AMI ridge proportion during the first four periods was lowered somewhat by the AMI project opening of two sloughs for the latter four post-AMI periods. Values of fractal dimension (f_d) and DCI were adjusted to account for AMI effects using the correlations of those parameters with p (Harvey and Choi, 2022a).

In summary, 24 sets of input data were prepared for simulation group 1. The simulations predicted overland flow through the DPM for eight flow periods by accounting for effects of various combinations of culvert discharge, landscape conditions such as slough area that was affected by the AMI project), and particular type of vegetation present in the sloughs (for example, dense spikerush with periphyton [DSP], dense spikerush [DS], and sparse spikerush [SS]). Each period had a different average culvert discharge, Q, through the S-152 culvert that spread radially with a spreading angle(θ) that was presumed to be constant across simulations. Water slope for the simulations was calculated from observed surface-water elevations during each period being simulated. Effective upscaled roughness for the simulations was estimated for each period based on observed vegetation-type landscape conditions. The measured average water depth in sloughs in the simulations for each period allowed us to select the appropriate depth-dependent estimates of ridge and slough vegetation roughness. The type of slough vegetation can vary over time, and we tested the influence of slough vegetation type by repeating each period’s simulation three times to reflect how slough vegetation influenced the effective roughness (table 9). The input parameters for all 24 group 1 simulations are provided in table 10.

Group 2 simulations of DPM high-flow scenarios.—The S-152 culvert discharge was held to a constant value of 7.38 m³/s (average of eight high-flow periods) for the purpose of investigating vegetation and landscape effects on overland flow. In contrast to group 1, where measured surface-water slopes were used as inputs for simulations, water slope needed to be estimated for the unobserved conditions for group 2. To predict water slope for pre- and post-AMI scenarios, we used the average of the four water slopes from table 8 as estimates for pre-AMI scenarios and the average of the four post-AMI cases from table 8 as estimates for post-AMI scenarios (table 11). For fully restored slough area scenarios, we estimated water slope by extrapolating from a linear regression that related measured water slopes at RS1D, RS2, and UB2 for eight periods (table 8) against the appropriate pre- and post-AMI ridge proportions (0.835 and 0.730, respectively). The resulting regression equation was used to estimate surface water slopes at RS1D, RS2 and UB2 based on best-fit regression slope coefficients 6.2434 E−05, 1.9799 E−05, and 1.0752 E−05, respectively, and intercepts 1.3376 E−04, 4.2418 E−05, and 2.3035 E−05, respectively. When the predicted slopes were multiplied by distances(r) the resulting water slope coefficients were 0.0837 for the flour pre-AMI scenarios and 0.0807 for the four post-AMI scenarios.

We compared group 2 scenarios to assess outcomes for several important variables affecting Everglades wetland health, that is, water depth and sediment transport conditions. For depth, it is well understood that water depth must remain high enough to preserve diverse ridge and slough wetland communities (for example, see Choi and Harvey, 2017). For sediment transport, we assessed outcomes to determine how the area of sediment entrainment from both ridges and sloughs changes, as well as how acres of sediment redistribution from sloughs to ridges (that is, sediment entrainment in sloughs and sediment deposition on ridges) changed. The thresholds for both sediment transport metrics were assessed with velocity criteria based on the experimental observations and modeling results for sediment transport in the Everglades by Harvey and others (2011) and Larsen and others (2009b). Slough velocities of 1.5–3.0 cm/s delimit a typical range wherein sediments are “redistributed” from slough to ridge, meaning that sediments are mobilized in sloughs and transported until they are deposited on ridges. Fine particulate sediments are mobilized by entrainment from their points of attachment to surfaces of floating vegetation and associated periphyton mats, or attachment to stems of rooted vegetation and associated epiphyton coatings (Harvey et al., 2011). At the higher end of the velocity range there also may be entrainment of flocculent organic sediments from the bed. Above a velocity of 3.0 cm/s, sediments are typically entrained in both slough and ridge, and below a velocity of 1.5 cm/s, sediments are not entrained in either ridge or slough (Larsen and others, 2009b).

The last of the group 2 scenarios anticipated how a fully restored slough area, which lowers effective roughness, interacts with total culvert discharge. We raised total culvert discharge by almost a factor of five, from the average observed value of 7.28 m³/s to 35 m³/s, in order to anticipate how much the culvert discharge would have to increase in order to sustain a slough water depth greater than 50 cm in a fully restored DPM where ridge proportion was decreased from its present-day value of 0.735 to 0.421. However, to complete that simulation, we needed an estimate for the water surface slope coefficient, which was uncertain. Since both our observations (table 8) and our governing equation tell us that surface-water slope depends on total culvert discharge, we estimated the water slope coefficient based on the extrapolation of a regression between measured slope coefficient and culvert discharges using all of the observed data for all eleven of the high flow periods (table 8). The estimated Slope Coefficient was 0.155 based on the estimation equation for Q (35 m³/s) which was 0.072+0.003×(35 m³/s−7.38 m³/s ) = 0.155, where 0.072 is the Slope Coefficient that we had estimated for fully restored slough area (Q=7.38 m³/s), and 0.003 is the proportionality between discharge through the S-152 culvert structure (Q) and 11 measured slope coefficients (table 8). Table 11 summarizes the measurements that were used to estimate inputs for the BioFRE base simulations at DPM.

Estimating the effective upscaled roughness for group 2 scenario simulation required an estimate of water depth in sloughs. For the pre- and post-AMI base case simulations, we used the average measured water depths for pre- and post-AMI bases cases as the estimate for water depth. For the unobserved bases cases for fully restored conditions, there are no estimates of water depths to use. Therefore, we began by using a guess for water depth that was needed to establish an initial estimate of effective roughness. Those estimates were both updated using an iteration approach that after 25 steps converged on a final solution for water depth and effective roughness. The resulting outcomes for a final simulated water depth had an iteration uncertainty of between 0.2 and 3.2 cm, which was judged to be satisfactory for our purposes in performing the base-case simulations.

Table 12 provides the input parameters for the group 2 scenario simulations of DPM high flow. Simulations 1–7 held culvert discharge constant to illustrate effects of the AMI project compared with full restoration of the slough area. Simulation 8 increased culvert discharge to answer the question, “how much more flow would be needed to avoid shallow water depths and overly dry conditions in a fully restored Everglades landscape?”

Group 3 sensitivity analysis of DPM simulations.—We computed sensitivity for DPM high flow simulations with a purpose to assess the relative importance of controlling factors, some management related and some related to landscape conditions, that together influence the pathways and speed of degradation and restoration processes in the Everglades. We used a local sensitivity approach in which the input variables Q and θ represented management related controls on overland flow, as opposed to the landscape variables ridge proportion, directional connectivity, and fractal dimension, microtopographic height, and vegetation roughness parameters n_s and n_r. The effective roughness, n_eff incorporates all controls imposed by landscape, microtopography, and vegetation conditions. Inputs were perturbed one at a time, with perturbations being small, usually an amount equal to 20 percent of the observed range of the input variable (table 12). Ridge proportion (p) was perturbed by an amount equal to 20 percent of its measured range in 23 Everglades landscapes; directional connectivity (DCI) and fractal dimension (f_d) were perturbed using a regression equation that predicted the change in those landscape metrics based on the 20 percent change in the ridge proportion previously described. Microtopography (Z_p) was perturbed by 20 percent of the range of area-averaged microtopographic height at 57 locations (Harvey and Choi, 2022b). For the base case using a culvert discharge of 7.38 m³/s, the discharge was perturbed by 20 percent of the range of selected 8 periods for restoring direction. For the base case using a discharge of 35 m³/s, the discharge was perturbed by 20 percent of the difference between 35 and 7.38 m³/s. For both base cases, the spreading angle was perturbed by 20 percent of the base case’s spreading angle (90 degrees).

For each perturbation we re-ran the simulation to quantify the relative change in several simulation outputs and computed dimensionless sensitivity of resulting outcomes as:

where “objective” is the value of a simulation output of interest. Here, the objectives of interest were the simulated water depth, sediment entrainment area, and sediment redistribution area. Computations of dimensionless sensitivity were compared and ranked by importance.

A local sensitivity analysis approach is valid for situations that are similar to the tested base case. For that reason, we performed sensitivity testing for two relevant endmember scenarios: the first scenario representing a poorly functioning but potentially restorable landscape. To do this, we used the pre-AMI simulation with dense spikerush and periphyton in sloughs (see scenario 1 from table 12) as the base case, and we refer to it as the “restoring” test because we perturbed inputs in a direction to improve wetland function. The second scenario represented a well-functioning, fully restored landscape that can be degraded by lowering of water depth and other factors. For that, we used the fully restored simulation (see scenario 8 from table 12) with dense spikerush and periphyton in sloughs as the base case, and we refer to it as the “degrading” test because we perturbed inputs in a direction that would degrade wetland function (table 13).

Results and Discussion

Simulating high-flow releases and overland flow through the DPM makes use of BioFRE equations that are appropriate for radial flow to estimate how point-source releases through culverts spread through downstream wetlands. The simulations can help anticipate future restoration outcomes for water availability and ecosystem health of the Everglades. Before discussing those uses, we evaluate the performance of BioFRE in simulating observed hydrologic conditions.

Introduction

A loss of high flows through a patterned landscape with microtopographic variation and ridge and slough vegetation communities is thought to degrade ecosystem functions, as well as diminish water storage during droughts and diminish the capacity for wetlands to convey flood waters without drowning tree islands or putting urban communities at risk. Modeling supports a hypothesis that high flows entrain organic sediments in sloughs and deposit them on ridges, a process described as sediment redistribution. The process has been investigated in small laboratory flumes and flumes installed in the Everglades (Larsen and others, 2009a; Harvey and others, 2011) and in tracer experiments in the field (Regier and others, 2018), however the flume work was primarily aimed at understanding entrainment of suspended sediment and phosphorus and not deposition. Furthermore, although the tracer work indicated that sediment redistribution is operative, the results are qualitative. Sediment transport is thought to be an important driver of phosphorus transport because much of it is associated with sediment in the Everglades. Confirmation of the magnitude of transfer of sediment and phosphorus from sloughs to ridges is needed at larger spatial scales within ridge and slough ecosystems.

DPM is a 10-square-kilometer area of Everglades wetlands that have been isolated by levees for nearly sixty years. The recent completion of the S-152 gated structures on the L67-A levee had permitted the release of large quantities of water (hundreds of cubic feet per second) that spread across the DPM wetlands and raise velocities substantially, by about two to ten times within a kilometer of the culverts. Here we report on measurements of suspended sediment and total phosphorus in sloughs and ridges and implement a simple flow model to assess the importance of sediment redistribution.

Methods

Methods of data collection and analysis are described in the “Suspended Sediment Sampling” and “Water Quality and Physiochemical Monitoring” sections. Here, we describe the modeling methods used to assess transport and redistribution of suspended sediment.

Model estimation of sediment redistribution.—We used the One-Dimensional Transport with Inflow and Storage (OTIS) model (Runkel, 1998) to simulate the transport and reaction of suspended sediment and total phosphorus during flow releases in 2015. Our modeling approach was designed to distinguish the fate of suspended sediment and total phosphorus flowing through a slough and ridge landscape. The input and downstream movement of water and suspended constituents were simulated before and during the operation of the S-152 culverts, which deliver high flow into the DPM. The OTIS simulation accounts for flow through a main flow path (slough) that exchanges water with secondary flow paths moving through ridges. It also accounts for removal of constituents by settling to the bed, by getting captured on vegetation stems, or by reactive transformation.

The OTIS governing equation is:

where the left side of the equation represents the changing concentration, C, in the primary flow paths through sloughs and the right-hand side of the equation represents effects of several processes including advection, longitudinal dispersion, lateral mixing, hydrologic exchange with secondary flow paths through ridges, and removal of constituents from flowing water, either by settling, capture on vegetation stems, or reactive transformation. The associated governing equation for concentration in the secondary flow paths through ridges is: where the variables of the above equations are defined as:

The application of the OTIS model for simulating suspended sediment concentration (SSC) and total phosphorus (TP) transport through sloughs, including exchange with ridges and settling and particle capture on vegetations stems, can benefit by first simulating a conservative solute tracer. The conservative tracer helps constrain the transport time and rate of water exchange between the main flow path through a slough and secondary flow paths through ridges. For the conservative solute tracer at DPM we simulated total dissolved solids as detected by measuring specific conductivity (SC). Because the SC of surface water entering the DPM through the S-152 culverts differed from the SC in the DPM wetland, and because total dissolved solids are transported without reaction, the SC signal could be used to quantify conservative transport processes, specifically, advection and longitudinal dispersion of flow through sloughs, and rate of water exchange between the sloughs and secondary flow paths across ridges. With those parameters constrained it was easier to model total phosphorus (TP) which undergoes the same conservative transport processes but also is affected by reactive processes. Specifically, TP concentrations are also affected by settling of phosphorus associated sediment particles and capture of particles on vegetation stems.

Three OTIS models were solved, one for SC and the others for SSC and TP, to determine the “best-fit” estimates of the parameters. The flow data, specific conductivity, suspended sediment, and total phosphorus data used in the calculations are available in Choi and others (2016), Harvey and others (2022a) and the model inputs and outputs are detailed in Choi and Harvey (2023).

The following sections describe two steps to solve OTIS models to arrive at final answers. In the first step, an OTIS model for SC was solved to match the measured time series of SC in the slough which produced the conservative transport parameters for a slough flow path interacting with more slowly moving flow paths though ridges. In the second step, the conservative transport parameters from the SC simulation were used as fixed parameters in the reactive transport simulations, from which the “best-fit” reactive parameters for SSC and TP were determined for slough flow paths interacting with ridge flow paths. The two steps are explained in detail below:

Step 1. Estimating conservative transport of constituents in ridges and sloughs.—The OTIS modeling approach helps interpret measured concentrations by distinguishing the effects of transport and reaction processes. Time-series data from the 2015 flow releases were used to prepare model inputs including water velocities, water depths, and SC, SSC, and TP to optimize parameters describing transport and reaction in the main flow zone through sloughs and secondary flow paths through the more densely vegetated ridges.

Conservative transport simulations were undertaken by modeling transport between sites Z51_USGS and RS1D in 2015. In 2015, there were two flow releases, and we simulated specific conductivity from the second release because there was a better distinction between the SC of the flow-release water and background values in the DPM. To account for slough sinuosity, the measured linear distance of the 2014 and 2015 reaches was increased by multiplying by 1.3. Discharge (Q) was estimated for each simulation as the product of measured velocity and depth at the upstream end of the reach. Lateral inflow (q_Lin) was set to zero for each simulation. After the flow release had steadied, the differences in measured velocities and depths at the upstream and downstream endpoints were used to estimate cross-sectional area (A) for each reach. We then created input files using the estimated Q and A and the adjusted distance along to run OTIS models that provided optimal estimates of the dispersion coefficient (D), the exchange rate (α) and storage cross sectional area (As) as outputs. Optimization was accomplished using the Parameter Estimation Software Tool (PEST) (White and others, 2020) in tandem with OTIS to objectively fit the measured SC data.

Step 2. Estimating SSC and TP removal rate constants in ridges and sloughs.—Reactive transport simulations were undertaken with the purpose of quantifying and comparing SSC and TP removal in sloughs and ridges. The best-fit conservative transport parameters (Q, A, D, and A_s) that had been determined from the conservative transport simulation (see section “Dissolved oxygen and temperature modeling”) were used as inputs to the reactive transport simulations. The reactive transport simulations were aimed at estimating SSC and TP mass removal rate constants (λ and λ_s ) for main flow paths (slough) and secondary flow paths (ridges), respectively. Like the conservative transport simulation, the reactive transport simulations were optimized using the PEST code. As a result, SSC and TP removal rate constants were estimated that achieved the “best fit” to measured SSC and TP concentrations.

Partitioning mass removal between slough and ridge required four simulations each for SSC and TP: (1) a simulation to quantify transport through the reach without removal (simulation A; described in section “Dissolved oxygen and temperature modeling”); (2) a simulation to quantify removal rate constants for both slough and ridge (simulation B); (3) a simulation to quantify the removal rate constant if removal only occurred in the slough (simulation C); and (4) a simulation to quantify the removal rate constant if removal only occurred in the ridge (simulation D). In summary, the four simulations were:

Mass removal fluxes of SSC and TP for slough and ridge were computed from differences between the conservative simulation in step 1 and the appropriate reactive simulations in step 2. First, the mass fluxes through the reaches were computed for each simulation from the simulated concentration curves at the downstream end of the reach multiplied by the average discharge for that period:

Mass fluxes for simulations C and D can be subtracted from simulation A to estimate removal in slough and ridge, respectively. Those estimates are biased slightly high because of assuming independence of removal in slough and ridge. However, dividing those quantities by simulation B accurately computes the fraction of SSC and TP mass removal in slough and ridge, respectively, during the period of interest. For example, (Mass flux__Sim-A – Mass flux__sim-C)/(Mass flux__Sim-A – Mass flux__sim-B) computes the fraction of mass removal that occurred in the slough. Likewise, (Mass__Sim-A – Mass__sim-D)/(Mass__Sim-A – Mass__sim-B) computes the fraction of mass removal that occurred in the ridge. Multiplying the fractions of removal by (Mass__Sim-A – Mass__sim-B) adjusts and finalizes the calculation of mass removal occurring in slough and ridge respectively during the period of interest.

Results and Discussion

We were interested in comparing mass removal that occurred during the initial pulse of a high-flow release with mass removal that occurred later during steady high-flow. The initial pulse lasted approximately 6 hours when SSC or TP concentration spiked, but concentrations tended to settle down during the steady flow period that followed, which lasted as long as the culverts were open, often 30 days or more. To accomplish that, we made mass removal computations for nine distinct periods immediately before and for several months after the flow releases in 2014 and 2015. The periods were selected to distinguish the different flow conditions that were well described by our measurements of SSC and TP concentrations.

The modeling characterized sediment and total phosphorus areal deposition in sloughs and in ridges caused by the 2015 flow release from the S-152 culverts. The areal deposition numbers characterize the average deposition in both sloughs and ridges between sites Z51-USGS and RS1-D. The instantaneous rate of areal deposition is compared for the initial pulse of high flow and a longer period of steady high flow. We also compared a 3-month set of conditions, where it was assumed that only 6-hr pulsed flows occurred every 14 days and allowed time for regrowth of transportable biomass (for example, Harvey and others, 2011), which was compared with deposition associated with a 3-month steady flow (tables 16 and 17).

Consistent with the sediment redistribution hypothesis, rates of sediment deposition on ridges were 3 to 6 times greater on ridges compared to sloughs during high-flow releases at DPM. Also, the short-term pulsed rates of deposition were higher than the longer-term steady-state rates in the sloughs, especially for the first of the single-day pulses (the second pulse was conducted a day later). The generally higher pulsed rates of sediment deposition reflect the time needed to build up the standing stock of the easily entrainable fine particulates associated with periphyton in the water column (Harvey and others, 2011), which comprises the major source of suspended sediment that settles to the bed (Noe and others, 2010).

Phosphorus deposition followed similar patterns to sediment deposition, being higher on ridges compared to sloughs. Ridge-slough differences during pulses were pronounced, with ridge deposition rates typically a factor of 50 greater than sloughs. The ridge-slough difference during the longer-term period of steady-state high flow was much less, by just a factor of 2, which was more similar to ridge-slough differences in sediment deposition. The steady-state rates of phosphorus deposition were more moderate, falling between the very low values in sloughs during pulses and the very high values for ridges during pulses. These findings suggest a higher efficiency of ridge capture of phosphorus during pulsed flows, however those pulsed-flow behaviors are unsustainable as reflected in the lower deposition rates of the second pulse compared to the first, as well as in the much lower steady-state rates. Like sediment deposition, the lower rates of phosphorus deposition over longer periods are presumed to be the result of the time needed for replenishing the standing stock of easily entrainable phosphorus in the water column (Noe and Childers, 2007).

Most of the rates of phosphorus deposition during high flows in these DPM studies are higher than published values for oligotrophic sloughs by approximately 10 to 20 times for the longer period of steady high flow (Noe and Childers, 2007), and more for pulsed flow, notably on ridges (however, slough deposition rates during pulsed high flow were comparable to the background low-flow values of Noe and Childers [2007]). Reasons for the high- and low-flow difference include the much higher loading of phosphorus that occurs with high flow, typically an order of magnitude higher for steady high flow, assuming that concentration changes little and the increase in transported phosphorus load is the same as the flow increase caused by the culverts. An order of magnitude higher phosphorus loading during steady-state high flows therefore approximately explains the order of magnitude higher phosphorus deposition rates for steady high flow (table 17). The spike in concentrations sometimes observed during pulsed flow indicates even higher load increases for a brief period of 6 hours or so which resulted in deposition rates on ridges as much as 60 times than literature values for low-flow conditions in sloughs. Those high values are for the edge of ridges and not central ridges, which agrees with modeling by Larsen and others (2009b) and Larsen and Harvey (2011) indicating that sediment redistribution is the result of entrainment in sloughs and deposition primarily on ridge edges. Another reason for higher phosphorus deposition values compared to low-flow literature values is the radial flow pattern that causes entrainment upstream and deposition downstream because of the logarithmic decrease in velocity as a function of distance. Our selected test reach encompassed conditions of entrainment in both slough and ridge at the upstream end and conditions near the threshold for entrainment in sloughs but deposition in sloughs at the downstream end. Thus, there are several reasons why our observations of higher phosphorus deposition rates are reasonable for DPM high-flow conditions with radial spreading compared to oligotrophic sloughs under low-flow conditions.

Introduction

A key aspect of Everglades restoration is the removal of levees to promote sheet flow through wetlands to better mimic the pre-drainage system. The DPM reintroduces sheet flow to an area of the central Everglades that has been isolated for nearly sixty years. An important objective of DPM is to test how levee gap design and the degree of backfill (no backfill, partially backfilled, and completely backfilled) influences sheet flow and conveyance of suspended sediment, total phosphorus, and total particulate phosphorus to downstream areas. In 2012, a 914-m section of the L67-C levee was bulldozed to grade with the surrounding wetlands, and the excess material was used to construct contrasting levels of backfill in three side-by-side treatments within the exposed L67-C canal. We found that backfill treatments differed in how they conveyed constituents, that is whether the different treatments passively conveyed constituents or whether treatments acted as a sink by retaining constituents as a result of deposition or biological uptake, or whether treatments acted as a source of constituents by deposition by entrainment from the canal bottom as water passed through them. We found that the treatment with no backfill (herein referred to as the no-backfill treatment) passively conveyed water and constituents, the partially backfilled treatment (herein, partial-backfill treatment) exhibited tendencies for both sources and sinks of sediment and phosphorus, and the completely backfilled treatment (herein, complete-backfill treatment) tended to retain total phosphorus and total particulate phosphorus. Low-flow and high-flow measurements at the L67-C gap also revealed an unintended effect: the unblocked portions of the L67-C canal adjacent to the gap became a major source of water flowing through the levee gap, which has the potential to transport undesirable levels of phosphorus to areas farther downstream.

Site Information and Methods

The DPM is in a small area in of the central Everglades where the wetlands have been isolated between the L67-A and L67-C levees for nearly sixty years, referred to as the pocket. The only sheet flow that occurred through the wetlands was caused by groundwater exchanges between the two levees and by local precipitation. The recent completion of the 914-m wide L67-C levee gap in August of 2012 permitted sheet flow from DPM to begin flowing directly through the gap and into WCA 3B. Flow through the gap was minimal at first because there was no new source of water to the DPM to sustain flow through the gap. Soon after, in November 2013, the S-152 gated structures on the L67-A levee were completed, which permitted the release of substantial quantities of water from WCA 3A into the DPM, spreading across the wetlands, and passing through the L67-C levee gap into WCA 3B.

The DPM high-flow experiments are being conducted in an area of 3×3 km that lies between two parallel levee-canal systems. At the upstream end of the DPM, a bank of 10 gated culverts was completed in 2013 to convey water from WCA 3A through the L-67A levee at rates up to 650 ft³/s in a radial pattern that spreads across the pocket and eventually passes through an approximately 914-meter gap in the L-67C levee that was bulldozed to the wetland grade to permit flow from the pocket into WCA 3B. The L-67C levee was flattened in 2012, at which time the canal portion directly adjacent to the levee degrade was backfilled to varying degrees in three 300-m treatment areas built to test the conveyance of water, sediment, and phosphorus across the 900-m levee gap into WCA 3B.

The three canal backfill treatments include an undisturbed no-backfill section (CB1), a partial-backfill section that was backfilled with rock to a depth of 30 cm on top of the pre-existing canal bottom, and a complete backfill section (CB3) that had been backfilled with rock to a grade resembling a 1- to 2-m deep slough in the surrounding wetlands (fig. 13). The height of the degraded levee surface is uneven, and it has two deep openings, one at the south end of CB3 and the other at the north end of CB1, where water flowing into WCA 3B is concentrated. There is also a relatively long portion of CB3 where the levee was not flattened sufficiently so water rarely overflows.

Here we report on a mass balance analysis of water and constituents entering and leaving the backfill treatments that was conducted during low-flow conditions in November 2015 and during high-flow conditions in November 2016. The purpose was to quantify the source of water and constituents entering the levee gap and to assess whether materials were passively conveyed through the gap, often referred to as “conservatively transported,” whether constituents settle or are removed by biological reactions (referred to as a “sink”), or whether constituents are entrained or released by biological or geochemical reactions (referred to as a “source”) during transport through the gap. Mass balances were computed in three more years (2019, 2020, and 2022); however, only water balances could be computed because suspended sediment and phosphorus concentrations were not measured in those years.

Field Measurements.—All field measurement locations are shown either in figure 1, which provides an overview, or in figure 14, which is a detailed inset of the L67-C canal gap area. Low-flow measurements were made on November 14, 2015, when the S-152 structure gates were closed and not releasing water into the DPM area. Inflows from the unblocked extensions of the L67-C canal just outside of the gap (CCN and CCS) were measured using two continuously recording SonTek Argonaut XRs. Sheet flow measurements during the low-flow campaign were only made at the three locations entering the levee gap (UB1, UB2, and UB3 sites) and three locations exiting (DB1, DB2, DB3 sites) the gap. The measurements were made using six SonTek ADVs, which were positioned at three fixed stations on the upstream and downstream sides of the levee gap.

High-flow measurements were made after flow had steadied following the opening of the S-152 structure gates in the fall of 2016, 2019, 2020, and 2022. Rather than just six measurements of sheet flow, many more measurements of sheet flow were made during the high-flow campaigns. During each flow release, approximately 87 measurements of sheet flow were made, which required roving teams to use hand-held SonTek Flowtrackers.

In 2016, the high-flow measurements were made on several days over two weeks that later were adjusted, as explained below, so that all flows were estimated on a single day, November 29, 2016. Flow adjustments were necessary to meet the assumption of the steady-state (non-changing) flow analysis. Sheet flow measurements entering and leaving the canal were made on November 26, 2016, and suspended sediment and phosphorus concentration measurements were made by SFWMD on November 14, 2016, and November 15, 2016, within the canal and downstream of the canal. Canal velocity cross sections were measured on November 17, 2016, using a roving SonTek acoustic Doppler current profiler (ADCP). However, because of short-term blocking of the S-152 culverts on that day, the flows were lower than the days when the other mass balance measurements were being made. The roving ADCP flows on November 17, 2016, were corrected to represent conditions on November 29, 2016, by using a relationship between S-152 culvert structure discharge and continuous measurements from the fixed ADCP station at cross-section CCN (fig. 14).

In 2019, the high-flow measurements were made on several days during a one-week period that were later adjusted, as explained below, so that all flows were estimated on a single day. Again, this was necessary to meet the assumption of the steady-state (non-changing) flow analysis. Sheet flow measurements entering and leaving the canal were made on October 15–16, 2019, discretely by SonTek FlowTracker. During October 8–10, 2019, canal velocity cross sections were measured discretely, using a roving SonTek ADCP and continuously using a SonTek Argonaut XR at the northern gap (CCN) and a SonTek IQ plus at the two downstream transects (CB2_S to CB3_N and CCS) (fig. 14). We accounted for the temporal changes in flow during the measurement period at the CCN location by replacing the water flux measured on October 8, 2019, with a water flux on October 16, 2019, estimated using a linear relationship between the flux at CCN and the discharge at the S-152 culvert. Furthermore, we used continuous discharge measurements made by a SonTek IQ plus on October 15–16, 2019, (the average of two days), during which times the S-152 culvert discharge was steady, in place of the discrete measurement on October 10, 2019, by the ADCP to better account for the temporal change at two other canal cross sections—CB2_S to CB3_N and CCS (fig. 14).

In 2020, the high-flow measurements were once again made on several days over a one-week period, with some flows being adjusted as explained below to meet the assumptions of the steady-state (non-changing) flow analysis. Sheet flow measurements around the canal were made on September 8–9, 2020, discretely by a SonTek FlowTracker. During September 1–3, 2020, canal velocity cross sections were measured discretely using a roving SonTek ADCP and measured continuously using a SonTek Argonaut XR at the northern gap (CCN) and Sontek IQ plus at the two downstream transects (CB2_S to CB3_N and CCS). We accounted for the temporal changes in flow during the measurement period at the CCN location by replacing the water flux measured on September 2, 2020, with a water flux on September 8, 2020, estimated using a linear relationship between the flux at CCN and the discharge at the S-152 culvert.

In 2022, the high-flow measurements were once again made on several days over a one-week period. Sheet flow measurements around the canal were made on January 25, 2022, and January 26, 2022, discretely by Sontek FlowTracker. Canal velocity cross sections were measured both discretely on January 24, 2022, to January 27, 2022, using a roving SONTEK ADCP. The temporal change in discharge at the S-152 culvert between the periods of sheet flow and cross section measurements was very small, therefore we did not apply any correction on the cross-section measurements.

Mass Balance Approach.—Mass balances were computed for a levee gap control volume, enclosing the 914-m opening in the levee and canal area between DPM and WCA 3B. Mass-balance computations assume steady-state conditions, that is, that no changes in storage occur in the control volume for the time period represented by the input data. Computations were made for water fluxes, suspended sediment fluxes, total phosphorus fluxes, and total particulate phosphorus fluxes entering and leaving the levee gap and backfill treatment control volumes. To compute flows entering and leaving, the local measurements of sheet flow were resolved in the direction normal to the gap and upscaled to account for the flux across an appropriate portion of the control volumes outer boundary. Likewise, canal fluxes were resolved in a direction parallel to the gap and upscaled to account for the flux across an appropriate endpoint of the canal gap.

The much higher number of measurements made during the high-flow campaigns allowed mass balances to be computed for each of the three backfill treatments, CB1 (no backfill), CB2 (partial backfill), and CB3 (complete backfill) (fig. 13). Furthermore, because of several localized areas of channelized flow across the degraded levee, it was beneficial to split the backfill treatments in half, which created a north and south section for each of the three backfill treatment, that is, representing six individual sub-sections of the entire levee gap (fig. 14).

To compute the mass balances for high-flow releases, the 87 sheet-flow measurements were therefore resolved into six fluxes entering a sub-section from the upstream (DPM) side and six fluxes leaving a sub-section on the downstream side to enter WCA 3B. All of the flux measurements were available to compute the balances except for two canal flux measurements, that is, the flux between CB1-south and CB2-north and the flux between CB2-south from CB3-north (fig. 14). Estimates of canal fluxes at those two locations were computed by difference, that is, the flux was computed that was needed to balance water mass exactly in the northern sub-section of each backfill treatment. That approach left a small water balance residual in the southern sub-sections of each backfill treatment which was deemed to be acceptable. Figure 14 names the fluxes that were used as inputs for the mass balance calculations.

Mass Balance Calculations.—The steady-state mass balance equation is:

where inflows include sheet flow from DPM wetlands directly upstream of the levee gap plus lateral canal flow from the unblocked extensions of the L67-C canal, and outflows are sheet flow to downstream wetlands in WCA 3B.The variable $U_{\epsilon }$ is the cumulative uncertainty associated with the mass balance measurements estimated by first-order uncertainty principles.

Measurement uncertainty was estimated separately for water, suspended sediment, total phosphorus, and total particulate phosphorus fluxes by assuming a priori a given value of uncertainty for water fluxes and a given uncertainty for constituent fluxes. We followed the general guidance of Winter (1981) and our experience in measuring Everglades fluxes to assign levels of uncertainty for the DPM levee gap study as 15 percent uncertainty for measurements of water fluxes and as 30 percent uncertainty for measurements of constituent fluxes.

Individual measurement errors propagate and are accumulated in the calculated mass balance uncertainty ($U_{\epsilon }$) shown in equation 36. To estimate the accumulated uncertainty, we assumed that individual measurement errors were independent of one another, which is referred to as the first-order uncertainty assumption. The cumulative uncertainty that results from summing measurements each with their own uncertainty is estimated as:

where m₁–m_n are the measurements of water flux, suspended sediment flux, total phosphorus flux, or total particulate phosphorus flux used in the mass balance, and u₁–u_n are the estimated uncertainties for each measurement, that is, 15 percent uncertainty for measurements of water fluxes and a 30 percent uncertainty for measurements of constituent fluxes.

Estimating sources or sinks of constituents within the levee gap.—If a calculated net flux for water or a constituent computed using equation 35 was larger than the associated computation of uncertainty using equation 36, then it was interpreted as evidence that the levee gap is a source (if the net flux is negative) or a sink (if the net flux is positive) for water or that constituent. Such results provide indicators of source-sink dynamics although they cannot identify what processes are responsible. Interpretations from mass balance analysis can be strengthened using additional measurements that independently support conclusions about the processes involved, for example, direct measurements of sediment settling or entrainment from the canal bed, and so forth.

To isolate individual effects of the canal-backfill treatment, the net-flux calculations were computed for six canal-backfill treatments (CB1-N, CB1-S, CB2-N, CB2-S, CB3-N, CB3-S) that subdivided the backfill treatment areas into north and south sections (fig. 15).

Results and Discussion

During the low-flow measurement period in November 2015, when the S-152 structure gates were closed, we observed 4.8 m³/s flow through the levee gap that had collected in the pocket because of groundwater seepage beneath the L-67A levee as well as local precipitation. At low flow, two-thirds (67 percent) of the water that entered the levee gap came from the unblocked L67-C canal and one third (33 percent) entered as wetland sheet flow from directly upstream (table 17).

During the high-flow measurement periods the total flow of water through the levee gap increased because of the added discharge from the S-152 gated culverts (table 18). Compared to the low flow discharge of 4.8 m³/s into the gap, total discharge through the gap increased by 25 percent to 6.1 m³/s in November 2016, again by 25 percent to 6.0 m³/s in October 2019, and by 102 percent to 9.6 m³/s in September 2020 (table 18).

Confounded effects of water sources, uneven flows, and backfill treatments.—The unblocked L-67-C canal accounted for slightly more than two thirds of water fluxes entering levee gap during 2015 low-flow and during the 2016 and 2019 high-flow measurement periods (table 19). During the 2020 high-flow period, when S-152 culvert discharge was higher, the unblocked L67-C canal accounted for a slightly smaller percentage of total flow, 62 percent (table 19).

Regardless of flow condition, the flow from the unblocked L67-C dominated water flow into WCA 3B (ranging from 61 percent to 70 percent) compared with smaller sources from wetland sheet flow (ranging from 30 percent to 39 percent) (table 19). Therefore, the unblocked L67-C canal interferes with the desired enhancement of sheet flow to improve direct wetland-to-wetland connectivity between the Water Conservation Areas in the central Everglades (Larsen and others, 2017b).

The source of DPM constituents entering the 900-m levee gap was dominated by delivery from the adjacent unblocked L67-C canal rather than from sheet flow (table 19). Approximately two-thirds of the water, suspended sediment, and total phosphorus with more than half of the total particulate phosphorus entering the levee gap from the L67-C canal. There was little change from contributions during low-flow background conditions versus during high-flow releases, with 67 percent versus 69 percent contributed for water, 67 percent versus 73 percent contributed for suspended sediment, and 63 percent versus 65 percent contributed for TP, respectively, during low flow versus high flow in terms of percent contributions from the canal. Accordingly, wetland sheet flow contributed less than half of the water and constituents entering the levee gap. For example, in the 2016 flow release, the contributions from the canal accounted for 31 percent, 27 percent, 35 percent, and 41 percent of the water, sediment, total phosphorus, and total particulate phosphorus that entered the levee gap (table 19).

Furthermore, the canal water flow from the unblocked L67-C canal enters the gap unevenly from the north end and the south end, with 96 percent entering from the north side of the gap (table 19). The dominance of the north side extension of the L67-C canal is explained by the easterly trending direction of DPM high-flow releases which, although spreading radially from the S-152 gated structures, trend eastward with a center of mass that reached the L67-C canal at a location several kilometers north of the levee gap. Therefore, sheet flow entering the levee gap might be somewhat higher than measured at DPM because of an unplanned consequence of the DPM design.

Drawing conclusions about canal backfill functions must confront challenges involving the uneven source contributions and flow distribution in the levee gap. A first challenge was that the unblocked L67-C canal, rather than wetland sheet flow, is the major supplier of water and constituents (approximately 94 percent) in the no-backfill treatment, the most northerly section of the levee gap closest to unblocked canal to the northeast that receives most of the flow. The other backfill sections received less canal flow, more than half in the partial-backfill treatment, and approximately a quarter of the flow entering the complete-backfill treatment (table 20).

A second challenge was that the magnitude of water and constituent fluxes passing through the levee gap into WCA 3B were much greater in the northernmost no-backfill treatment compared with the other treatments (table 21). Approximately sixty percent of the water and constituents entering WCA 3B passes through the no-backfill treatment on the northern side of the gap, one quarter of the water and constituents passes through the complete-backfill treatment on the southern side, and the smallest amount of water and constituents (less than one-fifth) passes through the partial-backfill treatment, which is located in the middle of the levee gap (table 21). As a result, the loading of constituents into WCA 3B was the highest by far since it was directly downstream of the no-backfill treatment rather than being evenly spread out through the entire 900-m levee gap.

The ramifications of uneven contributions of canal water in the three treatments are obscuring the effects of the different backfill functions. It would be easier to uniquely distinguish the effects of contrasting levels of canal backfill if treatments received similar mixes of wetland sheet flow and canal water, and if the total flow through each backfill treatment was closer to equal.

Evidence for canal backfill treatment sources and sinks of constituents in the levee gap.—We found that backfill treatments differed in the degree to which they passively conveyed constituents or whether the levee gap retained or was a source of constituents. Constituents might be retained because of settling of sediments or biological uptake of phosphorus. Or there might be a source of constituents because of the entrainment of sediments from the bed or a biogeochemical release of phosphorus from sediments as water passed through the gap (table 22).

The no-backfill treatment area, which had a very high flux as well as a high fraction of L67-C water passing through it, passively transported the materials entering it in both its northern and southern sections, causing significant loading of sediment, phosphorus, and particulate phosphorus into WCA 3B. (table 22, fig. 15).

The partial-backfill treatment area was a transition area where, in its northern half, suspended sediment was retained, probably by deposition, and particulate phosphorus was released from the canal bottom. In the southern section of the partial-backfill treatment, total and particulate phosphorus were retained. Sediment and phosphorus dynamics in the partial-backfill treatment could be the result of a slowdown of water as southerly flow from the unblocked L67-C canal crosses the ramp of backfill material in the partial-backfill treatment and causes deposition of sediment in its northern section and deposition of phosphorus in its southern section (table 22, fig. 15). That result also could reflect model errors that indicated a retention of water in the partial-backfill treatment (table 22), an unexpected result that is unlikely to have occurred.

Because of its southern position in the gap, the complete-backfill treatment was largely out of reach of inflows from the unblocked extension on the L67-C canal. Therefore, the complete-backfill treatment had lesser influence of canal input, receiving a quarter of its inflow from the unblocked extensions of the L67-C canal and three quarters of its inflow from wetland sheet flow (table 20). In its northern half, the complete-backfill treatment acted as a sink for phosphorus and particulate phosphorus, which is consistent with the observed presence of submerged aquatic vegetation in the complete-backfill treatment, which adds abundant surface area for microbial and algal communities with demand for phosphorus (fig. 15).

Thus, the three backfill treatments differed in how they conveyed constituents, however the observed effects are confounded to a degree with the uneven flow of water and constituents and the non-uniform canal versus sheet flow contributions. Nonetheless, it appears that the compete-backfill treatment, which best represents a backfilled canal receiving a dominant source of water from wetland sheet flow, did exhibit desirable retention of TP and TPP.

In conclusion, The DPM levee gap mass balance provided important insights and supported the following conclusions in spite of experimental design difficulties: 1) the unblocked extensions of L67-C canal significantly diminished wetland to wetland connectivity by sheet flow, which is one of the major goals of removing Everglades levees, and 2) the compete-backfill treatment dominantly received sheet flow from directly upstream and exhibited the desirable outcome of retaining TP and TPP. Thus, mass balance analysis suggests that backfilling canals can favorably influence wetland to wetland connectivity and can also perform the ecosystem function of retaining phosphorus which slows its movement to downstream areas. Our mass balance results can be combined with independently collected information such as sediment and phosphorus settling measurements and measurements of canal bed and wetland phosphorus content to provide further insights and strengthen interpretations.

Introduction

Many experimental re-introductions of high flow to low-gradient floodplains have been short term and small-scale experiments in flumes (for example, Larsen and others, 2009a; Harvey and others 2011). Recently, however, a kilometer-scale experimental facility was established to increase floodplain flows in an area of the central Everglades for several months each year. The DPM included extensive measurements to characterize short- and -longer term hydraulic, sedimentary, and ecological outcomes.

Here, we report on ecosystem outcomes of restored high flows at DPM using aquatic metabolism as an integrated measure of ecosystem response to high flows. Metabolism quantifies the balance between production by photosynthesizing organisms and consumption of organic carbon by respiring plants, animals, algae, and bacteria and through decomposition pathways. Metabolism can characterize ecological outcomes such as changing capacities for processing carbon and nutrients that may also affect geomorphic outcomes in an organic floodplain landscape. We also assessed metabolism at multiple ridge and slough sites in the DPM experimental area located in the central Everglades.

Analysis Approach

We used a Bayesian state-space model in the open-source statistical program R, streamMetabolizer, version 0.12.0, (Appling and others, 2018a) to inversely estimate aquatic metabolism in the Everglades. The model inversely estimates the parameters of a one-station oxygen mass balance equation where the change in oxygen concentration is the sum of three volumetric process rates affecting dissolved oxygen concentration in the water column:

where d[O₂]/dt is the rate of change in O₂ in milligrams of O₂ per liter per day (mg O₂/L/d) with time, t; P_t (in g O₂/m³/d) is the rate of productivity that is modeled as mean daily gross primary productivity (GPP, in g O₂/m²/d) using functional relationships with mean water depth (in meters), and light expressed in units of photosynthetic photon flux density (PPFD, in micromoles of photons per square meter per day [μmol photons/m²/d]); R_t is the rate of respiration that is modeled to estimate mean daily ecosystem respiration (ER, in g O₂/m²/d) using a functional relationship with mean water depth (in meters); D_t (in g O₂/m³/d) is the rate of O₂ gas exchange with the atmosphere that is modeled to estimate mean daily gas exchange as the product of the temperature dependent gas exchange rate constant for oxygen, $K_{O_2}$ (in 1/day) and the saturation deficit (in g O₂/m³) between water column and the equilibrium oxygen concentration. To facilitate comparisons, the model estimates K₆₀₀ (in 1/d) by normalizing by the temperature-dependent Schmidt number for O₂. Net ecosystem productivity (NEP, in g O₂/m²/d) is the sum of GPP and ER and represents the net balance accounting for the difference between oxygen (and organic carbon) production and consumption. Further explanation of the detailed solution approach to estimate GPP, ER, and K₆₀₀ are available in Appling and others (2018a, b).

We set the specifications of streamMetabolizer to incorporate both observation and process error. Observation error results from errors in measurement associated with instruments or the measurement techniques; process error occurs when some physical or chemical process is not included (for example advection and dispersion of dissolved oxygen, DO) or is mis-specified (for example, assuming ER is constant within a day, and so forth) (Holtgrieve and others, 2010; Appling and others, 2018a). According to Appling and others (2018a), the inclusion of process error in the model removes a main source of autocorrelation of DO residuals that arises from the persistence of process effects (for example clouds that dampen GPP or afternoon wind that increases gas exchange).

Bayesian hierarchical partial pooling was used to constrain the across day variability of K₆₀₀ for days with similar discharge (Appling and others 2018a) by establishing an across-day piecewise linear relationship between ln(K₆₀₀) and ln(Q). Default values of the priors were used for partial pooling because results were not sensitive to a change in these initial values. Each low and high-flow period at a site was modeled separately because of the large differences in flow conditions.

State-space metabolism models can in some cases identify erroneous solutions with high covariance of ER and K₆₀₀ when GPP is low and K₆₀₀ is high (Payn and others, 2017; Appling and others, 2018a; Appling and others, 2018b). Generally, we did not have those issues with the Everglades sites; however, to ensure good performance we monitored for when the coefficient of determination, R², of daily ER estimates with daily K₆₀₀ estimates exceeded 0.5 for the modeled period; if so those models were re-run with a smaller ${\sigma }_{{\sigma }_{K_{600}}}$ to constrain K₆₀₀ pooling and increase warmup steps (Blaszczak and others, 2019). Furthermore, if poor convergence was indicated by the Rhat (the Gelman-Rubin convergence statistic) of standard deviation of K₆₀₀ or standard deviation of process error greater than 1.1, we increased the warmup steps and re-ran the model (Appling and others, 2018a; Arroita and others, 2019). Daily outputs of GPP, ER, and K₆₀₀ also were evaluated using four criteria. Flag 1 was applied for days with poor signal strength, identified by a daily coefficient of determination, R²_det, that was both less than zero and less than the 15th percentile of all daily R²_det for that modeling period. Flag 2 was applied for biologically unrealistic GPP, where GPP on that day is less than −0.5. Flag 3 was applied for biologically unrealistic ER, where ER on that day was greater than 0.5. Flag 4 was applied for days where modeled K₆₀₀ was greater 20, which is unrealistically high. After flagging daily output, we prepared censored output files that eliminated days with any flag.

In total, thirteen sites were modeled to produce metabolism estimates for 5,248 site-days of which 93 percent (4,876 site-days) passed quality assurance filtering. Ten of the 13 sites were used to compare annual and seasonal low-flow trends in sloughs and ridges after removing days affected by experimental high flow releases. That data set comprised 2,597 site-days from sloughs that were compared with 1,424 site-days from ridges. Before comparing slough and ridge sites, each site with multi-year data was summarized by averaging daily data from across years using a Julian calendar. High flow metabolism responses were evaluated at two key sites, a slough site (RS1D-S) and a ridge site (RS1-R). Sensor data are available in Harvey and others (2022a) and model inputs and outputs are detailed in Harvey and others (2022b).

Results and Discussion

The primary impacts of flow, based on the 1-station streamMetabolizer model, included changes in temperature and DO (fig. 18) and large increases in ER and gas exchange at sites nearest S-152 (fig. 19, table 24). During high-flow, ER at RS1 (−31.81 g O₂/m²/d) and Z5-3 (−34.75 g O₂/m²/d) sites were double those during low-flow, and 5- to 6-fold higher than at C2 (control). K₆₀₀ at RS1 and Z5-3B also increased 6.4- and 3.8-fold, respectively, over low-flow values. Using flow-period averages as replicates, t-tests were used to evaluate flow differences for C2 and RS1. These analyses indicated high- and low-flow values were not different at C2 but were significantly different at RS1 for ER (p<10⁻⁵), K₆₀₀ (p<10⁻⁵), and NEP (p<0.01). While high flow increased GPP at the intermediate-flow sites, Z5-3B (+68 percent) and Z4-1D (+58 percent), such results are tentative given the shorter periods of record (1 flow) for those sites.