Modeling surface and ground water mixing in the hyporheic zone using MODFLOW and MT3D

Advances in Water Resources 29 (2006) 1618–1633 www.elsevier.com/locate/advwatres

Modeling surface and ground water mixing in the hyporheic zone using MODFLOW and MT3D Laura K. Lautz a

a,*

, Donald I. Siegel

b,1

Department of Forest and Natural Resources Management, 206 Marshall Hall, One Forestry Drive, SUNY, College of Environmental Science and Forestry, Syracuse, NY 13210-2778, USA b Department of Earth Sciences, Syracuse University, 204 Heroy Geology Lab, Syracuse, NY 13244, USA Received 5 August 2005; received in revised form 1 November 2005; accepted 2 December 2005 Available online 2 February 2006

Abstract We used a three-dimensional MODFLOW model, paired with MT3D, to simulate hyporheic zones around debris dams and meanders along a semi-arid stream. MT3D simulates both advective transport and sink/source mixing of solutes, in contrast to particle tracking (e.g. MODPATH), which only considers advection. We delineated the hydrochemically active hyporheic zone based on a new definition, specifically as near-stream subsurface zones receiving a minimum of 10% surface water within a 10-day travel time. Modeling results indicate that movement of surface water into the hyporheic zone is predominantly an advective process. We show that debris dams are a key driver of surface water into the subsurface along the experimental reach, causing the largest flux rates of water across the streambed and creating hyporheic zones with up to twice the cross-sectional area of other hyporheic zones. Hyporheic exchange was also found in highly sinuous segments of the experimental reach, but flux rates are lower and the cross-sectional areas of these zones are generally smaller. Our modeling approach simulated surface and ground water mixing in the hyporheic zone, and thus provides numerical approximations that are more comparable to field-based observations of surface–groundwater exchange than standard particle-tracking simulations.  2005 Elsevier Ltd. All rights reserved. Keywords: Hyporheic zone; MODFLOW; MT3D; Groundwater flow model; Debris dam

1. Introduction The hyporheic zone lies at the interface between surface and ground waters, in near-stream sediments. Hyporheic flow paths, centimeters to tens of meters in length, are nested within the groundwater flow system and begin in the stream, returning to the channel within the reach of interest [10,11]. Hyporheic exchange increases the contact stream water solutes have with subsurface sediments and microbial communities. By redirecting water from the stream to the near-stream environment, hyporheic mixing moves solutes to alternating oxic and anoxic environments, *

Corresponding author. Tel.: +1 315 470 4765. E-mail addresses: [email protected] (L.K. Lautz), [email protected] (D.I. Siegel). 1 Tel.: +1 315 443 3607. 0309-1708/$ - see front matter  2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.advwatres.2005.12.003

geochemically active sediments, and microorganisms, and increases the potential for nutrient retention and enhanced microbial activity [5]. It is important to identify the most important characteristics of streams and their adjacent subsurfaces that drive hyporheic exchange. With this information the degree of hyporheic interaction in a variety of geomorphic and hydrologic settings, and at difference scales, can be compared and contrasted. Several hydrologic and geomorphic factors have been shown to influence hyporheic interaction, including debris dams [15], stream meanders [29], variability of streambed topography along pool-riffle sequences [10], the hydraulic conductivity of stream sediments [18], and stream water velocity [4]. Much of the current understanding of hyporheic exchange comes from empirical observations of what effect hyporheic interactions have on solute transport in experi-

L.K. Lautz, D.I. Siegel / Advances in Water Resources 29 (2006) 1618–1633

mental reaches [2,4,18]. In field experiments, tracers are injected into the reach of interest and the arrival of the solute is monitored downstream. The shape of the tracer breakthrough curve is used to infer the degree of temporary hyporheic storage along the reach (i.e. [28]). Although in-stream tracer tests characterize the degree to which hyporheic interactions occur along individual, experimental reaches, they only capture the net effect of those interactions on stream solute transport. A hyporheic zone’s reach-average residence time, size or flux rate, determined from these tracer experiments, does not indicate actual values for any particular point along the stream. Furthermore, in-stream tracer tests do not directly identify where hyporheic exchange takes place or what drives the exchange. Because results of in-stream tracer tests are empirical, rather than process-based, it can be difficult to generalize results to other reaches or other hydrologic or geomorphic conditions. Computer simulation experiments of groundwater flow are a way to bridge the gap between field observations and general characterization of stream hyporheic zone processes. Groundwater flow models simulate the spatial variability of movement of surface water into the subsurface from a process-based perspective, in contrast to inverse modeling of lumped exchange parameters from in-stream tracer tests [20]. The results of such models, which simulate processes that drive hyporheic exchange, need to be paired with empirical observations of the net effect of hyporheic interactions along experimental stream reaches. Although strongly affected by stream channel characteristics, hyporheic interactions can be simulated using numerical groundwater flow models because hyporheic flow paths are nested within the local groundwater flow system [10,11] and boundary conditions within the model can simulate driving forces in stream channels [10,29]. Groundwater flow models further allow simulation and quantification of the gross effect of hyporheic interactions, including where the exchange takes place, and the residence times and flux rates in discrete zones along reaches. Modeling results calibrated on data from well-instrumented sites can be generalized to other sites of interest and used as a heuristic and sometimes predictive tool. MODFLOW, a well-established US Geological Survey computer code that solves the groundwater flow equation [9,16], has been successfully used to simulate water flow through the hyporheic zone. For example, Wroblicky et al. [29] simulated the lateral extent of the hyporheic zone and the flux rates through the hyporheic zone in two dimensions along two first-order stream channels. Storey et al. [23] simulated subsurface flow within a single riffle along a gravel bed stream in Southern Ontario. Kasahara and Wondzell [14] estimated hyporheic fluxes and residence times along mountain streams in the Cascades and examined the relative influences of different geomorphic features there. These studies simulated water transport through the hyporheic zone using the MODPATH particle-tracking

1619

package [21], and the hyporheic zone was defined as those flow paths that exit and re-enter the stream channel [14,23] within 10 days or less [29]. MODPATH predicts the travel time of an advectively transported parcel of water along a flow path. Although MODPATH can be used to assess individual water parcel source areas and average travel times from those sources it cannot directly determine the percentages of water at a given point that are derived from multiple source areas. Nor can it reflect the mixing of surface and ground waters in the near-stream environment or the dispersion of dissolved solutes as they travel through the subsurface. Particle-tracking models that have been used to delineate the hyporheic zone have relied on a purely hydrological definition of the hyporheic zone, defined as those flow paths exiting and returning to the channel within the reach of interest. This definition is most useful when considering the impact of hyporheic exchange on the transport of solutes down fluvial channels, where full recovery of an instream tracer is possible only if all hyporheic subsurface flow paths return the tracer to the stream channel. Alternatively, when considering the impact of hyporheic exchange specifically on subsurface water chemistry, defining the hyporheic zone based solely on flow path ignores the importance of subsurface residence times with respect to solute transformations and oxygen depletion. The hyporheic zone is an important component of the stream ecosystem where relatively high concentrations of dissolved oxygen are sustained because of recharging stream water. These high oxygen concentrations support microbial and invertebrate communities, and when the oxygen is depleted, closely spaced oxic and anoxic environments that can promote nutrient cycling develop. To model the impact of hyporheic exchange on subsurface water chemistry, it may be more appropriate to delineate the hydrochemically active hyporheic zone using a new definition—specifically as locations where geochemically reactive surface water solutes are present in high concentrations within a designated travel time. A modeling approach that combines advective transport with dispersion and mixing of water from different sources is a more effective method to simulate mixing of surface and ground waters in the hyporheic zone than purely advective modeling. Indeed, dispersion and mixing have explicitly been used in field experiments to define hyporheic zones as those places where greater than 10% of subsurface water is derived directly from the stream (in contrast to ground water contributions) [25]. By using computer simulations of groundwater flow and dispersive solute transport to model surface and ground water mixing, field measurements of hyporheic zone chemistry can be directly compared to model results. This comparison is not possible using solely MODPATH simulations. Groundwater flow modeling, using MODFLOW, has been paired with solute transport simulations, using MT3D, to simulate the movement of surface water solutes into the subsurface. Woessner [27] showed how variable

