Home Introduction Site Description Tracer Design & Methodology >Modeling Results Discussion Acknowledgments References Figures, Table & Equations PDF Version

4.1. Analysis Equations

[11] An extended version of the one-dimensional solute transport model by Runkel [1998], developed by Choi et al. [2000], was used to characterize the results of tracer experiments. We refer to our modified version of OTIS-P as the OTIS-2stor module. The governing equations are

  D

  D

  D

where

C main channel solute concentration [mg L^-1];
Q volumetric flow rate [m³ s^-1];
A main flow zone cross-sectional area [m²];
D_L longitudinal dispersion coefficient [m² s^-1];
q_Lⁱⁿ lateral inflow rate [m³ s^-1 m^-1];
q_L^out lateral outflow rate [m³ s^-1 m^-1];
C_L lateral inflow solute concentration [mg L^-1];
₁ storage zone 1 exchange coefficient [s^-1];
C_S1 storage zone 1 solute concentration [mg L^-1];
A_S1 storage zone 1 cross-sectional area [m²];
₂ storage zone 2 exchange coefficient [s^-1];
C_S2 storage zone 2 solute concentration [mg L^-1];
A_S2 storage zone 2 cross-sectional area [m²];
t time [s];
x distance [m].

[12] The governing equations of the OTIS-2stor module were solved by modifying the USGS numerical code OTIS-P (One-dimensional Transport with Inflow and Storage with Parameter Estimation) [Runkel, 1998] to solve the extended governing equations for a flow system with two-storage zones. Since the model extension is not the primary contribution of the present paper, the conceptual basis of the model is only briefly reviewed here. The reader is referred to Choi et al. [2000], Harvey and Wagner [2000], and Runkel [1998] for more detailed discussions of the model's conceptual basis, its relation to earlier models, and solution techniques. A thorough documentation of the OTIS-2stor module is in preparation.

[13] Both the OTIS-2stor module and its predecessor are often used for the analysis of tracer experiments by applying them in the inverses sense, i.e., model parameters are adjusted using statistical optimization to determine the values of the parameters that best simulate the measured tracer data. Average velocity and cross sectional area of flow are two of the primary parameters of interest. Both a mean velocity, V, and an associated cross sectional area, A, are identified through optimization of the parameters of the model. Typically the modeling shows that not all parts of the water column participate in downstream advection, due to the existence of zones with stagnant or very slowly flowing water referred to as "storage" zones. The zone with cross sectional area A is therefore typically smaller than measurements of the total cross sectional area of surface water. We therefore refer the modeled zone where downstream advection occurs as the "main flow zone" in order to distinguish it from the total cross sectional area.

[14] In addition to quantifying advection we were also interested in quantifying parameters that describe longitudinal mixing of solutes. The first of those is the longitudinal dispersion parameter, D_L, which characterizes the relatively fast components of mixing that arise due to velocity variations in the main flow zone. To be characterized by the longitudinal dispersion parameter, mixing processes must achieve equilibrium (i.e., all tracer must experience the full range of flow velocities in the main flow zone) during the tracer experiment. Slower processes of longitudinal mixing are better characterized as "storage exchange." Conceptually, storage is the result of water being exchanged between the main flow zone and "storage zones," where water is stagnant or very slow moving relative to the main flow zone. The storage parameters and A_s characterize these slower components of mixing that do not achieve a state of equilibrium mixing during the experiment [Harvey and Wagner, 2000]. The residence time distribution for fluid and solute entering the storage zone is exponential, with an average residence time equal to A_s/(A). An additional parameter of interest is the depth-averaged velocity of tracer, V', which is slower than velocity in the main flow zone. Depth-averaged velocity can be estimated from the model-estimated parameters as V x A/(A + A_s) where V is the velocity for the main flow zone, while A and A_s are the cross sectional areas of the main flow zone and storage zone, respectively.

[15] Using two storage zones instead of just one as prescribed by Choi et al. [2000] allows a broader range of storage timescales to be characterized, due to the model's inclusion of a second storage zone with a different (usually a much longer) mean residence time. Therefore, for the twostorage model there are four storage parameters (₁, A_S1, ₂, and A_S2) that must be estimated to characterize exchange between the main flow zone and storage environments. Each storage zone has a unique fluid residence time, A_S1/(₁A) and A_S2/(₂A). Also, V' is now calculated as V x A/(A + A_S1 + A_S2) where A_S1 and A_S2 are the cross sectional areas of storage zones 1 and 2, respectively.

4.2. Approach for Estimating Model Parameters

[16] Our purpose in using the two-storage-zone model was to test the hypothesis that the major mechanisms of storage in the Everglades wetlands could be identified and quantified through tracer experimentation and inverse modeling. We anticipated that the mixing that resulted from velocity variations within the macrophyte stems would be characterized by the longitudinal dispersion parameter, while mixing that resulted by exchange between the main flow zone and zones of much more slowly moving water in floating vegetation and in peat pore water could be characterized by the two storage zones, respectively.

