This work demonstrates the development, calibration, and use of the three-dimensional Environmental Fluid Dynamics Code (EFDC) hydrodynamic and mass transport model to simulate mixing patterns and study the environmental controls on the residence time distributions of a 1.3 ha agricultural wetland in central Iowa. Incorporating time-varying flow boundary conditions and atmospheric forcing, the model was calibrated against observed state variables, including water temperatures, basin hydraulic characteristics, and dye concentrations monitored at the outlet for six tracer studies conducted under varying flow and atmospheric conditions when submersed aquatic vegetation was mostly absent from the basin. EFDC reasonably reproduced observed basin internal hydraulics, temperatures, and mass transport dynamics, with mean absolute relative errors ranging from 0.02 to 16.3%. Sensitivity analyses suggest that wind shear exerts the greatest control on the modeled, and by extension observed, RTD for this system, primarily affecting measures of short-circuiting and, to a lesser degree, basin-wide mixing, particularly in the absence of atmospheric thermal forcing. Thermal forcing was found to significantly influence short-circuiting and mixing during warmer periods, with this effect being highly influenced by wind. Transient flows nominally influenced most RTD characteristics, save for mean and median residence times.

  • Three-dimensional numerical hydrodynamic, mass, and heat transport models can be used in wetland design.

  • A fully parameterized and calibrated model of an agricultural wetland shows these models can represent wetland mixing under a range of conditions.

  • Wind forcing has a significant effect on mixing in this wetland, and by extension can for similar systems.

  • Thermal forcing strongly affects mixing in this wetland.

Ab

vertical turbulent mass diffusivity (m2 s−1)

Ah

horizontal turbulent (eddy) viscosity (m2 s−1)

Aho

background horizontal turbulent viscosity (m2 s−1)

Av

vertical turbulent (eddy) viscosity (m2 s−1)

Avo

background vertical turbulent viscosity (m2 s−1)

b

buoyancy flux (-)

average pool (basin) width (m)

C

constituent concentration

Cd

discharge coefficient (m3/2 s−1)

C(o,m)

observed (o) or modeled (m) tracer concentration (μg L−1)

f

Coriolis parameter (s−1)

fe

Coriolis acceleration (m s−2)

fr

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

g

local acceleration due to gravity (m s−2)

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)

Hk

interpolated, kriged, water depth (m)

Hk(max)

maximum interpolated, kriged, water depth (m)

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

Kx(o,m)

observed (o) or modeled (m) longitudinal dispersion coefficient (m2 s−1)

L

wetland pool (basin) length (m)

MAE

mean absolute error

MARE

mean absolute relative error

MIA

modified index of agreement

Mr(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

mx,y

dimensionless coordinate transformation factors

p

water pressure (mbar)

Pex(o,m)

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

observed (o) or modeled (m) average of inflow and outflow discharges (m3 s−1)

Qo(o,m)

observed (o) or modeled (m) outflow discharge rate (m3 min−1)

sources and sinks of heat in the heat transport equation

Qu

source and/or sink of horizontal turbulence

Qv

source and/or sink of vertical turbulence

Rc

constituent source or sink rate (g m−2 s−1)

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

tf

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

Tw

water temperature (°C)

u

x-component of water velocity (m s−1)

V(o,m)

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

v

y-component of water velocity (m s−1)

local wind speed (m s−1)

w

vertical water velocity (m s−1) in dimensionless sigma coordinates

w*

physical vertical velocity in Cartesian coordinates (m s−1)

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−1)

ζ(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−3)

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.

Study area

The wetland studied in this work is a shallow 1.34-ha surface water impoundment located in northern Story County, Iowa (Figure 1), developed in 2005 as part of the Iowa CREP to intercept agricultural runoff from a 300-ha farmed, tile-drained watershed. The inlet to the wetland basin is a shallow and narrow stream channel that gradually merges with a deeper and wider, wind-exposed, and seasonally vegetated pool of variable depth (Table 1). Water outlets from the system over a 4 m wide concrete broad-crested weir built into an earthen dam.
Table 1

Estimated bulk morphometric characteristics from each bathymetric survey

Survey date(s)Survey points (m) (m) (m)
April, 2010 449 0.59 1.96 0.012 
April, 2011 789 0.52 1.81 0.010 
Survey date(s)Survey points (m) (m) (m)
April, 2010 449 0.59 1.96 0.012 
April, 2011 789 0.52 1.81 0.010 

, H(k)max, and are the mean and maximum wetland water depth relative to the crest of the outlet control structure, and the initial water depth respectively.

Figure 1

Site map and bathymetry derived from the April 2010 bathymetric survey.

Figure 1

Site map and bathymetry derived from the April 2010 bathymetric survey.

Close modal

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).

Environmental monitoring and field surveys

Bathymetric surveys

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 m2 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).

Meteorological and hydrological monitoring

Local meteorological data, including wind speed (, m s−1), and direction, hourly rainfall (m h−1), air temperature (°C), percent relative humidity, and incident solar radiation (W m−2), 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).

Field-scale tracer studies

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

StudyDate/time startDate/time endDur. (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 
StudyDate/time startDate/time endDur. (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.

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).

Numerical modeling

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).