1620

L.K. Lautz, D.I. Siegel / Advances in Water Resources 29 (2006) 1618–1633

streambed hydraulic conductivity and topography can influence mixing of surface and ground water in the streambed by constructing a generic hydrologic model of a sand and gravel streambed, assigning surface water boundaries a constant concentration of 100 mg/l and groundwater an initial concentration of 0 mg/l. He delineated the 10% surface water mixing zone based on the steady-state 10 mg/l isoline. Cardenas et al. [3] used the same technique to examine the impact of a sinusoidal stream stage pattern on hyporheic zone development in highly heterogeneous streambed sediments. Both looked at bedform-scale hyporheic zone development under hypothetical boundary conditions and the top boundary of their models represent the stream-bed interface. Here, we expand the application of their technique to the reach scale and include analysis of both lateral and vertical hyporheic zone development. We also calibrate the model to actual field conditions, including stream stage and subsurface water chemistry. In this paper, we report the results of using two numerical modeling packages, MODFLOW and MT3D, to quantify the degree of hyporheic interaction along an experimental reach of Red Canyon Creek, Wyoming. Our modeling experiments were done specifically to identify

the locations of maximum hyporheic interaction along our experimental reach in the creek’s watershed, and to examine the geomorphic and hydrologic features driving the exchange at those sites, such as debris dams and meanders. In this study, we quantified the degree of hyporheic interaction based on the simulated flux of water through the hyporheic zone and the direct extent to which surface water derived solutes are present in the subsurface. 2. Study site The study site is a reach of Red Canyon Creek, a second order meandering stream between 1750 and 2000 m (5750 and 6550 ft) above sea level along the southeastern flank of the Wind River Range of the Rocky Mountains (Fig. 1). Red Canyon was cut along the bedding plane between the eastward dipping Triassic age Chugwater Red Beds and the underlying Permian age Phosphoria Formation. The 80 km2 Red Canyon Creek watershed contains several ephemeral streams and the perennial Cherry Creek, which joins Red Canyon Creek near its confluence with the Little Popo Agie River at the base of the watershed.

Fig. 1. Map of the Red Canyon Creek study site, showing the locations of water table wells, piezometers, and in-stream mini-piezometers used for model calibration, and the locations of dams. The inset to the upper left shows the location of the site within the Red Canyon watershed and the inset map, bottom left, is an enlargement of the near-stream transect shown on the main map.

L.K. Lautz, D.I. Siegel / Advances in Water Resources 29 (2006) 1618–1633

The western side of the creek valley is bounded by a series of at least five sand and gravel terraces overlying the dip slope of the Phosphoria Formation. Approximately 35 cm of precipitation falls annually in this southwestern Wyoming semi-arid ecosystem region [19]. Average monthly stream discharge rates in the region can fluctuate by greater than a factor of 10 over the course of the year due to spring snowmelt. Peak flows occur during May and slowly decrease during baseflow recession. The study reach is located at the base of the watershed, downstream of where Cherry Creek joins Red Canyon Creek (Fig. 1). Here, Red Canyon Creek traverses a series of low relief meadows towards the creek’s confluence with the Little Popo Agie River. The study reach contains several small meanders and one severe meander, approaching cut-off. The meadow adjacent to the study reach is underlain by 6–7 m of unconsolidated alluvium. Silty sands extend from the surface to a depth of about 2 m where we find a discrete layer of sand and gravel. This highly permeable unit is underlain by silty sands and bedrock. The study site lies above the Chugwater Formation, which has a very low hydraulic conductivity relative to the unconsolidated alluvial deposits along the valley bottom. The overall slope of the channel across the study site is less than 2%. Along unobstructed reaches the slope is even smaller, about 0.3%, and interrupted by sudden drops in stream stage over beaver dams and small man-made dams. The stream is approximately 1.5–3.0 m wide and 0.15–0.35 m deep. Stream discharge at the study site during July 2004 was approximately 0.25 m3/s. 3. Methods 3.1. Field measurements A 320 m reach of Red Canyon Creek was selected for this study. During the summers of 2000 through 2004, the meadow adjacent to the study reach was instrumented with over 30 water table wells, piezometers and in-stream mini-piezometers as a part of hydrogeologic investigations at the University of Missouri-Columbia’s Branson Geology field camp (Fig. 1) [15]. Eight water table wells and 14 piezometers were distributed across the meadow to provide a good horizontal and vertical spatial distribution of hydrologic data. Water table wells and piezometers were installed using a Geoprobe (Kejr, Inc., Salina, Kansas), which is a hydraulically powered drill rig that uses both static force and percussion to advance sampling tools into the subsurface. Wells and piezometers were constructed of 2 cm diameter PVC casing to depths of between 2.1 and 6.1 m (7 and 20 ft), with the bottom screened at an interval of 0.6–3.0 m (2–10 ft), depending on well design. Water table wells were screened at the water table and piezometers were screened at intervals below the water table. Eleven in-stream mini-piezometers were installed along the study reach to obtain hydrologic data in a variety of streambed settings, including along pools and riffles, at

1621

Fig. 2. Histogram showing distribution of hydraulic conductivity measurements for the silty sand unit and the sand and gravel unit.

the immediate upstream and downstream ends of the severe meander, and immediately preceding and following flow obstructions. In-stream mini-piezometers were constructed of 2.5 cm outer diameter PVC pipe with slit screens along the bottom 10 cm. Mini-piezometers were installed by hand to a depth of between 25 and 40 cm in the streambed. Locations and casing elevations of all water table wells, piezometers, and mini-piezometers were surveyed during June 2004 using a Trimble total station. The total relief across the meadow is less than 2 m, so this high-resolution casing elevation data was necessary to accurately contour the meadow topography. Detailed longitudinal profiles of the streambed and water surface were created from measurements of the depth to the streambed and water surface from the surveyed mini-piezometer casings. Hydraulic head measurements were taken at all water table wells, piezometers, and mini-piezometers on June 23rd, 2004. Hydraulic conductivity values for the silty sand unit and the sand and gravel unit were measured at several wells and piezometers on-site using variable head tests [13] (Fig. 2). Water samples were collected from three monitoring wells, 13 piezometers and nine mini-piezometers on July 5th, 2004. Samples were pumped through silicone tubing to polypropylene sample bottles, via a peristaltic pump, and then refrigerated prior to analysis. Analyses for sodium and calcium were done using a Direct Current Plasma Emission Spectrometer (Beckman). 3.2. Groundwater flow modeling Groundwater flux at the field site was calculated using MODFLOW, a US Geological Survey block-centered finite-difference computer code that solves the groundwater flow equation [9,16]. We used the Visual MODFLOW preand post-processors (Waterloo Hydrogeologic, Inc.) and the model was run to steady-state. Our three-dimensional model domain was 100 m by 200 m, which we simulated using a 0.5 m uniform horizontal grid spacing sufficient to capture sub-meter scale hypor-

