Numerical simulation of the residence time distribution dynamics of a restored agricultural wetland under varying environmental conditions

A_b

vertical turbulent mass diffusivity (m² s⁻¹)

A_h

horizontal turbulent (eddy) viscosity (m² s⁻¹)

A_ho

background horizontal turbulent viscosity (m² s⁻¹)

A_v

vertical turbulent (eddy) viscosity (m² s⁻¹)

A_vo

background vertical turbulent viscosity (m² s⁻¹)

b

buoyancy flux (-)

average pool (basin) width (m)

C

constituent concentration

C_d

discharge coefficient (m^3/2 s⁻¹)

C_(o,m)

observed (o) or modeled (m) tracer concentration (μg L⁻¹)

f

Coriolis parameter (s⁻¹)

f_e

Coriolis acceleration (m s⁻²)

f_r

fraction of shortwave radiation absorbed at the water surface (-)

g

local acceleration due to gravity (m s⁻²)

g_(o,m)

observed (o) or modeled (m) volume-based RTD function (-)

h

local water depth below a reference (m)

H(t)

modeled time-varying total water depth (m)

average pool depth as determined from interpolated, kriged, depths (m)

H_k

interpolated, kriged, water depth (m)

H_k(max)

maximum interpolated, kriged, water depth (m)

average initial depth of water above the crest of the outflow structure (m)

K_x(o,m)

observed (o) or modeled (m) longitudinal dispersion coefficient (m² s⁻¹)

L

wetland pool (basin) length (m)

MAE

mean absolute error

MARE

mean absolute relative error

MIA

modified index of agreement

M_r(o,m)

observed (o) or modeled (m) tracer mass recovery (kg)

observed (o) or modeled (m) first moment of the volume-based RTD

observed (o) or modeled (m) second central moment of the volume-based RTD

m_x,y

dimensionless coordinate transformation factors

p

water pressure (mbar)

Pe_x(o,m)

observed (o) or modeled (m) Péclet number (-)

observed (o) or modeled (m) average of inflow and outflow discharges (m³ s⁻¹)

Q_o(o,m)

observed (o) or modeled (m) outflow discharge rate (m³ min⁻¹)

sources and sinks of heat in the heat transport equation

Q_u

source and/or sink of horizontal turbulence

Q_v

source and/or sink of vertical turbulence

R_c

constituent source or sink rate (g m⁻² s⁻¹)

RRMSE

relative root mean square error

Δt

time step in real-time (min)

t_{(o,m)(i,p,50)}

initial (i) and peak (p) arrival, and median (50) detention times of the observed (o) and modeled (m) RTD

t_f

time from injection when the tracer study was ended (min)

T_w

water temperature (°C)

u

x-component of water velocity (m s⁻¹)

V_(o,m)

observed (o) or modeled (m) wetland basin volume (m³)

v

y-component of water velocity (m s⁻¹)

local wind speed (m s⁻¹)

w

vertical water velocity (m s⁻¹) in dimensionless sigma coordinates

w*

physical vertical velocity in Cartesian coordinates (m s⁻¹)

x

longitudinal direction in the momentum and transport equations (m)

y

lateral direction in the momentum and transport equations (m)

Y_(o,m)

observed (o) or modeled (m) state variables of interest

z

stretched vertical sigma coordinate in the momentum and transport equations

Greek Letters

η

free surface displacement (m): H-h

λ

vertical light extinction coefficient (m⁻¹)

_ζ(o,m)

observed (o) or modeled (m) dimensionless flow-weighted time (-)

_ζ(i,p,50)

initial (i) and peak (p) arrival, and median (50) dimensionless flow-weighted detention times of the observed (o) or modeled (m) RTD

ρ_w

density of water (g m⁻³)

dimensionless variance of observed (o) or modeled (m) RTD

τ

dummy variable of integration

Field-scale tracer studies that yield outlet residence time distributions (RTDs) are a primary technique for measuring dispersion in flow-through basins such as constructed surface water treatment wetlands. However, these tests are conducted after a system has been constructed and is in operation and can only give a post-development snapshot of mixing behavior under the environmental and flow conditions present during the study. An alternative approach is to employ multi-dimensional numerical hydrodynamic and mass transport models to simulate the mixing behavior of these systems under realistic time-varying environmental conditions (for example, Badrot-Nico et al. 2009, and Min & Wise 2009, among others). Such models can help reveal the complex physical processes that work to yield an observed RTD for a system for a set of environmental conditions, as well as be used to systematically evaluate how environmental factors affect RTD temporal features and inform how wetland basins might be modified or designed to promote beneficial mixing and mass transport patterns (Persson et al. 1999; Persson 2000; Jenkins & Greenway 2005).