EFDC solves the three-dimensional turbulence-averaged equations of momentum (Equations (1)–(7)) using a Cartesian or orthogonal curvilinear grid system in the horizontal, with a non-dimensional sigma-stretch coordinate system in the vertical. The governing equations are solved numerically using second-order accurate spatial finite differences and a second-order accurate three-time level finite difference scheme for time integration (Hamrick 1992; Hamrick & Wu 1997). The three-dimensional momentum and continuity equations in curvilinear coordinates are defined as follows (Hamrick & Wu 1997; Ji 2008):
(1)
(2)
(3)
(4)
(5)
(6)
(7)
where Equations (1) and (2) are the momentum equations in the x and y directions, Equation (3) is the hydrostatic pressure approximation, and Equations (4) and (5) are the continuity and depth-integrated continuity equations, respectively. Equation (7) defines the vertical velocity (w; m s−1) in sigma-stretched coordinates as a function of the velocity in non-transformed vertical coordinates (w*). Equation (6) is used to convert between curvilinear and Cartesian coordinates.

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−1); 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*; Av is the vertical turbulent (eddy) viscosity (m2 s−1); g is acceleration due to gravity (9.81 m s−2), and b is the buoyancy flux (-), defined as the difference between water density, ρw (g m−3), and a reference density, ρw(0). The coefficients mx and my are used to transform the governing equations into their Cartesian equivalents (Ji 2008). The Qu and Qv terms in Equations (1) and (2) represent sources and sinks of horizontal and vertical turbulence, respectively, including the effects of wind shear. Coriolis acceleration, fe, is incorporated into the x and y momentum equations and is determined by the Coriolis parameter, f (9.761 × 10−5 s−1 for this site), and local accelerations induced by grid curvature.

EFDC couples the continuity and momentum equations to a generic three-dimensional transport equation for dissolved conservative and reactive scalar substances (Hamrick & Wu 1997; Ji 2008):
(8)
where C is the spatially and temporally varying constituent concentration; Ah and Ab are the horizontal (x and y) eddy viscosity (m2 s−1), and vertical turbulent mass diffusivity (m2 s−1) respectively; and Rc is a lumped first-order constituent reactive source and sink (loss) rate parameter (s−1). The constituent loss rate parameter Rc for each simulation was estimated as the quotient of the fractional tracer recovery and the test duration (Table 2).
The equation governing three-dimensional heat transport is defined as follows (Hamrick & Wu 1997; Ji 2008):
(9)
where Tw is water temperature, and represents inlet, outlet, atmospheric, and sediment sources and sinks of heat to the model domain (Cole & Wells 2003; Ji 2008).

Model development

Model domain

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 m2, 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).

Boundary conditions

Time series of estimated average hourly inflow rates (m3 s−1) 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.

For all simulations, the model outflow boundary condition was designated as a control structure, with outflow rates (m3 s−1) calculated from the discharge equation for a horizontal broad crested weir:
(10)