1622

L.K. Lautz, D.I. Siegel / Advances in Water Resources 29 (2006) 1618–1633

heic flow paths. The model has four layers, with variable vertical grid spacing, for a total model depth of between 6.3 m and 7.5 m. Model topography (the top of the first layer) was imported from elevation data for the well casings and a topographic map of the site, generated by hand from survey data. The elevation of the top of the second model layer coincides with the elevation of the streambed along the reach. This configuration allowed us to place mini-piezometer screens in the model immediately below the base of the stream and outside of cells assigned constant head boundary conditions. The elevations of the top and bottom of the third model layer coincide with the elevations of the top and bottom of the sand and gravel unit measured from soil borings at each well installation. As a first approximation, hydraulic conductivity was modeled as isotropic and heterogeneous (Table 1). The top two layers and the bottom layer of the model contain the silty sand unit and were assigned a hydraulic conductivity of 3 · 104 cm/s. The third layer contains the sand and gravel unit and was assigned a hydraulic conductivity of 103 cm/s. These hydraulic conductivity values are consistent with the geometric mean of measurements taken at the site via variable-head tests (Fig. 2). The groundwater recharge rate in this semi-arid climate is minimal during summer. In June 2004, precipitation, as rainfall, was 43 mm (1.7 in) at the study site. Assuming a modest 5% of precipitation was not lost to evapotranspiration at this time, a recharge rate of 2.2 mm/month was assigned to the top of the model (Table 1). Sensitivity analysis suggests that the model is not sensitive to the recharge rate. The eastern and western boundaries of the model domain are parallel to sharp breaks in valley topography at the base of stream terraces. These boundaries also are parallel to the general direction of groundwater flow inferred from water table maps [15]. Little, if any, groundwater flux occurs across these flowlines, so they were assigned the Neumann no-flow boundary condition [1]. The northern and southern boundaries of the model are perpendicular to the inferred direction of groundwater flow. These boundaries cross meadows similar to that found in the center of the study site. The flux rate, per unit area, across these boundaries was assumed to be the same as the flux rate, per unit area, across the central portion of the study site. We assumed a rate of 3 · 103 m3/day per square meter as a first approximation based on the modeled flux rate across the central portion of the model. To Table 1 Hydrologic parameters used in the MODFLOW/MT3D simulation Parameters

Value

Horizontal grid size K (silty sand unit) K (sand and gravel unit) Effective porosity (ne) Recharge (R) Channel gradient

0.5 m · 0.5 m 0.0003 cm/s 0.001 cm/s 0.30 2.2 mm/month 0.02

simulate the flux boundary conditions, a series of pumping wells were simulated to create a constant flux boundary [1]. The bottom of the model corresponds to bedrock and was assigned the no-flow boundary condition. We modeled the stream using a spatially variable constant head boundary condition in the uppermost layer of the model. The location and width of the constant head boundary were designated based on an aerial photograph of the site, coupled with elevation data from direct surveys that were imported into Visual MODFLOW. The average gradient of the stream was only 0.3% along unobstructed segments and just under 2% overall. The constant head values were assigned as linear gradients between surveyed stream surface elevations (at the locations of in-stream mini-piezometers) and in-stream flow obstructions. Sudden drops in stream stage over flow obstructions, including small beaver dams and a man-made dam, were simulated with sudden drops in constant head, corresponding to the actual changes in water elevation at the locations of these structures. 3.3. Solute transport modeling To simulate surface and ground water mixing in the hyporheic zone, we used MT3D, a solute transport package included with Visual MODFLOW Pro (Waterloo Hydrogeologic, Inc.). MODFLOW solves for the distribution of hydraulic head within the model domain and, from these results, the velocity components of flow are calculated. MT3D uses the velocity values as inputs to solve the transport equation. Macroscopic heterogeneities are the primary reason solutes do not travel along purely advective flow paths in field experiments [31]. As solutes are transported through the subsurface, small-scale heterogeneities in the porous media cause flow velocity variations, which cause longitudinal and transverse spreading of solutes. To practically simulate the effect of these small-scale heterogeneities, a dispersive transport calculation is used [31]. The MT3D solution reflects concentration changes that are caused by both dispersion and sink/source mixing [30,32]. In Visual MODFLOW, several solution methods are available to solve the solute transport equation and they can be categorized as particle-based methods, finite-difference methods and a total-variation-diminishing (TVD) method. Briefly, the particle-based methods, Method of Characteristics (MOC), Modified Method of Characteristics (MMOC) and Hybrid Method of Characteristics (HMOC), solve the advection term using conventional particle tracking based on a mixed Eulerian–Lagrangian approach and solve the dispersion and source/sink mixing terms using finite difference [30,32]. The MOC method is virtually free of numerical dispersion, so is most appropriate for problems with sharp concentration fronts, but requires a large amount of computing memory. MMOC is more computationally efficient than MOC, but can introduce some numerical dispersion in problems with sharp

L.K. Lautz, D.I. Siegel / Advances in Water Resources 29 (2006) 1618–1633

concentration fronts. HMOC combines the strengths of MMOC and MOC, providing more accurate solutions in problems with both sharp and non-sharp concentration fronts. Finite-difference methods, such as the Upstream FiniteDifference Method (UFD), and TVD are typically more computationally efficient than particle-tracking methods, but can result in significant numerical dispersion, particularly for problems that are advection-dominated and have sharp concentration fronts [26,30,32]. To maximize computational efficiency and prevent numerical dispersion in our solution, we used the Hybrid Method of Characteristics (HMOC) solver. To assess the impact of solver selection on model results, we ran the model using four other solvers included with Visual MODFLOW (MOC, MMOC, UFD and TVD). A more detailed discussion of the concepts and fundamental ideas behind each of these methods can be found in [31] and in the MT3D User’s Manual [30]. The solute transport package uses effective porosity to convert Darcy velocities to linear velocities and uses longitudinal, horizontal, and vertical dispersivities for the finitedifference solution to the dispersion term. We selected a porosity value of 0.30, consistent with literature values for moderate to poorly sorted sands and gravels [7]. Longitudinal dispersivity is scale-dependent and, by convention, can be considered approximately as one-tenth of the flow path length. Hyporheic flow paths are expected to be centimeters to tens of meters in length [11], so a longitudinal dispersivity value of 0.1 m was chosen for this study as a first approximation. Although smaller longitudinal dispersivities (0.01–1 cm) are expected in laboratory experiments, values of the order of magnitude we selected (0.1 m) are needed to account for dispersion at the field scale, which is caused by macroscopic heterogeneities, rather than pore-scale dispersion [31]. This value is consistent with the most reliable longitudinal dispersivity data for similar scales presented by Gelhar et al. [8]. Longitudinal dispersivity is typically greater than transverse and vertical dispersivity [7,22]. We assumed transverse dispersivity is one-tenth the longitudinal dispersivity and vertical dispersivity is one-hundredth the longitudinal dispersivity. There are two sources of water to the hyporheic zone— surface water and ground water. To identify the percent mixing of these two sources, we assigned an arbitrary concentration of 100 ppm to the surface water and 0 ppm to the ground water [3,27]. All constant head boundary cells, representing the creek, were designated as 100 ppm constant concentration boundaries. All other cells in the model were assigned an initial concentration of 0 ppm and other boundaries are convective flux boundaries. The transport model was run at a 0.5-day time interval for a total of 10 days. Near-stream flow paths along the experimental reach are long (with travel times on the order of days to years), so we could not logistically do in situ field subsurface tracer experiments to fully calibrate our MT3D modeling compo-