[17] The volumetric flow rate of surface water in the channel (Q) and inflow and outflow fluxes (q_Lⁱⁿ, q_L^out) were calculated using estimates of (1) injection pumping rate, (2) Br concentration in the injection solution barrel, (3) background concentration of Br tracer, and (4) concentrations of Br in surface water during the injection a short distance downstream of the injection and at a distance of 6.8 m. All of the other parameters were estimated by the nonlinear least squares optimization described previously. The optimized parameters included, depending on the particular simulation, A, D_L, A_S1, and A_S2, and exchange coefficients, ₁ and ₂. Model runs began at the start of the injection, but, due to practical and theoretical reasons, only tracer measurements during the last 4 hours of the tracer plateau and the 27 hours following the end of the injection (i.e., the tail of the breakthrough curve) were used for parameter calibration. The reason was that the background tracer concentration is known with considerably more certainty compared with the plateau tracer concentration. As a consequence, the task of estimating storage parameters on the tail of the experiment (compared with on the rising limb of the experiment) is more reliable. Note that this argument only applies if a plateau tracer concentration is reached during the experiment, because that condition allows volumetric flow rate and inflow to be calculated independently of the model optimization. If volumetric flow rate and inflow are not known from plateau tracer data (as is typical for a "pulse" tracer breakthrough that never achieves plateau), then all the data from the rising and falling limb of the tracer experiment must be used in the statistical optimization.

4.3. Optimization and Uncertainty of Estimated Parameters

[18] The "best fit" values of the transport parameters were estimated by an inverse approach that uses generalized nonlinear least squares, together with other statistical criteria, to objectively search for a set of parameters that minimize the differences between model calculations and observations. An optimization routine that has frequently been used in the past for that purpose is called STARPAC [Donaldson and Tryon, 1990]. Its use with stream tracer applications is thoroughly documented by Wagner and Harvey [1997], Runkel [1998], and Harvey and Wagner [2000]. The reliability of the parameter values that result from optimization are a function of several criteria. First is the choice of the objective function itself (along with other quantitative criteria that help determine whether convergence on a solution has occurred). STARPAC uses a formulation of the residual sum of squares (RSS) which can be weighted (to account for unequal variances associated with observations). This weighting scheme treats the problem of unequal variances of observations across the full range of values of the tracer observations, which helps to emphasize the valuable information about storage parameters contained within the lowest magnitude concentrations on the tail of the breakthrough curve [Wagner and Harvey, 1997]. To make the test for convergence even more rigorous, we always use the established convention of changing the values of the starting parameters and rerunning even the models that successfully converged to ensure that resulting parameter estimates are not affected by choice of starting parameter values. A second important indicator of the reliability of parameters determined by inverse modeling is the standard deviations (reported here as coefficients of variation) that are associated with the estimation of the "optimal" parameter values. These uncertainties are a function of the sensitivity of the model output with respect to changes in the parameter values, as well as the number of data collected and assumptions about the precision of those data [Wagner and Harvey, 1997]. In our experience, if all coefficients of variation are small relative to the parameter estimates (we usually judge them favorably if CVs are less than 0.5), the parameters can be considered to be reliably estimated. A final method of judging reliability of estimated parameters is from a calculation of the experimental Damkohler number, Da1 = ((1 + A_S/A)L)/V, where L, the length of the experimental reach, is the only parameter not previously defined. The Damkohler number expresses the relative importance of different solute transport mechanisms (i.e., downstream solute advection compared with solute mass transfer into and out of storage zones). Wagner and Harvey [1997] demonstrated that parameter reliability was greatest (i.e., uncertainty was lowest) in situations where Da1 was between 0.01 and 100. Da1 values much less than 0.01 indicate situations where relatively little of the tracer mass interacted with storage zones compared to what was advected downstream. Sensitivity to determine storage parameters in that situation is usually limited. In contrast, Da1 values much greater than 100 indicate situations where it is likely that too much interaction occurred between the tracer and storage zones. In those situations the storage parameters are likely to be unidentifiable because of nonuniqueness issues (i.e., tracer data can be simulated equally well either by adjusting the longitudinal dispersion coefficient alone, storage parameters alone, or all of these parameters simultaneously).

4.4. Comparison of Modeling Efficiency Using Zero, One, and Two Storage Zones

[19] To further assess the validity of using a tracer model with storage zones, the results obtained using the OTIS-2stor module (with two storage zones) were compared with results obtained using the original OTIS model (using first one storage zone and then zero storage zones). As a means to compare the efficiency of the three models in describing the tracer data, we calculated the model selection criterion (MSC), which is based on the Akaike information criterion [Akaike, 1974]. The MSC computes the fraction of the total variance explained by the model but applies a penalty based on the number of optimized parameters. The model with the highest MSC is generally considered to be the most appropriate for describing the experimental data because of the higher information content [Koeppenkastrop and DeCarlo, 1993]. The MSC was calculated as,

  D

where i indicates the ith observation of tracer concentration, C_i and c_i are the observed and simulated tracer concentrations, respectively, at the ith observation time, w_i is the weight factor, p is the number of optimized parameters, and m is the total number of observations used in the optimization.