Using multi-dimensional hydrodynamic and mass transport models to study mixing processes in surface water flow-through basins is particularly promising as an aid for the informed design of nutrient-reduction wetlands in agricultural watersheds. These systems, such as those developed under the Iowa Conservation Enhancement Program (Iowa CREP; https://iowaagriculture.gov/crep), are intended to operate in a natural state and, as such, are subject to ambient atmospheric conditions and receive unregulated time-varying hydrological and mass loads (Crumpton et al. 2020). Their unique design considerations support the need for tools that can be used to help ensure that wetland systems developed in agricultural landscapes are designed to promote beneficial mixing patterns. Consideration of design factors that mitigate the potential negative effects of ambient environmental conditions on mixing, such as from external wind forcing, unregulated flow-through rates, and/or internal thermal stratification, can help to ensure that these systems are built with prior knowledge of the potential range of mixing conditions for a given site.

However, despite their growing adoption as part of the wetland design process (Persson 2000; Min & Wise 2009), hydrodynamic and mass transport models have not been applied as a general design tool for wetland development in agricultural landscapes. Before these models can be adapted as a design tool, their ability to replicate observed mixing, hydrodynamics, and other system processes relevant to mass transport needs to be assessed. Such an assessment would necessarily involve calibration to and comparison with results from multiple observed field tracer studies conducted under a range of flow and atmospheric conditions and evaluation of the ability of the employed model to replicate observed RTD characteristics.

In this study, we describe the development and calibration of the three-dimensional Environmental Fluid Dynamics Code (EFDC; Hamrick 1992) numerical hydrodynamic and mass and heat transport model to six tracer studies conducted under time-varying hydrological and ambient atmospheric conditions on a single restored agricultural wetland located in central Iowa. In addition to assessing the suitability of the EFDC model for simulating wetland hydraulic and mixing processes under natural conditions and for replicating observed RTD features for these tracer studies, the calibrated models were used to systematically test the relative influences of the primary environmental effects of wind forcing, ambient hydraulic conditions, and water column thermal dynamics on the RTD characteristics of this system.

Like most Iowa CREP wetlands, this site operates in an unmanaged state and receives time-varying influent flow and nutrient mass loads and is directly exposed to local ambient atmospheric and environmental conditions. Submersed aquatic vegetation is highly prevalent in the wetland from late spring through late fall and is typically absent from the system before and after this period (Green 2016; Eeling 2019).

The bathymetry of the wetland was surveyed in early April of 2010 and 2011. Repeated surveys were performed to document the change in basin bathymetry after a large region-wide flood event in August 2010. Each survey involved taking between three and five replicate point measurements of local water depths at locations (n > 400) throughout the wetland using a surveyor's staff (±3.05 cm). The coordinates of each point were recorded using a sub-meter handheld Global Positioning System unit (GeoTrimble Corp.; Sunnyvale, CA). Measured depths were adjusted downward by the observed water level above the crest of the outflow structure on the day of each survey and merged with local Light Detection and Ranging (LiDAR) point data, vertically offset by the estimated pool boundary elevation (Iowa LiDAR Mapping Project, Accessed 2010, http://www.geotree.uni.edu/lidar/). A 1 m² resolution grid representing basin bottom elevations relative to the elevation of the outlet weir crest was derived by kriging of the merged point datasets (Green 2016).

Volume-depth and volume-area curves were developed from the derived bathymetry grid by incrementally integrating over the basin volumes and planar areas below and up to one-half meter above the weir crest reference elevation. Each depth-volume and depth-area curve was fitted to an appropriate nth-order polynomial to permit estimation of time-varying system water volumes and surface areas from measured water depths above the weir (Green 2016). These values were subsequently used for estimating hydraulic residence times and volumetric inflow rates. Spatial interpolation was performed using the gstat package in the R statistical computing environment (Pebesma 2004).

Local meteorological data, including wind speed (⁠⁠, m s⁻¹), and direction, hourly rainfall (m h⁻¹), air temperature (°C), percent relative humidity, and incident solar radiation (W m⁻²), sampled at 1-h intervals, was acquired from the Gilbert, Iowa station (site number A130219) of the Iowa State University Agricultural Climate Monitoring Network (ISU-AG; https://mesonet.agron.iastate.edu/agclimate/) located approximately four linear miles from the wetland. Air pressure data (mmHg) was obtained from the St. Cecilia, Ames, Iowa SchoolNet Weather Monitoring Network station (site number SAMI4; https://mesonet.agron.iastate.edu/schoolnet/). Measured wind speeds were extrapolated to 10 m above the wetland basin using a power-law vertical wind speed profile suitable for lightly vegetated regions (Irwin 1979).

Time-varying wetland water elevations were monitored at 5-min intervals at a single stilling well near the basin outlet (as noted in Figure 1) using a pressure transducer (Solinst, Georgetown, ON, Canada). Dynamic water surface depths above the crest of the outlet weir were estimated by offsetting the measured stilling well depths by the average water depth at the stilling well when the wetland water volume was at its full pool capacity. Continuously measured depths were used to estimate dynamic wetland pool volumes using the regression equations developed from the hypsographic curves.

Mean water velocity and depths were measured at 5 min intervals in the stream channel below the wetland outflow structure using a submerged area velocity meter and stage recorder (Solinst, Georgetown, ON, Canada). Manual discharge measurements at this location were made using the mid-section or 0.6 depth methods (Buchanan & Somers 1969) and a handheld side-looking two-dimensional Sontek Flow-tracker Doppler velocimeter (Sontek, San Diego, CA). Manual flow measurements were used to develop rating curves for the channel directly downstream of the wetland outlet control structure and to derive the coefficients for the discharge equation for the weir. Because of a lack of instrumentation in the inlet channel, time-varying inflow rates were estimated using measured volumetric outflow rates, estimates of dynamic pool volumes and water surface areas, and the reverse level-pool routing procedure described by Zoppou (1999).

Each tracer study (Table 2) entailed instantaneously injecting a known mass of rhodamine WT dye (RWT, 20% solution, Organic Dye Stuff Corporation, Providence, RI) across the wetland inlet channel and continuously measuring concentrations at the outlet. To minimize initial density effects, RWT was diluted with approximately five gallons of stream water before injection. Dye concentrations were measured directly upstream of the outflow structure every 5 min starting at least one-half hour before tracer injection using a Turner Designs Cyclops 7 fluorometer (Turner Designs, Sunnyvale, CA). The fluorometer was secured in a custom-made perforated and screened gray cylindrical plastic housing designed to minimize sunlight exposure and to prevent measurement interference by floating detritus and algae (Green 2016).

Table 2

Tracer study durations, mass recoveries, and mean flow and wind conditions

Study .	Date/time start .	Date/time end .	Dur. (days) .	Mass Inj.(kg) .	Mass rec. (kg) .	Mass rec. (%) .	(m3 min−1) .	(m s−1) .
TS1	2009-11-22 17:17	2009-12-01 06:42	16.6	0.113	0.092	81.4	1.3	3.5
TS2	2010-05-06 18:04	2010-05-13 22:50	7.2	0.113	0.11	97.3	3.1	4.0
TS3	2010-05-18 20:18	2010-05-25 15:28	6.8	0.113	0.076	67.3	3.1	2.9
TS4	2011-04-18 12:01	2011-04-21 03:10	2.6	0.0226	0.013	57.5	4.2	4.0
TS5	2011-05-22 18:09	2011-05-25 21:04	3.1	0.0226	0.026	115.0	6.3	3.5
TS6	2012-04-23 18:55	2012-05-03 12:34	9.7	0.0226	0.017	75.2	1.8	3.4

Study .	Date/time start .	Date/time end .	Dur. (days) .	Mass Inj.(kg) .	Mass rec. (kg) .	Mass rec. (%) .	(m3 min−1) .	(m s−1) .
TS1	2009-11-22 17:17	2009-12-01 06:42	16.6	0.113	0.092	81.4	1.3	3.5
TS2	2010-05-06 18:04	2010-05-13 22:50	7.2	0.113	0.11	97.3	3.1	4.0
TS3	2010-05-18 20:18	2010-05-25 15:28	6.8	0.113	0.076	67.3	3.1	2.9
TS4	2011-04-18 12:01	2011-04-21 03:10	2.6	0.0226	0.013	57.5	4.2	4.0
TS5	2011-05-22 18:09	2011-05-25 21:04	3.1	0.0226	0.026	115.0	6.3	3.5
TS6	2012-04-23 18:55	2012-05-03 12:34	9.7	0.0226	0.017	75.2	1.8	3.4

^a equals the average of the instantaneous observed inflow and outflow rates averaged over the course of the study.

^bEquals the average wind speed observed during the study.

View Large

Tracer response curve data conditioning entailed correcting for background fluorescence, removing erroneous concentration measurements resulting from instrument error, and accounting for temperature effects (Smart & Laidlaw 1977). The studies TS1, TS4, and TS5 were terminated before all of the dye had exited the system, and the tails of these response curves were extrapolated to measured background concentrations using an exponential decay profile fitted to the descending limb (Green 2016).

EFDC is a three-dimensional flow and mass transport model designed to simulate the hydrodynamics, mixing, and water quality processes of shallow surface water systems (Hamrick 1992; Hamrick & Wu 1997). This model has been used extensively in lake (c.f. Jin et al. 2000, 2002, 2007; Wu & Xu 2011), reservoir (c.f. Li et al. 2010; Zhang et al. 2013), riverine (Huang et al. 2008; Franceschini & Tsai 2010), and coastal and estuarine (Wool et al. 2003) hydrodynamic and biogeochemical modeling studies. However, EFDC, to our knowledge, has not been used to analyze the mixing and temperature dynamics of extremely shallow flow-through surface water constructed wetlands. EFDC can simulate eutrophication processes through a linkage to the HEM3D water quality model (Hamrick & Wu 1997) and temperature dynamics through a linkage to the CE-QUAL-W2 equilibrium heat exchange model (Cole & Wells 2003).

In the above equations, H is the total water depth (m); η is the free surface displacement depth (m); p is water pressure (mbar); u, v, and w are the horizontal, lateral, and vertical velocity components (m s⁻¹); x and y are the horizontal Cartesian coordinates (m) and z the stretched vertical sigma coordinate; h is the depth below a relative reference coordinate z*; A_v is the vertical turbulent (eddy) viscosity (m² s⁻¹); g is acceleration due to gravity (9.81 m s⁻²), and b is the buoyancy flux (-), defined as the difference between water density, ρ_w (g m⁻³), and a reference density, ρ_w(0). The coefficients m_x and m_y are used to transform the governing equations into their Cartesian equivalents (Ji 2008). The Q_u and Q_v terms in Equations (1) and (2) represent sources and sinks of horizontal and vertical turbulence, respectively, including the effects of wind shear. Coriolis acceleration, f_e, is incorporated into the x and y momentum equations and is determined by the Coriolis parameter, f (9.761 × 10⁻⁵ s⁻¹ for this site), and local accelerations induced by grid curvature.

The basin depth grids derived from each bathymetric survey defined the bathymetry of the model domain for the pre- and post-flood tracer study simulations. The computational domain was developed to feature an orthogonal curvilinear grid in the inlet channel merged with a quasi-Cartesian grid for the central pool (Green 2016). The model domain, used with both bathymetric survey results, consisted of 3,211 horizontal cells and 5 vertical layers, totaling 16,055 active grid cells. Average horizontal cell sizes for the pool and inlet channel sections were 2 and 1.1 m², respectively. This model grid configuration required a time step of 0.3 s to maintain numerical stability. The suitability of five layers to resolve vertical velocity, temperature, and mixing dynamics in this system is consistent with other numerical models of shallow surface water systems (e.g. Jin et al. 2007). Hydrodynamic roughness lengths, used to calculate internal friction losses of momentum in EFDC, were estimated from surveyed bathymetry and had little influence on simulation results (Green 2016).

Time series of estimated average hourly inflow rates (m³ s⁻¹) were supplied as an upstream volumetric flux boundary condition for each simulation. For the simulation TS1, continuous outflow monitoring was terminated on 2009-12-01. Outflow discharges after the termination of flow monitoring for this study were estimated from scaled discharge measured at the United States Geological Survey (USGS) monitoring station Ioway Creek at Ames (site number 05470500; https://waterdata.usgs.gov/nwis/inventory/?site_no=05470500), located approximately 20 miles from the wetland.

H is the modeled water depth in model cells directly upstream of the weir (m), and C_d is the discharge coefficient specific to the 4 m wide weir for this wetland (4.23 m^3/2 s⁻¹). The discharge coefficient for the wetland outflow structure was developed based on the weir dimensions and outflow discharge measurements obtained in the channel below the outlet (Crumpton et al. 2020).

Each tracer study was simulated by applying an instantaneous injection of dye of a known concentration at the upstream boundary of the model domain. This concentration was estimated from the mass injected for each respective field study (Table 2) and the instantaneous volumetric flow rate at the time of injection. Time series of inlet temperatures for upstream thermal forcing were estimated from a linear regression of air temperatures observed at the ISU-AG weather station against hourly averages of instantaneous observed inlet channel temperatures for the period late March through late May 2012 for this wetland (R² = 0.67; Green 2016). Linear regression between air and surface water temperatures has been used in other studies to reasonably approximate missing water temperature data (Pilgrim et al. 1998; Devkota et al. 2013). For each simulation, inflow temperatures were assumed uniform over the inlet channel width and depth.

For each calibration for each simulation, the intrinsic background vertical and horizontal turbulent mass transport parameters, A_vo and A_ho (Ji 2008) of the hydrodynamic and mass transport governing equations, and the fraction of shortwave solar radiation reaching the water surface (f_r; -) and the vertical light extinction coefficient (λ; m⁻¹), influencing the net water surface-atmosphere and water column-sediment heat fluxes, were systematically adjusted to obtain best fits between modeled and observed values for dye concentrations (see Green 2016 for more detail). These parameters were adjusted over the ranges 1 × 10⁻⁶–0.01 m² s⁻¹; 1 × 10⁻⁴–0.025 m² s⁻¹; 0.45–0.90 (-); and 0.1–0.9 m⁻¹, respectively. During the calibration process, greater emphasis was placed on obtaining reasonable fits for modeled dye concentrations and temperatures than for hydraulic variables, as the primary focus of this study was the replication of observed tracer response curves and internal temperature time series.

The dimensionless initial and peak arrival times are considered to be reasonable metrics of tracer short-circuiting in this wetland. Thackston et al. (1987) suggest that ‘significant’ short-circuiting is indicated by an initial dimensionless arrival time, ξ_i, of less than 0.2. Likewise, Persson et al. (1999) suggest that significant short-circuiting and diminished basin hydraulic efficiency are implied for dimensionless peak arrival times, ξ_p, of less than 0.75. We adopt these general conventions as well.

As shown in Table 4, basin water depths (referenced to the maximum basin depth) of best-fit simulations reasonably matched observed values for most simulations. Despite the comparatively large percent error and lower MIA shown for some of the simulations, MAE ranged from only 1 mm to nearly 1 cm. The estimates of MAE for modeled water depths are significantly less than the errors reported in other studies (Jin et al. 2000; Li et al. 2010). Basin volumes for each simulation period exhibited error levels similar to those reported for modeled depths (not shown). While the RRMSE and MARE estimates, ranging from 12.6 to 37%, indicate fairly significant overall error in modeled volumes for some model runs, the moderately small MAE (2.64–79.3 m³) and maximum difference estimates (11.8–266.8 m³) suggest reasonably accurate dynamic volume estimates over all simulations. Both modeled depths and volumes were highly sensitive to the value of A_vo, but were not affected by A_ho (data not shown). In general, higher values of A_vo resulted in greater depths and larger modeled volumes. The maximum–minimum prediction bounds for these modeled time-varying variables showed that unique A_ho and A_vo parameter combinations result in widely varying model responses.

Modeled outflow discharges corresponded closely with observed values and, in all cases, resulted in MARE and RRMSE of less than 5% and MIA greater than 0.8. The high degree of agreement between modeled and observed outflow discharges across all studies suggests that EFDC can realistically simulate outlet rates for this system under a range of flow conditions.

Time-varying temperatures resulting from each best-fit simulation reasonably matched temperatures observed at the outlet monitoring location and the stilling well, as indicated by relatively high MIA and relatively low MAE, MARE, and RRMSE for both monitoring locations (Table 4). However, the accuracy of the modeled temperature time series is distinctively variable between simulation periods, as is shown in Table 4 and Figures 2–4.

For all simulations, the model better reproduced temperatures in the outlet channel than in the stilling well. Despite their proximity, the difference in model accuracy between these two locations may be attributable to the isolating effects of the stilling well casing or to the uncertainty of the exact depth at which the pressure transducer thermistor was located within the stilling well. All fit statistics in Table 4 exhibited a similar discrepancy between the outlet channel and stilling well monitoring locations. A notable exception is the high degree of difference between outflow and stilling well temperatures for the study TS6. In this simulation, modeled temperatures within the stilling well were, on average, 2.32 °C less than observed temperatures, with this difference often exceeding 7 °C.

Modeled temperatures were highly sensitive to the magnitude of A_vo, and relatively insensitive to A_ho, f_r, and λ. Higher values of A_vo, corresponding to increased ambient background vertical mixing rates, tended to result in lower overall predicted temperatures and dampened diurnal temperature fluctuations at both monitoring locations. In contrast, A_ho was observed to exert little influence on modeled temperature patterns for all of the simulations, suggesting that horizontal mixing processes have little influence on the internal temperature dynamics of this wetland. The temperature calibration coefficients f_r and λ influenced the magnitude of modeled temperature curves but did not influence diurnal patterns at both monitoring locations. Of these two parameters, f_r maintained a stronger influence on model results. In general, decreases in f_r, representing diminished short-wave radiation absorption at the water surface, resulted in higher overall modeled temperatures; however, this effect was comparatively small, resulting in no more than a 0.5 °C maximum difference in modeled and observed temperatures at both locations over the tested range. All error statistics exhibited similarly minor changes over the tested range of this parameter.

In general, lower values of λ resulted in overall elevated temperatures at both monitoring locations, suggesting increased sediment bed heating from reduced short-wave radiation attenuation. The lack of data on the actual values of these parameters for this system over the simulation periods tested necessarily precludes inference about their validity and respective roles in influencing modeled temperatures beyond the qualitative descriptions given here. The best-fit temperature calibration coefficients are within the range of those reported by several authors (Jin et al. 2000, 2002; Li et al. 2010; Devkota et al. 2013; Jin & Ji 2013).

Due to an absence of instrumentation, inlet water temperatures were estimated from the linear regression model, as discussed previously, between air and water temperatures. The comparatively high correspondence between modeled and observed temperatures at both monitoring locations for all studies but TS6 indicates that the model is relatively insensitive to inaccuracies in inlet boundary temperatures. This finding suggests that the temperature dynamics of this system are driven primarily by in-pool heat exchanges between the water column and the atmosphere and, less importantly, between the water column and the wetland sediment, as also concluded by Sweeney et al. (2005), Badrot-Nico et al. (2009), and Eeling (2019). However, using this regression model to estimate inlet temperatures cannot be discounted as a potentially significant source of the error reported for the temperature simulation results. Regardless, considering the comparatively good fit between observed and modeled temperatures for every simulation but one, the total overall error from this approximation is likely small, but the exact error cannot be quantified at this time.

Tracer concentrations of the selected best-fit simulations, as observed at the monitoring point in the outlet channel, were in reasonable agreement with observed tracer concentrations, with average MARE, RRMSE, and MAE of 5.4%, 9.2%, and 0.16 μg L⁻¹, respectively (Table 4). Modeled tracer response curves for several studies – notably TS1, TS3, and TS5 – featured earlier arrival and peak times, and sharper fronts than corresponding observed curves. Additionally, in every case but TS3 and TS6, modeled peak concentrations were between 1 and 5 μg L⁻¹ less than observed (corresponding to percent differences between 7 and 20%). In the cases of these two exceptions, modeled peak concentrations were between 0.5 and 1 μg L⁻¹ greater than observed peak concentrations. Overall, modeled tracer response concentration curves more closely matched late-time concentrations in descending limbs than early-time concentrations up to and including the peak response time.

The calibration process revealed that the overall shape and initial and peak arrival times were highly sensitive to the magnitudes of both A_vo and A_ho. Optimized values of A_ho ranged from 2.5 × 10⁻² to 1 × 10⁻⁴ m² s⁻¹. In general, larger values of A_ho induced earlier initial arrival times and longer and smoother response curve tails. Optimized values of A_vo ranged from 2.5 × 10⁻⁴ to 1 × 10⁻⁶ m² s⁻¹, with smaller values resulting in earlier tracer initial and peak arrival times and sharper ascending limbs. That different combinations of A_vo and A_ho produced variable best-fit statistics for the simulations considered suggests that these parameters are likely sensitive to ambient hydraulic and atmospheric conditions and are thus moderately time-varying in most cases, at least for small surface water bodies such as Iowa CREP wetlands.

The sensitivity of the early time features of the modeled tracer response curves to the varying combinations of A_vo and A_ho tested in this work is illustrated in the minimum-maximum bounds of modeled dye concentrations for each simulation period, as shown in Figures 2–4. These bounds also show the relative insensitivity of response curve tails to the magnitudes of A_vo and A_ho. The shapes of these bounding curves suggest that model simulations intended to calibrate against initial and peak tracer arrival times – both metrics are commonly considered to be indicators of the degree of short-circuiting within a basin (Thackston et al. 1987; Persson 2000) – will be highly sensitive to values of these parameters. However, if simulations are intended to emphasize mean or median residence times over short-circuiting metrics, these parameters appear to be less important.

The TS5 tracer study was conducted when submersed aquatic vegetation was in its nascent growth phase. During this tracer study, the area-weighted average percent cover of submersed aquatic vegetation was 11% (Green 2016). The presence of vegetation in the basin during this tracer study may have influenced the discrepancy between the observed and modeled RTD characteristics for this study. The EFDC simulation for this tracer study does not consider the presence of vegetation and shows a sharply earlier peak arrival time of the tracer than what was observed, although the initial arrival time and mean detention time roughly agree.

As part of a broader study of the RTD characteristics of constructed agricultural wetlands, Green (2016) found that peak arrival times were delayed for tracer studies conducted under vegetated conditions for this system as well as others and observed that submersed vegetation at all growth phases could trap tracer, resulting in an increase in the mean tracer detention time above the mean hydraulic residence time, as well as an increase in general dispersion, as reflected in and K_x. In our opinion, the exclusion of vegetation in the model structure for TS5 is the likeliest explanation for the significant under-prediction of K_x for this study.

Green (2016) showed that environmental effects such as wind speed, wind direction, and ambient flow conditions can have significant but varying influences on considered tracer response curve characteristics for this wetland. Applied wind shear at the water surface and influent volumetric flow rates supply the primary energy sources available for mixing to this system, with the effects of wind being particularly important in the absence of vegetation. Wind shear applied to the water surface can act as a stirring mechanism, enhancing horizontal and vertical turbulent mixing rates and increasing internal rates of local advection (Fischer et al. 1979).

An additional environmental effect that may strongly influence mixing and RTD characteristics is temperature-induced vertical and horizontal thermal heterogeneity. For instance, Macdonald & Ernst (1986) and Sweeney et al. (2005) report that vertical thermal stratification may be a significant cause of short-circuiting in moderately shallow flow-through basins. What is less clear is whether vertical and horizontal variances in internal basin temperatures can affect bulk mixing and RTD characteristics for extremely shallow systems that are also highly wind-exposed. The influence of temperature-induced density differences on the hydraulic behavior of Iowa CREP wetlands may be significant, as these systems receive the majority of their flows from upland subsurface drainage tile networks, which during mid-spring to late summer tend to export water with temperatures that can be significantly lower than ambient water column temperatures (Green 2016; Eeling 2019).

To isolate the relative influences of environmental effects on the varying observed RTD characteristics of this wetland, we conducted a series of sensitivity analyses with each of the calibrated best-fit EFDC simulations. These sensitivity analyses involved performing additional simulations that systematically and alternately excluded wind, atmospheric radiative, and inlet thermal forcing from the model. Additionally, to isolate the influence of observed time-varying flow rates on modeled RTD characteristics, we conducted alternate steady flow simulations for which the inflow boundary discharge was set to the flow rate observed at the time of dye injection for each tracer study. Seven sensitivity scenarios were performed for each modeled tracer study (Table 5).

Modeled tracer response curves for these sensitivity runs were subsequently analyzed using the presented RTD analysis techniques. The derived RTD statistics were compared with corresponding statistics for the best-fit simulations. The presumption underlying this analysis was that the calibrated modeled temperatures and tracer response curves reasonably represent the interior tracer dispersion, velocity, and temperature dynamics of this system over the periods studied despite the comparatively small number of calibration locations for temperature fit assessment and the linear regression used to estimate inlet boundary temperatures.

The exclusion of thermal effects with the retention of wind forcing, regardless of the specification of time-varying flow boundary conditions (scenarios 6 and 7), was shown to have only moderate influence on each of the considered RTD statistics, except for ξ₅₀ and ⁠. For these scenarios, these statistics, while exhibiting only nominal average changes, demonstrated significant percent changes, ranging from −84 to 56% and −50 to 40% for ξ₅₀ and ⁠, respectively, with moderate differences between steady and unsteady flow conditions.

When temperature and wind forcing were both removed from the simulations (scenarios 3 and 5) ξ_i and ξ_p increased, on average, by nearly 200 and 240%, respectively, above the best-fit cases, corresponding to ranges of 62–402% and 8–644%. There were no significant differences in the responses of these statistics to flow conditions. Percent changes in and ξ₅₀ ranged between −68 and 22%, and −32 to 33%, respectively, corresponding to mean changes of −3 and 3%, with nominal differences in responses between steady and unsteady flow conditions. In contrast, and K_x decreased on average by approximately 55 and 46%, respectively, corresponding to responses ranging from −65 to −33% and −88 to 27%, with exhibiting a broader response range under steady flows, and K_x exhibiting a wider range under unsteady flow conditions. The Péclet numbers for these scenarios increased, on average, by nearly 100% above the best-fit case, corresponding to increases of between 37 and 163%, and with only marginal differences between steady and unsteady flow conditions. Mean increases in ξ_i and ξ_p under these scenarios are between 70 and 200% greater, and and K_x are lower on average by around 250 and 400% than those observed for the no-wind cases discussed previously (scenarios 2 and 4), regardless of flow conditions.

These results suggest that wind forcing exerts the greatest influence on most considered RTD statistics, particularly on measures of short-circuiting and, to a lesser degree, indicators of basin-scale mixing. While temperature dynamics in isolation have, when compared to the effect of wind, a moderate influence on tracer short-circuiting, as can be seen from scenarios 6 and 7, this effect, when considered in concert with wind forcing, can have a marked influence on measures of internal mixing and dispersion rates in this system. This can be seen when comparing the successive exclusion of wind and temperature dynamics in scenarios 2 and 4, and 3 and 5 (Figure 6). In these scenarios, both short-circuiting and mixing are considerably diminished in the absence of wind forcing alone, but more significantly diminished in the absence of both thermal and wind forcing. This result suggests an interaction between these effects on observed short-circuiting and mixing.

Wind forcing has long been considered a potential driver of mixing in shallow flow-through basins, such as the one featured in this work. For instance, Thackston et al. (1987) and Watters et al. (1973) noted that wind-driven surface currents accompanied by bottom-layer return currents, when present, may cause these systems to become more fully mixed by increasing both turbulent diffusion and shear dispersion. Indeed, Watters et al. (1973) demonstrated in wind-tunnel experiments the tendency for increasing wind shear to increase longitudinal dispersion rates in enclosed basins, with dispersion rates also strongly correlated to the direction of the applied wind (a factor not considered in this study). Their results confirm the theoretical work of Wu (1969), which is also supported here. Whether apparent increases in dispersion resulting from applied wind shear can be attributable to wind-induced differential advection or increases in the overall degree of turbulence is uncertain at this time.

The underlying mechanisms resulting in increased dispersion and short-circuiting from temperature variations – whether independent of or in concert with wind forcing – are poorly understood. However, other researchers have noted that these basins tend to stratify daily during warm periods, a process that is typically characterized by daytime stratification followed by nighttime convective mixing, with some systems exhibiting longer-duration quasi-stable stratification over successive daytime-nighttime periods (Gu et al. 1996; Eeling 2019). This is especially true if influent water temperatures are significantly lower than receiving water temperatures, as occurs in the wetland studied in this work generally from late spring through early fall (Eeling 2019). In this case, colder and denser influent water may plunge into the basin, causing some tracer to remain confined within slower, deeper, and cooler layers of the pool, which are generally isolated from the effects of wind-induced stirring. This deeper-lying influent might be subjected to nighttime convective vertical mixing, which transports bottom water to the faster-moving, wind-exposed upper layers of the water column.

While thermal stratification was not explicitly measured during this study – but is explicitly represented in the calibrated simulations – the recent monitoring work of Eeling (2019) suggests that this is likely the case for tracer studies conducted when temperature differences between the influent water and the pool were greatest (i.e. all studies but TS1). It is possible that in the absence of this effect (as with scenarios 3, 5, 6, and 7), water that is transported primarily along the bottom of the basin might not become vertically mixed enough during transit to be significantly affected by the faster-moving water in the upper layer of the water column, as is indicated by scenarios 6 and 7.

Indeed, Macdonald & Ernst (1986) attributed observed short-circuiting in maturation ponds to vertical stratification, a conclusion also made by Pedahzur et al. (1993) for tracer studies conducted in a shallow stabilization pond. Sweeney et al. (2005) also found, through multi-dimensional heat and mass transport modeling of a waste stabilization lagoon, that diurnal temperature inhomogeneity, particularly vertical stratification, in these systems may have a large influence on their hydraulic and mixing characteristics, particularly short-circuiting. These researchers observed that short-circuiting might be most pronounced in the areas of these types of flow-through basins that are more frequently stratified, which in the wetland studied in this work would correspond to the remnant deeper channel that traverses the basin length from inlet to outlet, as can be shown in Figure 1.

Neither temperature dynamics nor wind effects appear to have had a strong influence on the median and mean residence times, suggesting that mean basin flow rates mostly influence these RTD characteristics. This finding contrasts with the observations of Bentzen et al. (2008), who found through numerical simulation of wind effects on the RTD characteristics of a shallow highway detention pond that both wind speed and direction strongly influenced mean tracer detention times in their study system.

This study demonstrates the efficacy of the EFDC model in simulating the hydraulic, temperature, and tracer transport dynamics of a constructed agricultural wetland over a range of flow and environmental conditions. Of the state variables modeled, time-varying water depths, wetland volumes, and outflow discharges featured the lowest overall errors, with mean absolute relative errors ranging from approximately 0.02–0.43% for water elevations, 0.03–2.01% for basin volumes, and 0.05–3.1% for outflow rates. For dynamic temperature and tracer concentrations, the differences between modeled and observed values were reasonably low, ranging from 0.3 to approximately 6.2% for tracer concentrations and 0.6–16% for temperatures, supporting the suitability of EFDC as a general modeling tool for evaluating the temperature and mixing dynamics of these types of very shallow surface water systems.

In addition to accurately simulating the raw, untransformed, dynamic state variables of interest in this system, the EFDC model was shown to be reasonably capable of replicating most of the considered calculated volume-based RTD curve statistics for each tracer study, save for one conducted during a period in which submersed aquatic vegetation was beginning to become established but still of low density, but which occupied a significant area of the wetland. Whether the presence of vegetation explains the relatively significant differences between modeled and observed RTD statistics for this study is unknown, but prior research into the effects of vegetation on tracer transport suggests that aquatic vegetation can strongly influence RTD characteristics in surface water wetlands.

While the model was seemingly capable of replicating some dynamic variables for the tracer studies, the variable vertical and horizontal turbulent mixing coefficients (A_vo and A_ho) were found to produce the best correspondence between observed and modeled tracer response curves, which varied significantly between tracer studies, differing in some cases by several orders of magnitude. Further testing of these variables on the modeled response curves provided a large range of potential predicted responses for tracer concentrations. These findings suggest that any application of this model to test mixing and temperature dynamics in shallow flow-through basins, absent a set of corresponding field-scale tracer studies and continuous in-pool temperature measurements, must carefully consider the values of vertical and horizontal turbulent mixing rate coefficients. Any conceptual test application of the model to a non-tested system should account for a range of possible tracer response curve responses as dictated by a span of A_vo and A_ho values of several orders of magnitude. Before adopting EFDC as a wetland design tool, future work should also focus on determining the relative influence of ambient wind, temperature, and flow conditions on dictating the approximate magnitude of these primary model calibration coefficients.

The sensitivity analyses suggest that ambient flow conditions, on average, maintain minimal impact on most of the RTD characteristics considered, save for the mean and median residence times. However, for some tracer tests, early and late-time changes in flow rates were shown to have an impact on median and mean detention times and estimates of short-circuiting despite the purported ability of the volume-based RTD function to remove these effects through the implementation of the flow-weighted time scheme. Neither early nor late-time changes in flow rates were shown to have an appreciable influence on observed RTD variances and derived measures of bulk basin mixing.

Environmental effects testing using the calibrated EFDC models indicated that wind and temperature dynamics strongly influence RTD characteristics and derived basin mixing indices on average. The relative contribution of these effects on the ultimate shape of the RTD observed at a system outlet, however, appears to be variable and is likely dictated by factors such as the total degree of vertical thermal stratification and horizontal temperature variation within the basin, time-varying temperature differences between influent and ambient pool water, and the timing, magnitude, and direction of wind forcing. Both wind and temperature were shown to have only a nominal influence on calculated median and mean detention times. The influence of horizontal and vertical temperature inhomogeneity on RTD development in flow-through surface water wetlands should be explored more in future research.

This work was partly funded by the Iowa Department of Agriculture and Land Stewardship and by the Department of Ecology, Evolution, and Organismal Biology at Iowa State University. Gratitude is given to Dr Greg Stenback for his insightful review of this manuscript and for providing flow and water depth data. Additional thanks are given to Claudette Sandoval-Green, Ian Ellickson, and several undergraduate students at Iowa State University who assisted with collecting field data. Additionally, we are grateful to the wetland owners for permitting us to study this site.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.