1623

nent. However, we used the ratio of Na:Ca in surface and subsurface water samples as a surrogate for the extent of hyporheic zone interaction. Along Red Canyon Creek, cation exchange in the subsurface replaces calcium in surface water with sodium, causing Na:Ca ratios to increase as water moves along near-stream flow paths. In areas with rapid influxes of surface water, calcium has not been replaced with sodium and water has relatively low Na:Ca ratios. In contrast, groundwater has higher Na:Ca ratios as a result of greater cation exchange. For this reason, we can use Na:Ca ratios to identify sites of surface and groundwater mixing in near-stream sediments along Red Canyon Creek. 3.4. Quantifying hyporheic exchange We delineated the hydrochemically active hyporheic zone within our model domain based on a new definition, specifically as near-stream subsurface zones receiving a minimum of 10% surface water [24] (analogous to a concentration of 10 ppm) within a 10-day travel time period [29]. Results of the MT3D simulation were used to delineate these regions. Model cells with a final concentration of 10 ppm or higher were considered part of the hyporheic zone [27]. The volume of the hyporheic zone was determined by summing the saturated volume of the hyporheic zone cells within the model. We quantified flux rates of water into and out of the hyporheic zone using cell-by-cell flux rates exported from MODFLOW. We measured the flux of water in and out of the constant head cells versus the distance along the reach. For comparison to the MT3D simulation, we used MODPATH to simulate subsurface flow paths originating in the stream. We introduced one particle into each constant head boundary cell, for a total of 3447 particles, and tracked these particles until they returned to the stream or left the model domain through the northern constant flux boundary. Although the model allows for the addition of more particles, one per cell was adequate to identify detailed distribution of flow paths. 3.5. Sensitivity analysis To assess the potential impact of numerical dispersion on the solute transport model results, we did both grid size and solver sensitivity analyses. We also assessed the model sensitivity to hydraulic conductivity (K), anisotropy (Kh:Kv), effective porosity (ne), dispersion (D) and recharge (R) because of the inherent difficulty in obtaining accurate estimates of these parameters. We were also interested in the impact of advective versus dispersive processes on the extent of the hyporheic zone. For the sensitivity analysis, values for K, ne, D and R were increased by 20% and 50%. K values for all layers were increased by the same proportion simultaneously. Anisotropy, the ratio of horizontal K to vertical K, was increased from 1 (isotropic) to 5 and then 10. We looked

1624

L.K. Lautz, D.I. Siegel / Advances in Water Resources 29 (2006) 1618–1633

at how the total hyporheic volume in the model changed under the calibrated model conditions and with incremental changes in the assumed model parameters. We also looked specifically at the size of the hyporheic zone at three transects across the stream—at a man-made dam, a debris dam and a severe meander—under the calibrated model conditions and with the same incremental changes. 4. Results and interpretation 4.1. Groundwater flow simulation The model accurately simulated the distribution of hydraulic head across the meadow (Fig. 3). For the 33 head calibration points, the maximum residual is 0.26 m and the absolute residual mean is 0.11 m. The normalized residual mean square error of the model is 6.47%. The observed versus predicted head values have a correlation coefficient of 0.993. In general, the model-simulated heads are slightly higher than the observed heads. Days before the field observations of stage and head were made, a precipitation event slightly elevated stream stage along the reach (by approximately 2 cm), and so the stream stage used for the model simulation is slightly higher than the stage of the channel at equilibrium with the groundwater. No change in head at the wells or piezometers was observed in relation to the event. As a result, head elevations simulated by the model are slightly too high, as they are at equilibrium with a slightly elevated stream stage. The magnitude and direction of vertical hydraulic gradients observed in the streambed calibrate well to model predicted gradients (Fig. 4). In situ field hydraulic gradients were calculated from stream water elevations and measures of hydraulic head at in-stream mini-piezometers, screened between 25 and 40 cm below the streambed. Simulated vertical gradients are calculated based on the difference between constant head values for the stream boundary and an interpolation of head between constant head nodes and underlying neighboring nodes. The calibration indi-

Fig. 3. Observed versus predicted hydraulic head values for all water table wells, piezometers, and in-stream mini-piezometers.

Fig. 4. Observed versus predicted vertical hydraulic gradients for all instream mini-piezometers.

cates the model accurately simulates hydraulic gradients over distances of only centimeters. If anything, the model underestimates streambed gradients, and therefore, our model conservatively simulates hyporheic exchange. As a first approximation, we assumed hydraulic conductivity was isotropic at the study site. This assumption provided the best calibration to observations of hydraulic head and streambed gradients. Increases in anisotropy (Kh:Kv = 5 and 10) produce a poorer fit to field data (Fig. 5), increasing the normalized root mean error of hydraulic head and stream gradient predictions.

Fig. 5. Sensitivity of model calibration to increasing anisotropy.

L.K. Lautz, D.I. Siegel / Advances in Water Resources 29 (2006) 1618–1633

Model results also calibrated well to total streambed flux rates across the site, measured from consecutive stream discharge measurements obtained by dilution gauging on June 30th, 2004. During the dilution experiments, plateau concentrations at consecutive points along the reach were not statistically different from one another (p = 0.15), indicating no change in the stream discharge rate along the reach (e.g. net groundwater inflow is zero). The model calculated a net stream inflow rate of 0.025 l/s (0.01% of the total discharge rate in the stream)—an inflow rate that is essentially zero, below the detection limits of dilution gaging. Water table contours generated from the groundwater flow simulation show that water generally flows across the experimental reach meadow, discharging from the stream at the southern end of the stream reach and reentering the stream at its northern end (Fig. 6). This flow pattern is consistent with previous observations of site hydrology [15]. Over the entire modeled reach, the creek loses 0.28 l/ s to the subsurface and gains 0.30 l/s from the subsurface. There is negligible net change in stream flow along the reach. Along the most upstream (southeastern) segment of the reach, the creek loses water to the predominantly northern groundwater flow path through the meadow. Immediately downstream, the flow pattern reverses, and the stream gains water from the upstream meadow (simulated as the flux boundary) and from flow temporarily diverted out of the creek by a dam. This pattern of flow reversal is repeated

1625