H is the modeled water depth in model cells directly upstream of the weir (m), and Cd is the discharge coefficient specific to the 4 m wide weir for this wetland (4.23 m3/2 s−1). 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 (R2 = 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.

Model calibrations

Each simulation was calibrated against tracer concentrations observed at the wetland outlet, water temperatures observed at the basin stilling well and outlet channel, dynamic basin volumes and water levels estimated from measured water surface elevations (WSE) at the basin stilling well, and outflow discharge. The accuracy of each simulation was assessed using several performance statistics, including the Modified Index of Agreement (MIA; Willmott et al. 1985), the relative root mean square error (RRMSE; Ji 2008), the mean absolute error (MAE; Legates & McCabe 1999), and the mean absolute relative error (MARE; Ji 2008):
(11)
(12)
(13)
(14)
where Yo and Ym are the time-matched discrete observed and modeled variables of interest, respectively. Overbars in Equations (11) and (12) represent time averages. The MIA statistic was selected in lieu of the Nash–Sutcliffe efficiency index (Nash & Sutcliffe 1970), as the latter tends to emphasize a lack of fit in peak regions and other rapidly changing parts of the time series curve (Legates & McCabe 1999).

For each calibration for each simulation, the intrinsic background vertical and horizontal turbulent mass transport parameters, Avo and Aho (Ji 2008) of the hydrodynamic and mass transport governing equations, and the fraction of shortwave solar radiation reaching the water surface (fr; -) and the vertical light extinction coefficient (λ; m−1), 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−6–0.01 m2 s−1; 1 × 10−4–0.025 m2 s−1; 0.45–0.90 (-); and 0.1–0.9 m−1, 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.

Residence time distribution analysis

A critical part of this work involved comparing observed and modeled RTD statistical features to assess the suitability of the EFDC model as a tool for evaluating the bulk hydraulic characteristics of constructed surface water flow-through wetlands during the basin design phase before construction or reconfiguration. To facilitate this comparison, the observed and modeled best-fit (calibrated) tracer response curves were converted to volume-based RTD functions using the transformation (Zuber 1986; Werner & Kadlec 1996):
(15)
where g(o,m) (t) is the dimensionless volume-based RTD function for the observed (o) or modeled (m) tracer response curves; V(o,m)(t) is the observed or modeled time-varying system volume (m3); C(o,m)(t) is the observed or modeled effluent tracer concentration (μg L−1). Because flow conditions varied during each tracer study and corresponding simulation, the volume-based RTD function was used in lieu of more standard steady-state transformations.
The modeled and observed mass recovery (kg), Mr(o,m), was estimated from:
(16)
where the last term in Equation (16) is the discrete time step approximation of the integral, and tf represents the time from injection (t = 0) at which the study was terminated. The time step, Δt, was set to 5 min for both observed and modeled time series.
All tracer response curve characteristics were evaluated against a dimensionless flow-weighted time ξ, used to represent the elapsed time from tracer injection on a common scale, and is defined as follows (Zuber 1986; Werner & Kadlec 1996):
(17)
where τ is a time-like dummy variable of integration, t is the time since injection, and Qo(o,m)(t) is the observed or modeled outflow discharge. The quantity within the integral and summation represents the turnover time and is the inverse of the instantaneous hydraulic residence time.
Dimensionless temporal features of each modeled and observed RTD, including the initial arrival time of tracer (ξi), peak arrival time (ξp), and median detention time (ξ50), were obtained directly by expressing ξ explicitly as a function of t (Zenger 2003; Green 2016):
(18)
where ξ is a monotonically increasing function of normal time. The modeled and observed flow-weighted time ξ was expressed as a cubic spline interpolant of the measured time from injection (t = 0), represented in standard time units. Spline interpolants were developed using the stats package in the R statistical computing environment (R Development Core Team 2008).

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.

Each observed and corresponding modeled RTD curve was analyzed using the method of moments to calculate primary RTD temporal statistics, including the centroid (the first moment about the origin; ), the temporal variance (the second moment about the centroid; ), and the normalized temporal variance (), defined, respectively, as follows (Kadlec 1994; Kadlec & Wallace 2008):
(19)
(20)
(21)
where, for simplicity, ξ in Equations (19) and (20) represents both observed and modeled estimates obtained from Equation (17).
The normalized variance of each observed and modeled RTD was used to approximate the observed and modeled system Péclet number, Pex(o,m), using (Levenspiel 2011):
(22)
The modeled and observed bulk rates of dispersion for each tracer response curve were estimated from the Péclet number and bulk approximations of basin morphometric characteristics using Kadlec & Wallace (2008):
(23)
where Kx(o,m) is the bulk quasi-longitudinal dispersion coefficient for the tracer study (m2 min−1), is the time average of the arithmetic mean of the modeled or observed inflow and outflow rates over the duration of each study, and , , and L are the average basin width (∼44 m), average basin depth (m), and the total length of the basin (∼327 m) along its centerline from the point of injection to the monitoring location, respectively.

Model calibration results

Time series plots of modeled, observed, and the minimum and maximum prediction bounds for all considered dynamic variables are given in Figures 24. Best-fit calibration parameters and fit statistics are given in Tables 3 and 4, respectively. The minimum and maximum prediction bounds in Figures 24 represent the prediction envelopes of modeled time-varying temperature, dye, and hydraulic variables over the tested parameter ranges in Table 4. These envelopes indicate the potential range of simulated state variables over the tested parameter ranges, which, in the case of the primary turbulence calibration parameters, the background horizontal and vertical turbulent viscosity, Aho and Avo, spanned between 4 and 5 orders of magnitude.
Table 3

Parameter values for each ‘best-fit’ calibrated simulation

StudyAvo (m2 s−1)Aho (m2 s−1)fr (-) (m day−1)Rc (day−1)
TS1 1.0 × 10–4 0.025 0.45 0.20 0.051 
TS2 1.0 × 10–4 0.01 0.90 0.45 0.008 
TS3 2.5 × 10–4 0.001 0.45 0.50 0.225 
TS4 1.0 × 10–6 0.01 0.90 0.45 0.425 
TS5 1.0 × 10–5 0.01 0.90 0.45 
TS6 7.5 × 10–4 0.0001 0.90 0.45 0.082 
StudyAvo (m2 s−1)Aho (m2 s−1)fr (-) (m day−1)Rc (day−1)
TS1 1.0 × 10–4 0.025 0.45 0.20 0.051 
TS2 1.0 × 10–4 0.01 0.90 0.45 0.008 
TS3 2.5 × 10–4 0.001 0.45 0.50 0.225 
TS4 1.0 × 10–6 0.01 0.90 0.45 0.425 
TS5 1.0 × 10–5 0.01 0.90 0.45 
TS6 7.5 × 10–4 0.0001 0.90 0.45 0.082 
Table 4

Simulation performance indices and statistics

StudyTS1TS2TS3TS4TS5TS6Mean
Water Surface Elevation (WSE) 
MIA 0.32 0.76 0.40 0.62 0.85 0.70 0.61 
MAE (m) 2 × 10−3 4 × 10−3 5 × 10−3 9 × 10−3 5 × 10−3 2.3 × 10−3 4.6 × 10−3 
MARE (%) 0.10 0.05 0.10 0.43 0.07 0.02 0.13 
RRMSE (%) 18.50 13.2 32.4 25.6 12.5 11.4 18.9 
ΔMax (m) 5 × 10−3 0.012 0.015 0.017 0.012 0.015 0.01 
Volume 
MIA 0.36 0.74 0.32 0.53 0.80 0.97 0.62 
MAE (m341.7 2.6 79.3 152 50.6 26.5 58.8 
MARE (%) 0.51 0.03 0.95 2.01 0.67 0.36 0.76 
RRMSE (%) 36.6 12.6 45.0 31.8 12.7 12.6 25.2 
ΔMax (m311.80 148.1 270.3 266.80 123.72 115.6 156.05 
Outflow Discharge 
MIA 0.86 0.96 0.95 0.95 0.97 0.97 0.94 
MAE (m3 s−11.0 × 10−5 6.0 × 10−4 1.8 × 10−4 4.0 × 10−4 2.0 × 10−3 2.8 × 10−4 5.8 × 10−4 
MARE (%) 0.05 1.09 0.35 0.60 3.07 0.95 1.02 
RRMSE (%) 4.8 3.1 2.6 3.4 3.2 2.1 3.2 
ΔMax (m3 s−13 × 10−3 0.01 1.8 × 10−3 0.014 0.025 7.1 × 10−3 0.01 
Outlet Temperature 
MIA 0.84 0.83 0.87 0.80 0.79 0.82 0.83 
MAE (oC) 0.65 0.07 0.76 0.25 0.76 0.42 0.49 
MARE (%) 11.44 0.63 4.17 4.31 4.35 2.98 4.65 
RRMSE (%) 9.5 6.8 9.3 9.4 11.0 15.1 10.2 
ΔMax (oC) 1.82 1.36 1.27 1.53 1.31 7.06 2.39 
Stilling Well Temperature 
MIA 0.78 0.79 0.86 0.77 0.71 0.67 0.76 
MAE (oC) 0.75 0.38 0.58 0.27 0.95 2.32 0.88 
MARE (%) 12.8 3.80 3.31 4.58 5.60 16.34 7.74 
RRMSE (%) 17.8 9.2 8.6 16.3 16.2 28.5 16.1 
ΔMax (oC) 3.03 0.64 3.38 0.75 1.94 8.74 3.08 
Tracer Concentration 
MIA 0.79 0.87 0.83 0.95 0.75 0.88 0.85 
MAE (μg L−10.39 0.22 0.10 0.002 0.21 0.02 0.16 
MARE (%) 6.21 4.25 4.16 0.30 15.6 1.85 5.40 
RRMSE (%) 11.7 8.9 13.1 4.0 9.5 8.1 9.2 
ΔMax (μg L−15.16 9.10 6.24 0.67 2.52 1.19 4.15 
StudyTS1TS2TS3TS4TS5TS6Mean
Water Surface Elevation (WSE) 
MIA 0.32 0.76 0.40 0.62 0.85 0.70 0.61 
MAE (m) 2 × 10−3 4 × 10−3 5 × 10−3 9 × 10−3 5 × 10−3 2.3 × 10−3 4.6 × 10−3 
MARE (%) 0.10 0.05 0.10 0.43 0.07 0.02 0.13 
RRMSE (%) 18.50 13.2 32.4 25.6 12.5 11.4 18.9 
ΔMax (m) 5 × 10−3 0.012 0.015 0.017 0.012 0.015 0.01 
Volume 
MIA 0.36 0.74 0.32 0.53 0.80 0.97 0.62 
MAE (m341.7 2.6 79.3 152 50.6 26.5 58.8 
MARE (%) 0.51 0.03 0.95 2.01 0.67 0.36 0.76 
RRMSE (%) 36.6 12.6 45.0 31.8 12.7 12.6 25.2 
ΔMax (m311.80 148.1 270.3 266.80 123.72 115.6 156.05 
Outflow Discharge 
MIA 0.86 0.96 0.95 0.95 0.97 0.97 0.94 
MAE (m3 s−11.0 × 10−5 6.0 × 10−4 1.8 × 10−4 4.0 × 10−4 2.0 × 10−3 2.8 × 10−4 5.8 × 10−4 
MARE (%) 0.05 1.09 0.35 0.60 3.07 0.95 1.02 
RRMSE (%) 4.8 3.1 2.6 3.4 3.2 2.1 3.2 
ΔMax (m3 s−13 × 10−3 0.01 1.8 × 10−3 0.014 0.025 7.1 × 10−3 0.01 
Outlet Temperature 
MIA 0.84 0.83 0.87 0.80 0.79 0.82 0.83 
MAE (oC) 0.65 0.07 0.76 0.25 0.76 0.42 0.49 
MARE (%) 11.44 0.63 4.17 4.31 4.35 2.98 4.65 
RRMSE (%) 9.5 6.8 9.3 9.4 11.0 15.1 10.2 
ΔMax (oC) 1.82 1.36 1.27 1.53 1.31 7.06 2.39 
Stilling Well Temperature 
MIA 0.78 0.79 0.86 0.77 0.71 0.67 0.76 
MAE (oC) 0.75 0.38 0.58 0.27 0.95 2.32 0.88 
MARE (%) 12.8 3.80 3.31 4.58 5.60 16.34 7.74 
RRMSE (%) 17.8 9.2 8.6 16.3 16.2 28.5 16.1 
ΔMax (oC) 3.03 0.64 3.38 0.75 1.94 8.74 3.08 
Tracer Concentration 
MIA 0.79 0.87 0.83 0.95 0.75 0.88 0.85 
MAE (μg L−10.39 0.22 0.10 0.002 0.21 0.02 0.16 
MARE (%) 6.21 4.25 4.16 0.30 15.6 1.85 5.40 
RRMSE (%) 11.7 8.9 13.1 4.0 9.5 8.1 9.2 
ΔMax (μg L−15.16 9.10 6.24 0.67 2.52 1.19 4.15 

Note. ΔMax represents the maximum difference between modeled and observed values.

Figure 2

Time series plots illustrating the correspondence between the modeled and observed state variables for the studies TS1 and TS2 conducted over the periods 2009-11-22 through 2009-12-01 and 2010-05-06 through 2010-05-13, respectively. For each panel, the solid black and thick dashed lines represent the observed and modeled variables, respectively. The lighter dashed lines show the maximum–minimum prediction envelope for each variable. Initial tracer injection times are shown as thick vertical dashed lines in the bottom panels labeled t0.

Figure 2

Time series plots illustrating the correspondence between the modeled and observed state variables for the studies TS1 and TS2 conducted over the periods 2009-11-22 through 2009-12-01 and 2010-05-06 through 2010-05-13, respectively. For each panel, the solid black and thick dashed lines represent the observed and modeled variables, respectively. The lighter dashed lines show the maximum–minimum prediction envelope for each variable. Initial tracer injection times are shown as thick vertical dashed lines in the bottom panels labeled t0.

Close modal
Figure 3

Time series plots illustrating the correspondence between the modeled and observed state variables for the studies TS3 and TS4 conducted over the periods 2010-05-18 through 2010-05-25 (left panel) and 2011-04-18 through 2011-04-21 (right panel), respectively.

Figure 3

Time series plots illustrating the correspondence between the modeled and observed state variables for the studies TS3 and TS4 conducted over the periods 2010-05-18 through 2010-05-25 (left panel) and 2011-04-18 through 2011-04-21 (right panel), respectively.

Close modal
Figure 4

Time series plots illustrating the correspondence between the modeled and observed state variables for the studies TS5 and TS6 conducted over the periods 2011-05-22 through 2011-05-25 (left panel) and 2012-04-23 through 2012-05-03 (right panel), respectively.

Figure 4

Time series plots illustrating the correspondence between the modeled and observed state variables for the studies TS5 and TS6 conducted over the periods 2011-05-22 through 2011-05-25 (left panel) and 2012-04-23 through 2012-05-03 (right panel), respectively.

Close modal

WSE, volume, and discharge

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 m3) and maximum difference estimates (11.8–266.8 m3) suggest reasonably accurate dynamic volume estimates over all simulations. Both modeled depths and volumes were highly sensitive to the value of Avo, but were not affected by Aho (data not shown). In general, higher values of Avo resulted in greater depths and larger modeled volumes. The maximum–minimum prediction bounds for these modeled time-varying variables showed that unique Aho and Avo 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.

Temperatures and tracer concentrations

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 24.

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 Avo, and relatively insensitive to Aho, fr, and λ. Higher values of Avo, 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, Aho 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 fr and λ influenced the magnitude of modeled temperature curves but did not influence diurnal patterns at both monitoring locations. Of these two parameters, fr maintained a stronger influence on model results. In general, decreases in fr, 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−1, 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−1 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−1 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 Avo and Aho. Optimized values of Aho ranged from 2.5 × 10−2 to 1 × 10−4 m2 s−1. In general, larger values of Aho induced earlier initial arrival times and longer and smoother response curve tails. Optimized values of Avo ranged from 2.5 × 10−4 to 1 × 10−6 m2 s−1, with smaller values resulting in earlier tracer initial and peak arrival times and sharper ascending limbs. That different combinations of Avo and Aho 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 Avo and Aho tested in this work is illustrated in the minimum-maximum bounds of modeled dye concentrations for each simulation period, as shown in Figures 24. These bounds also show the relative insensitivity of response curve tails to the magnitudes of Avo and Aho. 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.

Residence time distribution analysis results

Modeled and observed RTD moments and temporal features were in reasonable agreement (Figure 5; Green 2016). For most studies, the dimensionless initial time of tracer arrival, ξi, was reasonably replicated, with minimum and maximum percent differences between modeled and observed values of 0.07 and 33.2%, respectively. Modeled dimensionless peak arrival times generally fell to within 10% of observed ξp, except for TS5, which peaked 50% earlier. Modeled values of ξp were observed to consistently be under-predicted with increasing values of observed counterparts (Figure 5), suggesting that the EFDC model of this system tends to over-predict internal velocities under lower flow rates. Overall, modeled and observed median and mean dimensionless residence times and normalized variances were in reasonably good agreement, showing little to no bias, with percent differences ranging from 0.74 to 22% and 1.4 to 35%, respectively. This pattern is also observed for estimates of modeled and observed Péclet numbers and longitudinal dispersion rates. The exception again is TS5, for which the simulated values for ξ50, , , and Kx were significantly lower than observed values.
Figure 5

Comparison of observed (o) and modeled (m) calculated RTD characteristics and mixing parameters. One-to-one lines are represented as dashed lines.

Figure 5

Comparison of observed (o) and modeled (m) calculated RTD characteristics and mixing parameters. One-to-one lines are represented as dashed lines.

Close modal

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 Kx. In our opinion, the exclusion of vegetation in the model structure for TS5 is the likeliest explanation for the significant under-prediction of Kx for this study.

Assessment of environmental effects on RTD features

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).

Table 5

Tested sensitivity analysis scenarios

ScenarioApplied environmental effects
Applied wind series; Steady flow equal to the flow at dye injection;
Applied atmospheric radiative forcing 
No applied wind series; Steady flow equal to the flow at dye injection;
Applied atmospheric radiative forcing 
No applied wind series; Steady flow equal to the flow at dye injection;
No applied atmospheric radiative forcing 
No applied wind series; Dynamic flow boundary conditions; Applied atmospheric radiative forcing 
No applied wind series; Dynamic flow boundary conditions;
No applied atmospheric radiative forcing 
Applied wind series; Dynamic flow boundary conditions; No applied atmospheric radiative forcing 
Applied wind forcing; Steady flow equal to the flow at dye injection;
No applied atmospheric radiative forcing 
ScenarioApplied environmental effects
Applied wind series; Steady flow equal to the flow at dye injection;
Applied atmospheric radiative forcing 
No applied wind series; Steady flow equal to the flow at dye injection;
Applied atmospheric radiative forcing 
No applied wind series; Steady flow equal to the flow at dye injection;
No applied atmospheric radiative forcing 
No applied wind series; Dynamic flow boundary conditions; Applied atmospheric radiative forcing 
No applied wind series; Dynamic flow boundary conditions;
No applied atmospheric radiative forcing 
Applied wind series; Dynamic flow boundary conditions; No applied atmospheric radiative forcing 
Applied wind forcing; Steady flow equal to the flow at dye injection;
No applied atmospheric radiative forcing 

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.

Sensitivity analysis results

As shown in Figure 6, wind and thermal forcing, separately and in concert, can strongly influence considered short-circuiting indices (ξi and ξp) and some measures of bulk basin mixing (, Kx, and Pex). The effect of wind on RTD development is evident for scenarios 2 and 4. In these scenarios, the absence of wind forcing with applied atmospheric forcing is shown to have the strongest influence over the initial tracer arrival times, resulting in average increases of these RTD characteristics of 112 and 116%, corresponding to minimum and maximum increases in ξi of 11–212%, and 13–242% for steady and unsteady flow, respectively. A similar effect is observed for ξp, with average increases of this statistic of 73 and 93%, corresponding to increases above the baseline condition of 0–150% for steady flow and 1–235% for unsteady flow. In the absence of wind alone, dimensionless variances and estimated longitudinal dispersion coefficients decreased by 12–18% and 0–11%, respectively, with responses ranging from +14 to −52%, and +57 to −84% with only nominal differences between transient and steady flow conditions. These observed changes in and Kx corresponded to an increase of 16 and 25% for in average system Péclet numbers for transient and steady flow conditions, with responses ranging from +84 to −6%, and +53 to −10%, respectively. In general, when considered independently, wind forcing was observed to have a minor influence on the magnitude of the median and mean detention times.
Figure 6

Box plots showing the percent deviation from the baseline best-fit simulation for each considered RTD characteristic and derived mixing metric for each tested scenario.

Figure 6

Box plots showing the percent deviation from the baseline best-fit simulation for each considered RTD characteristic and derived mixing metric for each tested scenario.

Close modal

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 ξ50 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 ξ50 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 ξ50 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 Kx 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 Kx 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 Kx 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 (Avo and Aho) 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 Avo and Aho 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.

Badrot-Nico
F.
,
Guinot
V.
&
Brissaud
F.
(
2009
)
Fluid flow pattern and water residence time in waste stabilization ponds
,
Water Science & Technology
,
59
(
6
),
1061
.
doi:10.2166/wst.2009.087
.
Bentzen
T. R.
,
Larsen
T.
&
Rasmussen
M. R.
(
2008
)
Wind effects on retention time in highway ponds
,
Water Science & Technology
,
57
(
11
),
1713
.
doi:10.2166/wst.2008.267
.
Buchanan
T. J.
&
Somers
W. P.
(
1969
)
Discharge measurements at gaging stations
. In:
Techniques of Water-Resources Investigations of the United States Geological Survey
.
USGS Numbered Series 03-A8
.
U.S. Govt. Print Office
,
Cole
T. M.
&
Wells
S. A.
(
2003
)
CE-QUAL-W2: A two-Dimensional, Laterally Averaged, Hydrodynamic and Water Quality Model, Version 3.1. Instruction Report EL-03-1
.
Vicksburg, MS, USA
:
US Army Engineering and Research Development Center
.
Crumpton
W. G.
,
Stenback
G. A.
,
Fisher
S. W.
,
Stenback
J. Z.
&
Green
D. I. S.
(
2020
)
Water quality performance of wetlands receiving nonpoint-source nitrogen loads: nitrate and total nitrogen removal efficiency and controlling factors
,
Journal of Environmental Quality
,
49
,
1
10
.
https://doi.org/10.1002/jeq2.20061
.
Devkota
J.
,
Fang
X.
&
Fang
V. Z.
(
2013
)
Response characteristics of the Perdido and Wolf Bay system to inflows and sea level rise
,
International Journal of Environment and Climate Change
,
3
(
2
),
229
256
.
doi:10.9734/BJECC/2013/3516
.
Eeling
J.
(
2019
)
Temporal and Spatial Dynamics of Submersed Macrophytes and Water Temperatures in Shallow Flow-Through Wetlands
.
Master's thesis
,
Ames, IA, USA
:
Iowa State University
.
Fischer
H. B.
,
List
J. E.
,
Koh
C. R.
,
Imberger
J.
&
Brooks
N. H.
(
1979
)
Mixing in Natural Waters
.
San Diego, CA, USA
:
Academic Press
.
Franceschini
S.
&
Tsai
C. W.
(
2010
)
Assessment of uncertainty sources in water quality modeling in the Niagara river
,
Advances in Water Resources
,
33
,
493
503
.
doi:10.1016/j.advwatres.2010.02.001
.
Green
D. I.
(
2016
)
Environmental Factors Affecting the Residence Time Distribution Dynamics of Constructed Agricultural Wetlands
.
PhD dissertation
,
Ames, IA, USA
:
Iowa State University
.
Gu
R.
,
Luck
F. N.
&
Stefan
H. G.
(
1996
)
Water quality stratification in shallow wastewater stabilization ponds
,
Journal of the American Water Resources Association
,
32
(
4
),
831
834
.
Hamrick
J. M.
(
1992
)
A Three-Dimensional Environmental Fluid Dynamics Computer Code: Theoretical and Computational Aspects
.
Gloucester Point, VA, USA
:
The College of William and Mary, Virginia Institute of Marine Science
,
Special Report 317
, pp.
63
.
Hamrick
J. M.
,
Wu
T. S.
, (
1997
)
Computational design and optimization of the EFDC/HEM3D surface water hydrodynamic and eutrophication models
. In:
Delich
G.
&
Wheeler
M. F.
(eds.)
Next Generation Environmental Models and Computational Methods
,
Philadelphia
Society of Industrial and Applied Mathematics
, pp.
143
156
.
Huang
W.
,
Liu
X.
&
Chen
X.
(
2008
)
Numerical modeling of hydrodynamics and salinity transport in little manatee river
,
Journal of Coastal Research
,
2008
,
13
24
.
doi:10.2112/1551-5036-52.sp1.13
.
Irwin
J. S.
(
1979
)
A theoretical variation of the wind profile power-law exponent as a function of surface roughness and stability
,
Atmospheric Environment
,
13
(
1
),
191
194
.
doi:10.1016/0004-6981(79)90260-9
.
Jenkins
G. A.
&
Greenway
M.
(
2005
)
The hydraulic efficiency of fringing versus banded vegetation in constructed wetlands
,
Ecological Engineering
,
25
(
1
),
61
72
.
doi:10.1016/j.ecoleng.2005.03.001
.
Ji
Z.-G.
(
2008
)
Hydrodynamics and Water Quality: Modeling Rivers, Lakes, and Estuaries
.
West Sussex, England
:
John Wiley & Sons
.
Jin
K.-R.
&
Ji
Z.-G.
(
2013
)
A long term calibration and verification of a submerged aquatic vegetation model for Lake Okeechobee
,
Ecological Processes
,
2
(
1
),
1
13
.
doi:10.1186/2192-1709-2-23
.
Jin
K.-R.
,
Hamrick
J. H.
&
Tisdale
T.
(
2000
)
Application of three-dimensional hydrodynamic model for Lake Okeechobee
,
Journal of Hydraulic Engineering: (ASCE)
,
https://doi.org/10.1061/(ASCE)0733-9429(2000)126:10(758)
.
Jin
K.-R.
,
Ji
Z.-G.
&
Hamrick
J. H.
(
2002
)
Modeling winter circulation in Lake Okeechobee, Florida
,
Journal of Waterway, Port, Coastal, and Ocean Engineering
,
128
(
3
),
114
125
.
doi:10.1061/(ASCE)0733-950X(2002)128:3(114)
.
Jin
K.-R.
,
Ji
Z.-G.
&
James
R. T.
(
2007
)
Three-dimensional water quality and SAV modeling of a large shallow lake
,
Journal of Great Lakes Research
,
33
(
1
),
28
45
.
doi:10.3394/0380-1330(2007)33[28:TWQASM]2.0.CO:2
.
Kadlec
R. H.
(
1994
)
Detention and mixing in free water wetlands
,
Ecological Engineering
,
3
(
4
),
345
380
.
doi:10.1016/0925-8574(94)00007-7
.
Kadlec
R. H.
&
Wallace
S.
(
2008
)
Treatment Wetlands
, 2nd edn.
Boca Raton, FL, USA
:
CRC Press
.
Legates
D. R.
&
McCabe
G. J.
(
1999
)
Evaluating the use of ‘goodness-of-fit’ measures in hydrologic and hydroclimatic model validation
,
Water Resources Research
,
35
(
1
),
233
241
.
doi:10.1029/1998WR900018
.
Levenspiel
O.
(
2011
)
Tracer Technology: Modeling the Flow of Fluids
.
Berlin, Germany
:
Springer Science & Business Media
.
Li
Y.
,
Acharya
K.
,
Chen
D.
&
Stone
M.
(
2010
)
Modeling water ages and thermal structure of Lake Mead under changing water levels
,
Lake and Reservoir Management
,
26
(
4
),
258
272
.
doi:10.1080/07438141.2010.541326
.
Macdonald
R. J.
&
Ernst
A.
(
1986
)
Disinfection efficiency and problems associated with maturation ponds
,
Water Science and Technology
,
18
(
10
),
19
29
.
Nash
J. E.
&
Sutcliffe
J. V.
(
1970
)
River flow forecasting through conceptual models part I – a discussion of principles
,
Journal of Hydrology
,
10
(
3
),
282
290
.
Pebesma
E. J.
(
2004
)
Multivariable geostatistics in S: the gstat package
,
Computers & Geosciences
,
30
(
7
),
683
691
.
doi:10.1016/j.cageo.2004.03.012
.
Pedahzur
R.
,
Nasser
A. M.
,
Dor
I.
,
Fattal
B.
&
Shuval
H. I.
(
1993
)
The effect of baffle installation on the performance of a single-cell stabilization pond
,
Water Science and Technology
,
27
(
7–8
),
45
52
.
Persson
J.
(
2000
)
The hydraulic performance of ponds of various layouts
,
Urban Water
,
2
(
3
),
243
250
.
doi:10.1016/S1462-0758(00)00059-5
.
Persson
J.
,
Somes
N.
&
Wong
T.
(
1999
)
Hydraulics efficiency of constructed wetlands and ponds
,
Water Science and Technology
,
40
(
3
),
291
289
.
Pilgrim
J. M.
,
Fang
X.
&
Stefan
H. G.
(
1998
)
Stream temperature correlations with air temperatures in Minnesota: implications for climate warming
,
Journal of the American Water Resources Association
,
34
(
5
),
1109
1121
.
doi:10.1111/j.1752-1688.1998.tb04158.x
.
R Development Core Team
(
2008
)
R: a Language and Environment for Statistical Computing
.
Vienna, Austria
:
R Foundation for Statistical Computing
.
Smart
P. L.
&
Laidlaw
I. M. S.
(
1977
)
An evaluation of some fluorescent dyes for water tracing
,
Water Resources Research
,
13
(
1
),
15
33
.
doi:10.1029/WR013i001p00015
.
Sweeney
D. G.
,
Nixon
J. B.
,
Cromar
N. J.
&
Fallowfield
H. J.
(
2005
)
Profiling and modelling of thermal changes in a large waste stabilization pond
,
Water Science and Technology
,
51
(
12
),
163
172
.
Thackston
E. L.
,
Shields
D.
Jr.
&
Schroeder
P. R.
(
1987
)
Residence time distributions of shallow basins
,
Journal of Environmental Engineering
,
113
(
6
),
1319
1332
.
doi:10.1061/(ASCE)0733-9372(1987)113:6(1319)
.
Watters
G. Z.
,
Mangelson
K. A.
&
George
R. L.
(
1973
)
The Hydraulics of Waste Stabilization Ponds
.
Final Report to the Office of Water Resources Research on Project A-008-Utah
.
Logan, UT, USA
:
Utah State University
.
Werner
T. M.
&
Kadlec
R. H.
(
1996
)
Application of residence time distributions to stormwater treatment systems
,
Ecological Engineering
,
7
(
3
),
213
234
.
doi:10.1016/0925-8574(96)00013-4
.
Willmott
C. J.
,
Ackleson
S. G.
,
Davis
R. E.
,
Feddema
J. J.
,
Klink
K. M.
,
Legates
D. R.
,
O'Donnell
J.
&
Rowe
C. M.
(
1985
)
Statistics for the evaluation and comparison of models
,
Journal of Geophysical Research: Oceans
,
90
(
C5
),
8995
9005
.
doi:10.1029/JC090iC05p08995
.
Wu
J.
(
1969
)
An estimation of wind effects on dispersion in wide channels
,
Water Resources Research
,
5
(
5
),
1097
1104
.
doi:10.1029/WR005i005p01097
.
Wu
G.
&
Xu
Z.
(
2011
)
Prediction of algal blooming using EFDC model: case study in the Daoxiang Lake
,
Ecological Modelling
,
222
(
6
),
1245
1252
.
doi:10.1016/j.ecolmodel.2010.12.021
.
Zenger
K.
(
2003
)
Modelling, Analysis and Controller Design of Time-Variable Flow Processes
.
PhD dissertation
,
Espoo, FInland
:
Helsinki University of Technology
.
Zhang
C.
,
Gao
X.
,
Wang
L.
&
Chen
Y.
(
2013
)
Analysis of agricultural pollution by flood flow impact on water quality in a reservoir using a three-dimensional water quality model
,
Journal of Hydroinformatics
,
15
(
4
),
1061
1072
.
doi:10.2166/hydro.2012.131
.
Zoppou
C.
(
1999
)
Reverse routing of flood hydrographs using level pool routing
,
Journal of Hydrologic Engineering
,
4
(
2
),
184
188
.
Zuber
A.
(
1986
)
On the interpretation of tracer data in variable flow systems
,
Journal of Hydrology
,
86
(
1–2
),
45
57
.
doi:10.1016/0022-1694(86)900
.
This is an Open Access article distributed under the terms of the Creative Commons Attribution Licence (CC BY-NC 4.0), which permits copying, adaptation and redistribution for non-commercial purposes, provided the original work is properly cited (http://creativecommons.org/licenses/by-nc/4.0/).