along the most downstream (northern) segment of the reach: initially the creek loses water to the downstream meadow and then the flow pattern changes, gaining water from the upstream meadow. Although the stream water enters the subsurface and returns to the stream along these flow paths, the distance traveled (approximately 150 m) and the long residence time in the subsurface (several years) precludes considering these pathways as part of the hyporheic zone. Smaller flow systems are embedded within this overall groundwater flow setting. The stream repeatedly alternates between gaining water and losing water to the subsurface, as indicated by the alternating ‘‘v’’ pattern of equipotential lines immediately adjacent to the stream (Fig. 6). Also, flow budgets for individual constant head cells show flux rates alternate between negative, where water is entering the subsurface, and positive, where water is flowing into the channel (Fig. 7). Flux rates for individual constant head cells range from about 0.04 m3/day to about 0.04 m3/day. The individual cell flux rates confirm that stream inflows and outflows are balanced along the reach and that the stream repeatedly alternates between gaining and losing. The reversals of stream bed flux are strongly associated with the presence of dams within the creek. These conditions are verified by our field observations of alternating positive and negative vertical hydraulic gradients in the streambed (Fig. 4). The pattern of equipotential lines and gradients indicate water is exiting the channel along both lateral and vertical flow paths.

Fig. 6. Model domain showing hyporheic zones along Red Canyon Creek. The constant head boundary is shown in black. Particle pathlines shown here are a select summary of actual pathlines generated by MODPATH.

1626

L.K. Lautz, D.I. Siegel / Advances in Water Resources 29 (2006) 1618–1633

Fig. 8. Na:Ca ratios measured in samples taken along a transect from the stream to P14. Fig. 7. Constant head cell flux rates plotted versus distance downstream. Negative flux rates indicate the stream is losing water to the subsurface. Positive flux rates indicate the stream is gaining water from the subsurface.

Particle tracking shows that the majority of particles simulated in the model exited and returned to the channel within the model domain (Fig. 6). Particle residence times ranged from several hours to 10 years. Using a purely hydrological definition of the hyporheic zone, we would conclude that the majority of water transported through the model domain is hyporheic water. This is in sharp contrast to our hydrochemical definition, which limits the hyporheic zone to regions immediately adjacent to the creek at points of high flux from the stream into the subsurface. 4.2. Solute transport simulation The solute transport simulation shows high concentrations of surface water (>10%) in the subsurface only at nine discrete locations along the creek. These locations correspond to places where flux rates show the creek is losing water to the subsurface near dams and highly sinuous segments of the creek. We used the model results to delineate the extent of the hyporheic zone along the creek (Fig. 6). To assess the validity of our solute transport simulation, we compared field measurements of Na:Ca ratios at sites located within and outside of the model hyporheic zone to assess the influence of stream water on subsurface water chemistry at these sites. We hypothesize that calcium-rich stream water (2.1 mmol/l), which comes from the dissolution of gypsum and calcium carbonate within the watershed, undergoes ion exchange as water moves through the subsurface. The ratio of Na:Ca, measured in water samples taken along a transect of wells extending from the stream to P14, significantly increases a short distance away from the stream (Fig. 8). The change in the Na:Ca ratio indicates the concentration of sodium increases more than the concentration of calcium, suggesting ion exchange is replacing dissolved calcium with sodium. Six of our in-stream mini-piezometers (MP-A through MP-F) fall within the hyporheic zone, as delineated within the model domain using the hydrochemical definition. The percent surface water at these locations, as predicted by the

model, ranges from 10.7% to 95.7%. All other mini-piezometers, monitoring wells and piezometers fall outside of the hyporheic zone (<10% surface water). Comparisons of Na:Ca ratios measured in samples taken from the stream, piezometers and wells show how rapid influxes of stream water influence the subsurface water chemistry at MP-A through MP-F. The stream water has a low Na:Ca ratio of 0.10, compared to the mean value in ground water, 0.18. The mean Na:Ca ratio at sites MP-A through MP-F, 0.14, is significantly different from the mean ground water ratio (p = 0.01) and falls between the stream and ground water values. The distribution of the ground water and hyporheic zone Na:Ca ratios are shown in Fig. 9. The lower Na:Ca ratios in water samples collected within the simulated hyporheic zone suggest that the model is accurately predicting the rapid influx of surface water within the hyporheic zone at these sites. The model results show nine hyporheic zones with >10% surface water within the 10-day model simulation (Table 2). At Zones 1 and 7, near the southern and northern flux boundaries of the model, respectively, the creek loses water to the subsurface, where groundwater then flows north into the respective downstream meadow. At Zone 8, a small beaver dam (drop in stream stage of 0.08 m) creates a small

Fig. 9. Histogram showing the distribution of Na:Ca ratios measured in water samples taken from sites located outside the simulated hyporheic zone (groundwater), inside the simulated hyporheic zone and the stream.

L.K. Lautz, D.I. Siegel / Advances in Water Resources 29 (2006) 1618–1633

1627

Table 2 Characteristics of the hyporheic zones identified along Red Canyon Creek Hyporheic zone

Volume (m3)

Average cross-sectional area (m2)

Stream stage drop (m)

Sinuosity

Hydrologic setting

1 2 3 4 5 6 7/8 9 Overall

35.90 61.92 31.17 60.26 28.27 30.53 73.69 37.44 487.37

0.8 2.0 1.2 2.0 1.0 2.0 2.4 1.3 1.2

N/A 0.34 0.15 0.27 N/A 0.13 0.08 N/A

1.3 1.4 1.4 1.2 1.7 1.2 1.1 2.1

Losing Gaining Balanced Balanced Losing Losing Losing Gaining

vertical gradient, that drives water into the subsurface. Zones 7 and 8 overlap, so it is not possible to identify the boundary between the two and we have combined them for this analysis. Zones 2, 3, 4, and 6 are located directly behind dams. The drop in stream stage over these dams ranges from 0.13 m to 0.34 m (Table 2). Zones 5 and 9 are located at points of high sinuosity (1.7 and 2.1, respectively) along the creek, with Zone 9 extending through the nearly cutoff meander bend. Simulated flux rates for individual constant head cells confirm that surface water is moving rapidly into the subsurface at these locations (Fig. 7). For the entire reach and each individual zone, we tabulated the total volume of the hyporheic zone (Table 2). The solute transport model has a very sharp concentration front and is therefore subject to numerical dispersion, which can affect the model results. To evaluate the sensitivity of model-derived hyporheic volumes to numerical dispersion, we did a grid size sensitivity analysis and a solver sensitivity analysis. The original model contains four layers, with the extent of the simulated hyporheic zone extending only into layers 1 through 3. Layer 1 contains the constant head boundary and head values are typically only slightly above the bottom of the layer. For this reason, we did not subdivide the top layer

of the model for our grid size sensitivity analysis. We refined layers 2 and 3 by first subdividing each layer in half, for a total of 6 layers, and then divided those layers in half again, for a total of 10 layers. For the 6 layer and 10 layer models, the total volume of the hyporheic zone changes by only 2.2% and 4.0%, respectively (Fig. 10). The model is relatively insensitive to changes in grid size and it is unlikely that our results are affected by numerical dispersion caused by grid resolution. The original model was solved using HMOC, which solves the advection term using conventional particle tracking based on a mixed Eulerian–Lagrangian approach and solves the dispersion and source/sink mixing terms using finite difference (see Section 3). Lagrangian methods are particularly suitable for solving advection-dominated problems because they are virtually free of numerical dispersion and Eulerian methods are better suited for dealing with dispersion [26,31]. HMOC combines the strengths of both of these methods. We reran the model using four other solvers, MOC, MMOC, UFD and TVD, to assess the sensitivity of model-derived hyporheic volume to solver selection. Using other particle-based methods (MOC and MMOC), the total hyporheic zone volume changes by only 0.3% and 1.2%, respectively (Fig. 10). MOC is the most computationally expensive solver, but has the advantage of

1628

L.K. Lautz, D.I. Siegel / Advances in Water Resources 29 (2006) 1618–1633

Fig. 10. Sensitivity of total hyporheic volume to increasing grid refinement and solver selection. MOC = Method of Characteristics, MMOC = Modified Method of Characteristics, UFD = Upstream Finite Difference, TVD = Total Variation Diminishing Scheme.

Fig. 11. Sensitivity of total hyporheic volume to changes in dispersion (D), hydraulic conductivity (K), anisotropy (Kh:Kv), effective porosity (ne) and recharge (R). Changes in parameter values are indicated on the plot, in the bars.

being virtually free of numerical dispersion. The HMOC results and MOC results are very similar, suggesting that the use of HMOC improved computational efficiency without introducing significant numerical dispersion. MMOC is also more computationally efficient than MOC, but can introduce some numerical dispersion, which may explain the 1.2% change in total hyporheic zone volume.

UFD and TVD improve computational efficiency of the solution, but can introduce significant artificial oscillation and numerical dispersion, particularly when applied to advection-dominated transport problems, such as in this study. Our sensitivity analysis shows that selection of UFD and TVD for the solver caused significant numerical dispersion; we observed changes in total hyporheic volume

Fig. 12. Plan views of the water table and cross-sectional views of equipotential lines under isotropic and anisotropic (Kh:Kv = 10) conditions. The crosssectional view is taken along the western side of the model, down the length of the stream.

L.K. Lautz, D.I. Siegel / Advances in Water Resources 29 (2006) 1618–1633

of 29.6% and 64.2%, respectively, when these solvers were used (Fig. 10). The mixed Eulerian–Lagrangian methods (MOC, MMOC and HMOC) introduce the least numerical dispersion and are most appropriate for problems with sharp concentration fronts, such as ours. These methods combine the strength of particle tracking, which is extremely effective for solving the advection term, and the finite-difference method, which is more numerically stable for solving the dispersion term [31]. Sensitivity analysis of the total hyporheic volume indicates that changes in anisotropy, hydraulic conductivity and effective porosity cause the largest changes in the hyporheic volume (Fig. 11). Notably, increasing the anisotropy (Kh:Kv) to 10 caused over a 50% decrease in the volume of the hyporheic zone. The effect of anisotropy on the hydrologic simulation is most apparent in areas close to the stream, as illustrated by plan and cross-sectional views of

1629

model results under isotropic conditions (Kh:Kv = 1) and anisotropic conditions (Kh:Kv = 10) (Fig. 12). In plan view, at points away from the stream, the water table contours are very similar, but the v’ing of contours around the channel is more pronounced under anisotropic conditions. In cross-section the equipotential lines appear similar at depth, but anisotropy leads to steeper gradients near the stream boundaries. The sensitivity analysis suggests that it is important to calibrate groundwater flow models to hydraulic head and hydraulic gradients in the streambed; inaccurate assumptions of anisotropic conditions can greatly limit simulated hyporheic zone volume. Changes in assumed dispersivity up to 50% caused no change in the hyporheic zone volume. Cross-sections for transects located at the man-made dam, a beaver dam, and the severe meander show the extent of the hyporheic zone (Fig. 13). Sensitivity analyses

Fig. 13. Model sensitivity to water and solute transport parameters, including hydraulic conductivity (K), effective porosity (ne), anisotropy (Kh:Kv) and dispersion (D), at three different transects. A to A 0 is a transect across the channel 2 m upstream of the man-made dam; B to B 0 is a transect 2 m upstream of a beaver dam; C to C 0 is a transect across a very tight meander bend (the transect crosses the channel at the upstream and downstream ends of the bend). The shaded gray areas show the location of the highly permeable sand and gravel unit. The dark black areas show the locations of constant head boundary conditions.

1630

L.K. Lautz, D.I. Siegel / Advances in Water Resources 29 (2006) 1618–1633

at these locations confirm that changes in modeled hydraulic conductivity and anisotropy generate the largest changes in the extent of the hyporheic zone. This result shows that the movement of surface water into the hyporheic zone is predominantly an advective process. The hyporheic zone located at the man-made dam (Transect A–A 0 ; Zone 2) shows the greatest sensitivity to changes in model parameters. This is because the transect at this location is

Fig. 14. Plot of hyporheic zone volumes versus cross-sectional areas. Each point is labeled with the corresponding zone number.

parallel to the predominant groundwater flow direction across the meadow, so there is a greater horizontal hydraulic gradient across the transect here. In contrast, the other two transects are perpendicular to the predominant groundwater flow direction. Hyporheic zones underlying longer segments of the creek will naturally be larger than those underlying shorter segments. For this reason, absolute hyporheic volumes need to be considered relative to the length of the stream contributing water to that zone. To obtain an average cross-sectional area for each zone, we have normalized the volumes of the hyporheic zone to the corresponding length of stream. The average cross-sectional area of the hyporheic zones ranges from 0.8 m2 to 2.4 m2, from about equal to over twice the average cross-sectional area of the stream (about 1.0 m2). For the reach as a whole, the average cross-sectional area of the hyporheic zone is 1.2 m2. In general, the hyporheic zones were confined to the upper silty sand unit in the model. No surface water concentrations reached 10% in the deepest silty sand unit at the base of the model, and only hyporheic zone 2 extended into the sand and gravel unit to some extent. Zones associated with dams generally have larger cross-sectional areas (Fig. 14). This indicates that hyporheic zones created by dams extend further into the subsurface regardless of their absolute size. To illustrate this point we have included a three-dimensional isosurface plot, which outlines the maximum extent of the hyporheic zone along the creek

Fig. 15. Three-dimensional isosurface map showing the results of the MT3D solute transport simulation. The shaded outline shows the maximum extent of the hyporheic zone along the reach. The hyporheic zone is defined here as subsurface areas with greater than 10% surface water at the end of a 10-day model simulation. Subset A shows the shape of the hyporheic zone adjacent to a beaver dam. Subset B shows the shape of the hyporheic zone across a severe meander.

L.K. Lautz, D.I. Siegel / Advances in Water Resources 29 (2006) 1618–1633

(Fig. 15). In subset A, the debris dam is creating a deep hyporheic zone, whereas in subset B, the gradient through the meander bend creates a shallow hyporheic zone. 5. Discussion By defining the hydrochemically active hyporheic zone with MT3D, we have identified areas where we can expect to find high concentrations of surface–water solutes within a short period of time. Field measurements of Na:Ca ratios in the subsurface confirm that water chemistry within the simulated hyporheic zone is significantly different from the groundwater due to the rapid influx of surface water. The influx of solutes from surface water at these locations may balance the consumption of those solutes by microbial and geochemical processes. These areas are likely to have high concentrations of dissolved oxygen and may be sites of enhanced nutrient cycling and geochemical transformations of solutes. Although previous studies have used MODFLOW on studies of hyporheic exchange [3,14,23,27,29], our study extends MODFLOW applications to include simulating mixing and dispersion of surface and groundwater solutes both vertically and laterally through the hyporheic zone at the reach scale. This modeling strategy predicts the concentration of surface water solutes in the hyporheic zone, making it possible to ultimately simulate the loss of solutes, through adsorption and biogeochemical processes, as they mix with groundwater in the near-stream environment. Models of only advection-defined flow paths cannot be extended to include potential biochemical processes, such as denitrification, which are the focus of much of the hyporheic research [6,12,17,24]. Our model shows that movement of solutes into the hyporheic zone is predominantly driven by advection. The development of the hyporheic zone in our model is most sensitive to hydraulic conductivity and anisotropy. For this reason, when simulating hyporheic exchange using MODFLOW, it is most important to calibrate the groundwater flow model to hydraulic gradients, particularly in the streambed where the exchange takes place. We have been able to achieve that calibration through measurements made at in-stream mini-piezometers. One can also calibrate the MT3D solute transport simulation to the concentrations of surface water solutes in the near-stream sediments. Here we used changes in Na:Ca ratios, caused by ion exchange, to verify that the model accurately identified zones with rapid influxes of surface water. MT3D provides a more realistic simulation of hyporheic exchange processes than particle tracking and the HMOC solver is a computationally efficient solution of the transport equation that does not introduce numerical dispersion. Defining the hyporheic zone based solely on hydrologic flow path, versus a definition based on concentrations of surface water within the subsurface, provides a very different assessment of hyporheic exchange along

1631

Red Canyon Creek. If we used a purely hydrological definition of the hyporheic zone using MODPATH, we would have included very long near-stream flow paths with residence times of up to several years within our hyporheic zone. Chemical transformations occurring over time along these flow paths, including consumption of oxygen, nutrient transformation and sorption of dissolved solids, can significantly change the chemistry of the water as it moves through the subsurface. Ion exchange in the subsurface along Red Canyon Creek is replacing calcium with sodium as surface water moves along subsurface flow paths. This causes significant changes in Na:Ca ratios as water moves from the stream, through the hyporheic zone, becoming ground water. By limiting the travel time, our definition of the hyporheic zone is limited to regions immediately adjacent to the stream channel. Here, the chemistry of the subsurface waters corresponds more with stream water chemistry. By incorporating both residence time and mixing of surface and ground waters in our simulation, we have more accurately simulated the distribution of hyporheic waters along Red Canyon Creek. The hyporheic zone at Red Canyon Creek particularly is not uniformly distributed—dams and meanders create localized areas of high surface water concentrations in the subsurface (Fig. 6). Transport of high concentrations of surface water into the subsurface in a short time has implications for nutrient processing and geochemical transformations in the hyporheic zone. Surface water is oxygenrich and its presence in the subsurface could create tightly spaced oxic and anoxic zones. Large changes in stream stage over dams have created large hydraulic gradients and streambed flux rates, driving surface water into the subsurface in large and deep hyporheic zones. This result suggests that placing small flow obstructions, such as beaver dams or log dams less than 0.5 m high, along gaining reaches is a viable strategy to enhance hyporheic exchanges. Hyporheic zones generated by dams are unique because much of the surface–groundwater exchange takes place along deep, vertical flow paths, rather than horizontal, shallow flow paths. Stream segments with tight meanders (high sinuosity) also create distinct hyporheic zones but these zones have smaller cross-sectional areas and do not extend far into the subsurface. Along gaining or static reaches, the hyporheic zone only occurs near dams or meanders and the zone size is constrained by the groundwater flow pattern towards the stream. In contrast, along losing reaches of the creek, hyporheic zones are almost ubiquitously present and larger than along gaining reaches. Along losing reaches, surface water enters the subsurface regardless of the in-stream features—dams and meanders simply enhanced hyporheic zone development. Although more surface water is present in the subsurface along losing reaches, the water is not quickly diverted back to the stream, but captured by longer near-stream groundwater

1632

L.K. Lautz, D.I. Siegel / Advances in Water Resources 29 (2006) 1618–1633

flow paths. This is particularly true for zones created at the upstream and downstream ends of meadows. A limitation of defining hyporheic zones based on surface water content is that the definition does not consider the ultimate fate of the hyporheic water. Although the water in these hyporheic zones does not rapidly return to the channel, these sites are still where surface water rapidly moves into the subsurface and may generate nearby zones of enhanced nutrient and geochemical transformations. Biogeochemical processing in the hyporheic zone occurs where there are rapid influxes of oxygen and surface water solutes, regardless of the whether those hyporheic waters rapidly return to the channel. The results of this model are consistent with results of an in-stream tracer test done to measure hyporheic exchange along this reach in June 2003 [15]. At that time, flow conditions were very similar to those observed in July 2004 for this study. From the in-stream tracer test, we concluded that dams positioned along gaining segments of the reach increased hyporheic interactions while those positioned along losing segments diverted water to the subsurface where it was lost to long time scale flow paths. In that study, we suggested that while dams enhance hyporheic interaction, their position is an important consideration for maximizing their effectiveness. Here we have shown that positioning dams along gaining or non-losing (balanced inflow and outflow) segments leads to larger, deeper hyporheic zones where flow is likely diverted back to the main channel. 6. Conclusions Groundwater flow and solute transport modeling, using MODFLOW and MT3D, is effective for simulating the development of the hydrochemically active hyporheic zone. The simulation, which solves for both advection, dispersion and mixing of solutes, allows us to identify percentages of surface water present in the subsurface within a set time period, in this case 10 days. Defining the hyporheic zone as subsurface areas comprised of 10% surface water within 10 days, we can use the model results to calculate a threedimensional volume of the hyporheic zone. We have used this type of analysis to identify how different stream features enhance the development of the hyporheic zone, from a physical process based perspective. Pairing MODFLOW with MT3D to simulate the hydrochemically active hyporheic zone allows for calibration to both hydrologic and chemical field observations. Because the development of the hyporheic zone is most sensitive to hydraulic conductivity and anisotropy, it is most important to collect field measurements of those terms. Contrasts between stream water and groundwater chemistry can be used to assess how well MT3D simulates surface and groundwater mixing in near-stream areas. For this reason, it is important to collect representative samples of water chemistry from groundwater wells, the stream and in the

near-stream zone, such as with mini-piezometers in the streambed. This modeling strategy is more effective than particle tracking for simulating surface–ground water interactions because model predictions can be easily compared to field observations at near-stream wells. In the field, it would be difficult to track the movements of an individual parcel of stream water through the subsurface, but relatively simple to monitor the mixing of groundwater with tracer-laden stream water (i.e. [25]). The modeling strategy used here presents an ideal tool that can be used to predict how surface and ground water will mix over time along subsurface flow paths. The results of this modeling can be used to guide monitoring well placement during tracer tests and be compared to field observations of tracer-laden stream water moving into the subsurface. Debris dams are a key driver of rapid surface water– groundwater exchange and their installation may be a viable strategy for influencing nutrient transformations along Red Canyon Creek. Debris dams located along gaining segments of the creek will create deeper hyporheic zones where water is eventually diverted back to the main channel. The flux of oxygen-rich surface waters into these deeper zones will likely create tightly spaced oxic and anoxic environments where nutrient transformations can take place. Debris dams also enhance hyporheic zone development along losing reaches, but the water diverted to the subsurface is likely lost to longer time-scale flow paths. Highly sinuous reaches also enhance hyporheic zone development, but these zones are generally shallower and smaller. We have used a 10-day travel time window to represent a reasonable time for surface water solutes, such as oxygen or nutrients, to remain in hyporheic water before being subject to chemical or biological reactions that remove or transform solutes. This reaction time is highly variable across different stream settings, so in the future, we suggest developing MODFLOW and MT3D models that incorporate decay functions for solutes as they move through the subsurface. For example, a decay rate representing oxygen depletion in the subsurface could be used to simulate the concentrations of oxygen in the hyporheic zone at steady state—where rates of oxygen depletion are balanced with oxygen replenishment from surface water. Currently, MT3D allows for the simulation of first-order decay reactions and sorption of solutes. RT3D, based on MT3D, allows for the simulation of more complex chemical reactions, including three-dimensional, multi-species reactive transport. We anticipate that MODFLOW and MT3D simulations can be paired with long-term tracer tests to inversely model the decay rates of nutrients in the hyporheic zone. Hydrologic parameters in the model, including hydraulic conductivity and effective porosity, can be calibrated to the concentrations of a conservative tracer in the hyporheic zone. Then solute transport parameters, including decay

L.K. Lautz, D.I. Siegel / Advances in Water Resources 29 (2006) 1618–1633

rates, can be calibrated to the concentrations of nutrient (or other solute) tracers in the hyporheic zone. This type of analysis is not possible using MODPATH simulations and will be the subject of future work. Acknowledgments The National Science Foundation funded this work. We would like to thank The Wyoming Nature Conservancy for allowing us experimental use of this research site. We would also like to thank the University of Missouri’s Branson Geology Field Camp for help with site instrumentation. References [1] Anderson MP, Woessner WW. Applied groundwater modeling: Simulation of flow and advective transport. Academic Press; 1992. [2] Bencala KE, Walters RA. Simulation of solute transport in a mountain pool-and-riffle stream with a kinetic mass transfer model for sorption. Water Resour Res 1983;19(3):732–8. [3] Cardenas MB, Wilson JL, Zlotnik VA. Impact of heterogeneity, bed forms, and stream curvature on subchannel hyporheic exchange. Water Resour Res 2004;40:W08307. doi:10.1029/2004WR003008. [4] D’Angelo DJ, Webster JR, Gregory SV, Meyer JL. Transient storage in Appalachian and Cascade mountain streams as related to hydraulic characteristics. J North Amer Benthol Soc 1993;12(3):223–35. [5] Findlay S. Importance of surface–subsurface exchange in stream ecosystems: the hyporheic zone. Limnol Oceanograp 1995;40(1):159–64. [6] Findlay S, Strayer D, Goumbala C, Gould K. Metabolism of streamwater dissolved organic carbon in the shallow hyporheic zone. Limnol Oceanograp 1993;38(7):1493–9. [7] Freeze RA, Cherry JA. Groundwater. Englewood Cliffs, NJ: Prentice-Hall; 1979. [8] Gelhar LW, Welty C, Rehfeldt KR. A critical review of data on fieldscale dispersion in aquifers. Water Resour Res 1992;28(7):1955–74. [9] Harbaugh AW, Banta ER, Hill MC, McDonald MG. MODFLOW2000, the US Geological Survey modular ground-water models; user guide to modularization concepts and the ground-water flow process. US Geological Survey; 2000. [10] Harvey JW, Bencala KE. The effect of streambed topography on surface–subsurface water exchange in mountain catchments. Water Resour Res 1993;29(1):89–98. [11] Harvey JW, Wagner BJ. Quantifying hydrologic interactions between streams and their subsurface hyporheic zones. In: Jones JB, Mulholland PJ, editors. Streams and ground waters. San Diego, CA: Academic Press; 2000. [12] Hill AR, Labadia CF, Sanmugadas K. Hyporheic zone hydrology and nitrogen dynamics in relation to the streambed topography of a N-rich stream. Biogeochemistry 1998;42:285–310. [13] Hvorslev MJ. Time lag and soil permeability in ground water observations. US Army Corps of Engineers Waterway Experimentation Station, Bulletin 36, 1951.

1633

[14] Kasahara T, Wondzell SM. Geomorphic controls on hyporheic exchange flow in mountain streams. Water Resour Res 2003;39(1):1005. doi:10.1029/2002WR001386. [15] Lautz LK, Siegel DI, Bauer RL. Impact of debris dams on hyporheic interaction along a semi-arid stream. Hydrol Process 20(1):183–96. [16] McDonald MG, Harbaugh AW. A modular three-dimensional finitedifference ground-water flow model. Reston, VA: US Geological Survey; 1988. [17] Morrice JA, Dahm CN, Valett HM, Unnikrishna PV, Campana ME. Terminal electron accepting processes in the alluvial sediments of a headwater stream. J North Amer Benthol Soc 2000;19(4):593–608. [18] Morrice JA, Vallett HM, Dahm CN, Campana ME. Alluvial characteristics, groundwater–surface water exchange and hydrological retention in headwater streams. Hydrol Process 1997;11(3):253–67. [19] National Oceanic and Atmospheric Administration NCDC. Normal Monthly Precipitation (Inches), 2004. [20] Packman AI, Bencala KE. Modeling surface–subsurface hydrologic interactions. In: Jones JB, Mulholland PJ, editors. Streams and ground waters. San Diego: Academic Press; 2000. [21] Pollock DW. User’s guide for MODPATH/MODPATH-PLOT, version 3, A particle tracking post-processing package for MODFLOW, the US Geological Survey finite-difference ground-water flow model. US Geological Survey; 1994. [22] Schwartz FW, Zhang H. Fundamentals of ground water. New York: Wiley; 2003. [23] Storey RG, Howard KWF, Dudley WD. Factors controlling rifflescale hyporheic exchange flows and their seasonal changes in a gaining stream; a three-dimensional groundwater flow model. Water Resour Res 2003;39(2):1034. doi:10.1029/2002WR001367. [24] Triska FJ, Duff JH, Avanzino RJ. The role of water exchange between a stream channel and its hyporheic zone in nitrogen cycling at the terrestrial–aquatic interface. Hydrobiologia 1993;251:167–84. [25] Triska FJ, Kennedy VC, Avanzino RJ, Zellweger GW, Bencala KE. Retention and transport of nutrients in a third-order stream in Northwestern California: hyporheic processes. Ecology 1989;70(6):1893–905. [26] Waterloo Hydrogeologic I. Visual MODFLOW v. 4.1. Professional Edition, User’s Manual, Ontario, Canada, 2005. [27] Woessner WW. Stream and fluvial plain ground water interactions: rescaling hydrogeologic thought. Ground Water 2000;38(3):423–9. [28] Workshop SS. Concepts and methods for assessing solute dynamics in stream ecosystems. J North Amer Benthol Soc 1990;9(2):95–119. [29] Wroblicky GJ, Campana ME, Valett HM, Dahm CN. Seasonal variation in surface–subsurface water exchange and lateral hyporheic area of two stream-aquifer systems. Water Resour Res 1998;34(3):317–28. [30] Zheng C. MT3D: A modular three-dimensional transport model for simulation of advection, dispersion and chemical reactions of contaminants in groundwater systems, 1990. [31] Zheng C, Bennett GD. Applied contaminant transport modeling. New York, NY: Van Nostrand Reinhold; 1995. [32] Zheng C, Wang PP. MT3DMS: A modular three-dimensional multispecies model for simulation of advection, dispersion and chemical reactions of contaminants in groundwater systems: documentation and user’s guide. Vicksburg, MS: SERDP-99-1, US Army Engineer Research and Development Center; 1999.