Abstract

[1] In-channel stream restoration structures readjust surface water hydraulics, streambed pressure, and subsurface hyporheic exchange characteristics. In this study, we conducted flume experiments (pool-riffle amplitude of 0.03 m and wavelengths of 0.5 m) and computational fluid dynamic (CFD) simulations to quantify how restoration structures impacted hyporheic penetration depth, Dp, and hyporheic vertical discharge rate, Qv. Restoration structures were channel-spanning vanes with subsurface footers placed in the gravel bed at each riffle. Hyporheic vertical discharge rate was estimated by analyzing solute concentration decay data, and maximum hyporheic penetration depth was measured as the interface between hyporheic water and groundwater using dye tracing experiments. The CFD was verified with literature-based flume hydraulic data and with Dp and Qz observations, and the CFD was then used to document how Dp and Qz varied with flume discharge, Q, ranging from 1 to 15 L/s (3E + 03 < Re < 5E + 04). Flume experiments and CFD simulations showed that restoration structures increased Qz and decreased Dp, creating a shallower but higher flux hyporheic zone. Qz had a positive linear relationship with Q, while Dp initially grew as Q increased, but then shrunk when a hydraulic jump with low streambed pressured formed downstream of the structure. The restoration structures created counter-acting forces of increased downwelling head due to backwater effects, and increased upwelling due to low streambed pressure and standing waves downstream of the structure.

1. Introduction

[2] River steering structures, such as full-spanning and half-spanning vanes, are used in river restoration projects to redirect the primary flow maximum velocity into the channel center and to protect the meander cut bank from erosion [Radspinner et al., 2010; Rosgen, 2001]. In addition to river restoration structures regulating surface water hydrodynamics such as super-elevation, shear, and vortices [Zhou and Endreny, 2012], the river restoration structure also regulates the mixing of surface water and groundwater in a process called hyporheic exchange [Kasahara and Hill, 2006]. The restoration structures change water surface profiles, near bed velocities, and substrate permeability, and these factors contribute to hyporheic exchange and create a unique mixing of light, temperature, oxygen, carbon, and nutrients in the subsurface hyporheic zone [Harvey and Wagner, 2000]. Because of these distinct physical and chemical gradients, the hyporheic zone is a valued ecological niche providing habitat for benthic and interstitial organisms, spawning grounds for certain species of fish, a rooting zone for aquatic plants, as well as a set of nutrient transformation processes regulating water quality [Stanford and Ward, 1988; Gibert et al., 1994].

[3] Flume studies over the past two decades have provided a glass-walled view into the subsurface to document how channel bedform regulates hyporheic exchange [Elliott and Brooks, 1997a; Marion et al., 2002; Packman et al., 2004; Thibodeaux and Boyle, 1987]. The flume studies of Elliott and Brooks [1997a] and Thibodeaux and Boyle [1987] mapped how bedform undulations induced changes in bed pressure and generated distinct regions and paths of downwelling and upwelling hyporheic flux. Elliott and Brooks [1997b] also formulated a pumping model to predict the pressure and flux patterns along the flume bedform and mass flux across the entire flume bed. O'Connor and Harvey [2008] revisited the Elliott and Brooks [1997b] and solute mass-flux analysis from flume experiments by Marion et al. [2002] and Packman et al. [2004] to derive an effective diffusion equation to predict solute exchange across the bed based on bedform geometry and surface hydraulics. In a similar effort, Mutz et al. [2007] created a vertical discharge model to quantify the total volumetric exchange rate between the surface and subsurface. The effective diffusion and vertical discharge models are typically tested against observations of solute concentration decay as solute-free subsurface water mixed with solute-enriched surface water. The flume experiments provide local-scale information, illustrated by flow path direction and depth at a specific point, using the flume's glass-walled view. The flume experiments also provide reach-scale information, illustrated by mass flux or effective diffusion across the flume bed and its multiple wavelengths, using bulk concentration data. We use the local-scale and reach-scale terms to differentiate the hyporheic impacts of an individual structure and/or bedform from the integrated hyporheic impacts of multiple bedforms and/or structures.

[4] River engineers and scientists have called for investigations into the influence of river restoration on hyporheic exchange at the local scale and reach scale [Boulton, 2007; Boulton et al., 2010; Hancock, 2002; Hester and Gooseff, 2010; Krause et al., 2011]. While field experiments do not have glass-walled viewing of the local exchange, Lautz and Fanelli [2008] used piezometer sampling to investigate the spatial biogeochemical patterns around a channel-spanning restoration structure in a pool-riffle river. The analysis indicated that the structures strongly interrupted the spatial patterns of various biogeochemical processes. Similar local-scale analysis of vertical hydraulic gradient and temperature time series showed how in-channel structures organize hyporheic flow, conducted by Kasahara and Hill [2006] on a log-weir and Crispell and Endreny [2009] on seven restoration structures built from boulders. The structures create a break in the bed topography that generates an upstream water surface profile convexity with hyporheic downwelling and a downstream water surface concavity with hyporheic upwelling, a pattern documented in natural step-pool systems by Harvey and Bencala [1993]. Hester and Doyle [2008] replicated this classic downwelling into upwelling pattern around half-spanning and full-spanning weirs and drop structures, showing how hyporheic exchange rates and paths intensified with head loss at the structure. At the reach scale, field experiments of hyporheic exchange have used solute injection and recovery with transient storage modeling (e.g., OTIS) to quantify average mass exchange for the entire streambed in natural [Bencala and Walters, 1983] and restored channels [Bukaveckas, 2007].

[5] Detailed flume experiments of restoration structure impacts on hyporheic exchange have provided insights to local-scale impacts and reach-scale impacts, but the two scales have not been combined into a single analysis. Such an analysis has the potential to provide mechanistic local-scale explanation for the reach-scale trends. In a local-scale study using a flume and a computational fluid dynamic (CFD) model, Endreny et al. [2011] mapped how an in-channel structure created complex pressure-head patterns and accompanying downwelling and upstream flow paths nested within a regional upwelling section. Sawyer et al. [2011] combined a flume and hydrodynamic model to examine the downwelling to upwelling transition around a spanning log, documenting how the hyporheic exchange rate decreased exponentially with distance from the log. Sawyer et al. [2011] also explored the reach-scale impacts by simulating periodic structures in the model. However, neither the Endreny et al. [2011] nor Sawyer et al. [2011] flume studies modeled reach-scale impacts of a sequence of restoration structures. By contrast, Mutz et al. [2007] used a flume and computational model at the reach scale (e.g., entire flume bed, multiple wavelengths), analyzing how solute concentration decay was affected by large woody debris. Despite the significant investigation into the impact of hyporheic flux around restoration and other in-channel structures, no research has examined the scaling between impacts at the local and reach scale. This research was designed to quantify how restoration structures influence local-scale and reach-scale patterns of hyporheic exchange, and to examine the how these patterns change with channel discharge.

2. Methods

2.1. Flume Experiments

[6] Flume experiments were conducted to map flow path and quantify local-scale penetration depth into the bed, and measure solute concentration decay to quantify reach-scale vertical exchange. Flume solute concentration decay experiments were used to detect the hyporheic exchange rate and to validate the reach-scale processes CFD model. The flume experiments were conducted in a 7 m long, 0.5 m deep, 0.3 m wide Armfield® S series recirculating flume with transparent walls in the James Hassett Laboratory at SUNY-ESF (Figure 1). The sediment bed had a porosity of 0.383 and an average depth of 0.3 m. The sediment was washed, well screened and had a D50 of 4.8 mm with very little variation in the size distribution (D10 ∼ 3 mm, D90 ∼ 7 mm). We used conservative tracer (sodium chloride) in the solute concentration decay experiments. The sediment was rinsed after each experiment to remove any trace of the solute used in the previous experiment. The longitudinal bed profile in the flume began with 1 m of plane bed section, followed by 5 m of 10 pool-riffle sequences (transversely uniform) with an amplitude (A) of 0.03 m and a wavelength (λ) of 0.5 m, finishing with 1 m of plane bed section. Our ratio of sediment depth to wavelength was 0.6 and based on an analysis by Cardenas and Wilson [2007] should be adequate to accommodate hyporheic exchange that does not hit against the flume's lower boundary. This bedform without-structures replicates bedform geometry used by Packman et al. [2004], a well-established reference, and the bedform matches natural channel geometry. Using Carling and Orr's [2000] empirical relation, A=a(λ2)b, our A and λ terms are within the range of a natural pool-riffle system with a = 0.076 and b =0.67. The average water depth (d) was controlled as 0.1 m in all runs, which made the riffle-riffle spacing five times the water depth (λ/d = 5), and the water depth about 3.3 times the riffle crest height (d/A ≈ 3.3). These ratios are comparable to the natural stream observations at bankfull flow conditions [Allen, 1982; Carling and Orr, 2000].

Figure 1.

Side view of the flume settings and bed geometry with and without restoration structures. Photos show the bedform and structure settings in (left) without water and (right) with water conditions, facing upstream.

[7] The flume was run with the above bedform and then with a treatment bedform of channel-spanning restoration vanes located at the top of each of the 10 riffles. The vanes were constructed using cobbles that had a D50 of 8 cm, spanning the channel perpendicular to flow and representing the center 33% of the cross-vane structure [MDE, 2000; Radspinner et al., 2010]. Each structure consisted of sill cobbles set on footers. The sills had an average height (h) of 3 cm above the bed surface and the footers were buried in the sediment with the maximum depth of 15 cm, about 5 h. Simulation Run 1 was with structures in place and Run 2 was without structures, and each run was replicated three times, removing and rinsing the sediment between runs. A discharge of 6 L/s was maintained by a centrifugal pump and flow valve and monitored with a digital flow meter mounted in the recirculating pipe. The bedform geometries in all the flume experiments were stable with no bedload sediment transport in the flume.

2.1.1. Reach-Scale Exchange (Qz) Detection

[8] A solute tracer of 200 g of NaCl solution was injected at the downstream flume reservoir at a constant rate during one flume circulation period (3 min) for all replications of Run 1 (with structures) and Run 2 (without structures). The constant rate injection rather than a rapid slug injection was used to establish a well-mixed homogenous initial surface water NaCl concentration of 0.194 g/L. The surface water volume was 1030 L for Run 1 and Run 2. The initial concentration was computed as the total mass of NaCl divided by the total volume of surface water (Vs), which was the sum of surface water in the straight flume channel, upstream and downstream reservoirs, and pipes. Seven HOBO® conductivity loggers were placed along the flume's 7 m length, each recording at a 5 s interval. HOBO loggers were placed in the upstream and downstream reservoirs, over two riffles, and over three pools. An average surface water concentration, Cs, was compiled from the temperature adjusted concentration data recorded by the seven loggers.

[9] The hyporheic exchange rate was quantified using a vertical discharge model. This model was modified based on the Mutz et al. [2007] model by adding a time-dependent variable to describe the increasing volume of hyporheic pore water, Vp. The model assumes the surface water and hyporheic pore water are each well mixed, the concentration of hyporheic pore water is initially zero, and hyporheic pore water concentration increases due to hyporheic exchange with NaCl enriched surface water. The model conserves mass and is expressed as,

CsVs+CpVp=Ceq(Vs+Vpmax)(1)

where Cs is the average concentration of NaCl in surface water, Vs is the volume of surface water, Cp is the average NaCl concentration in the pore water, Vp is the hyporheic volume of pore water, Ceq is the NaCl concentration in surface water and pore water when the mixing reaches an equlibrum state, and Vpmax is the maximum volume of pore water involved in the hyporhic exchange process. The vertical discharge can be expressed by investigating the NaCl concentration time series:

QzVp(t)VpmaxCp(t)−QzCs(t)=dCs(t)dtVs(2)

where t is the time, Qz is the vertical discharge, which is the combined downwelling of surface water and upwelling of pore water over a unit of time. The term QzVp(t)Vpmax accounts for the vertical discharge from pore water to surface water at t, creating available volume for NaCl enriched surface water. By inserting equation (1) into equation (2), the variable Vp(t) cancels and we get an equation that does not require monitoring the time rate of change in pore water volume,

where C0 is initial maximum surface water concentration, which is normalized to 1. Equation (4) describes the surface water NaCl concentration decay as function of time. Vpmax is the maximum hyporheic pore water volume,

Vpmax=VsC0Ceq−Vs(5)

[11] After inserting equation (5) to equation (4), the only unknown variable at the right side of equation (4) is the vertical discharge, Qz, which can be solved by fitting the Cs curves we obtained from the flume experiments.

2.1.2. Local-Scale Exchange (Dp) Detection

[12] Flume dye tracing experiments were used to characterize the maximum hyporheic penetration depth (Dp) with and without the structures. The results also provided supplementary evidence with CFD model verification. Nonreactive blue dye was injected at four locations along a bed surface wavelength from a pool bottom to a downstream riffle top in Run 1 (with structure) and Run 2 (Figures 2a and 2c). The dye was injected 2 cm below the bed surface to avoid the near-interface turbulence-induced diffusion [Packman et al., 2004]. The dye trajectory was delineated to record the hyporheic flow paths through the porous substrate (Figures 2b and 2d). The Dp was measured from the average level of water-sediment interface to the deepest location where dye was traced.

Figure 2.

Dye tracing experiment photos and digital delineated flow paths of dye movements (a and b) with and (c and d) without structures.

2.2. Computational Fluid Dynamics (CFD) Modeling

[13] We applied a commercial CFD modeling package Flow3D® to replicate the flume experiments and further investigate the hyporheic response of restoration structures under more extensive flow conditions. The CFD model used finite-difference methods to solve the transient velocity and pressure using conservation of mass and momentum equations for coupled surface and subsurface flows. Flow surface was treated as a free surface boundary by using the volume-of-fluid (VOF) [Hirt and Nichols, 1981] method, which computed the fraction of fluid as a variable (F) in each cell. The physical fluid properties were then weighted by the fluid fraction in the governing equations for the flow. The CFD model directly coupled surface and pore water as a seamless domain and simulated solute concentration mixing and the impacts of near-bed turbulence on flow paths. This approach has the advantage of providing a seamless model across the bed and contrasts with methods that parameterize the subsurface hyporheic flow with a bed pressure boundary condition [Cardenas and Wilson, 2007; Tonina and Buffington, 2007]. Turbulence was simulated using the renormalization group (RNG) k − ε model [Yakhot and Orszaq, 1986; Yakhot and Smith, 1992]. Similar in form to the traditional k − ε model, the RNG k − ε model applied the renormalization-group theory statistical technique to better represent the Reynolds stress and production of dissipation terms. The constants in the turbulent dissipation equation empirically found in the traditional k − ε model (e.g., the turbulent Prandtl numbers) were derived analytically in the RNG k − ε model.

[14] The CFD model has been verified for flume surface hydraulics around complex bedforms in our prior research [Zhou and Endreny, 2012]; however, we further verified the CFD for turbulent surface flow using experimental data from van Mierlo and de Ruiter's [1988] flume experiment “Run T6.” These data have been widely adopted for CFD turbulent model verifications [e.g., Cardenas and Wilson, 2007; Dimas et al., 2008; Yoon and Patel, 1996]. The Run T6 flume experiment investigated the flow surface elevation, velocity profiles, and turbulence fluctuations over a fixed, impermeable, repeating dune geometry. The CFD model represented the Run T6 experiment with a periodic boundary condition to connect the upstream and downstream ends of one bedform wavelength to represent the entire flume length. The geometry of the bedform and flow properties was imported to the CFD model following the experiment configurations (Table 1). The depth averaged velocity was controlled by an adjustable streamwise directed gravity so that the velocity equaled the flume observed experiments. The simulations were run until reaching steady state, which satisfied stability criteria with changes in total mass and mean kinematic energy less than 1%. The CFD simulation results agreed well with the Run T6 flume observations with the maximum relative error of the flow surface elevation prediction less than 0.1% and the goodness of fit (R2) of all streamwise velocity profiles over 0.85 (Figure 3).

[15] The CFD was verified for subsurface flow using the Elliott and Brooks' [1997a] “Run 9,” without subsurface circulation loops, which was consistent with our experiment. The Run 9 experiment documented the subsurface flow paths and the solute penetration fronts, in a recirculating flume with repeating fixed triangular dunes, based on dye movement. We simulated one dune with periodic boundary conditions and employed the identical bedform geometry and sediment properties in the CFD model following the Elliott and Brooks' experiment configurations (Table 1). The drag coefficient (DRG) that represented porous media flow resistance in the CFD model was calculated using the Kozeny-Carman relation based on the porosity and grain size:

DRG =180ν(1−nn)21d502(6)

where ν is the kinematic viscosity of water and n is the sediment porosity. A standard advection-dispersion model was coupled with the CFD model to further compare the CFD simulation results to the observed solute penetration fronts in the flume experiment. The simulated flow paths and solute penetration front agreed well with the observations (Figure 4). These verification results support use of the CFD model to simulate surface and subsurface flow dynamics with and without our restoration structures.

[16] The CFD model was parameterized to simulate our flume experiment Run 1 and Run 2. It represented the flume sediment surface roughness length as 0.05 m and the sediment porosity as 0.383. After a group of grid dependence test runs, the CFD model computational mesh was set to 5 mm in the streamwise direction (x), 5 mm in the vertical direction (z), and 1 grid cell in the transverse (y) direction to make this a 2-D simulation. The cross-vane restoration structure was represented by an impermeable square block half buried near the top of the riffle section and two fully buried rectangular footers underneath it. These footer blocks had an average diameter of 5 cm with a 5 mm space in between, reflecting the observed gaps in the flume settings (Figure 6a). The CFD model output reported the vertical flow velocity (Vz) at the surface of the sediment bed, with positive values indicating upwelling flows and negative values indicating downwelling flows. The Vz was sampled at a depth of two grid cells (i.e., 1 cm) below the sediment bed surface to avoid the near-bed turbulence. These local-scale Vz values were integrated along the entire flume bed to provide a CFD model simulated reach-scale hyporheic exchange rate, Qz. Assuming the upwelling area and downwelling area at the stream bed are equal, then:

Qz=12Asn|Vz|¯(7)

where |Vz|¯ is the averaged absolute vertical velocity at the streambed and As is the surface area of flume stream bed, 2.1 m2. The CFD simulated hyporheic flow velocity vectors were analyzed as streamlines to delineate actively connected hyporheic zones, and from these delineations the Dp was derived to compare with the dye tracing experiment results. In addition to the surface and subsurface flow paths verifications as presented above, the simulated variables Qz and Dp were verified with the flume observations in the result section. The CFD model was used to examine how Dp and Qz varied with changes in flume discharge. Eight additional simulations were completed for the Run 1 and Run 2 geometries with flow discharge of 1, 3, 10, and 15 L/s (Table 2) with the Reynolds number ranging 3000-50,000.

Table 2. Conditions of All the CFD Simulations

Run Number

Flow Discharge (L/s)

Structures

1Lw

1

With

1L

1

Without

3Lw

3

With

3L

3

Without

6Lw

6

With

6L

6

Without

10Lw

10

With

10L

10

Without

15Lw

15

With

15L

15

Without

3. Results

3.1. Flume Tracing Experiments

[17] The flume NaCl concentration decay experimental results indicated that restoration structures increased the hyporheic exchange rate. In the with-structures Run 1, the normalized average surface water concentration, C* (C* = Cs/C0), rapidly decayed and reached the equilibrium concentration of 0.91 at approximately 5000 s, while in the without-structures Run 2, C* had a slower concentration reduction rate and reached the equilibrium nearly three times later (Figure 5). The fitted equation (4) model for Qz provided a good match with the observed NaCl concentration decay curves from Run 1 and Run 2 with the goodness of fit (R2) of 0.984 and 0.976, respectively (Figure 5). The fitted Qz was 0.1 L/s for Run 1 (with structures) and 0.045 L/s for Run 2, indicating that the structures nearly doubled the hyporheic exchange rate. The Qz was approximately 2% of the flume discharge with structures and 1% of flume discharge without structures.

Figure 5.

Concentration decay curves recorded in flume experiments (cross plot for Run 1 with structures and dotted plot for Run 2 without structures), and fitted by vertical discharge model (gray line for Run 1 and black line for Run 2).

Figure 6.

CFD simulated hyporheic flow paths for (a) Run 1 with structures and (b) Run 2 without structures with white area showing the hyporheic zone and gray area showing the under flow zone.

[18] The flume dye tracing experimental results indicated that structures reduced the Dp and minimized the upstream-directed hyporheic flow beneath the pool area. The Dp of the delineated hyporheic flow paths in Run 1 (with structures) was 0.2 m, which was 5 cm shallower than the Run 2 flow paths (Figures 2b and 2d). The restoration structures resulted in a 20% reduction in the hyporheic zone. Dye released in the Run 1 scenario had the same upstream to downstream flow direction at all four release points, while dye released in Run 2, without the structures, had split flow directions beneath the pool. In all three replications of the Run 1 and Run 2 experiments, this difference in local-scale flow paths was observed, with reverse flow beneath the pool (Figure 2d).

3.2. CFD Modeling Results

[19] The CFD simulated Qz and Dp of Run 1 (with structures) and Run 2 were compared to flume observed results (Figures 6a and 6b). The simulated Qz for Run 1 was 0.1 L/s, which was within 5% of the flume derived Qz. In Run 2, the CFD model simulated Qz was 0.045 L/s, only 0.6% higher than the flume derived Qz. CFD simulated Dp values for Run 1 and Run 2 were 0.18 and 0.23 m, respectively, which were within 10% of the respective observed Dp values of 0.2 and 0.25 m. Considering the natural uncertainties of the flume tracing experiments (e.g., the variations of the structure cobble size, homogeneity of the bed sediment, NaCl concentration variation at the beginning of the injection, etc.), the CFD simulated variables Qz and Dp for Run 1 and Run 2 agreed with the flume derived values with an acceptable error range (<10%).

[20] The Qz with structures was greater than that without structures, while the Dp with structures was less than without structures on various discharge rates (Figure 7). The structures consistently influenced the hyporheic exchange rate and penetration depth at multiple discharge rates, even with significant variation in the water surface elevations. The Qz increased with increasing flow discharge. In simulations without structures, the Qz increased from 1.28 × 10−3 L/s to 0.31 L/s as the discharge rate increased from 1 to 15 L/s; while in simulations with structures, the Qz increased from 1.75 × 10−3 L/s to 0.5 L/s (Figures 7 and 8a), implying that in both scenarios the Qz were positively related to the discharge rate. Interestingly, Dp did not follow a similar monotonic increase with discharge. With structures in place, Dp values ranged from 0.162 to 0.183 m, and increased with Q for all flows except for the negative adjustment in Dp between 6 and 10 L/s. This nonlinear trend in changes between Dp and Q was just detectable in our simulations, with incremental changes in Dp limited to 1 grid cell depth (5 mm), and future simulations need to test if it is a significant trend. In simulations with and without structures, Dp initially increased, starting at 0.174 m with structures and at 0.204 m without structures, reached a maximum depth for simulations with 6 L/s, and then decreased or held constant by 10 L/s. The values of Dp then began increasing again for Run 1 and Run 2, as discharge increased from 10 to 15 L/s. The changing trend in Dp coincided with the emergence of hydraulic jumps, which was observed at 10 and 15 L/s with structures and at 15 L/s without structures.

Figure 7.

CFD simulated hyporheic flow directions (top) without structures and (bottom) with structures with flow discharge of 1, 3, 6, 10, and 15 L/s. The Qz and Dp values of each simulation were noted in each figure.

Figure 8.

(a) CFD simulated Qz and Dp with and without structures versus discharge rate. (b) CFD simulated Qz and Dp with and without structures versus the drop of HT.

[21] Given the total hydraulic head (HT) at the streambed is a significant driving factor for hyporheic exchange, we plotted the CFD model predicted HT for the five discharges with and without structures (Figures 9a and 9b). The “W” indicated with structure simulations. Plots of HT showed that the structures caused a significant drop in HT across the sill. As discharge increased from 1 to 15 L/s, the drop in HT increased from 0.0004 m (1%) to 0.093 m (22%) (Figure 8b). After the drop the HT gradually increased, except under a hydraulic jump (10 and 15 L/s with structures) when it experienced a sharp lift to the local maximum at 0.45λ and then decreased. Plots of Qz and Dp versus drop in HT showed Qz had a very strong positive linear relationship (R2 = 0.97) regardless of a hydraulic jump. There was not, however, a monotonic relationship between Dp and drop in HT (Figure 8b). At lower discharges, without a hydraulic jump, Dp was positively related to the head drop, and then Dp decreased with the emergence of the hydraulic jump. These results demonstrate how the structures influence the hydraulic jump, which in turn altered the HT values, redistributed hyporheic exchange patterns, and reduced flow penetration depths.

Figure 9.

The CFD simulated total hydraulic head (HT) distribution at streambed (a) without structures and (b) with structures.

4. Discussion

[22] This research demonstrated how a fully coupled surface-subsurface CFD model could represent observed components of hyporheic flow at the grid scale, local scale (e.g., bedform facet), and reach scale (e.g., entire flume). The ability for the CFD to represent the observed phenomena from grid cell to entire computational mesh enables river scientists, engineers, and managers to use a single, integrated model to assess a restoration impacts at various scales. For example, the model could examine hyporheic penetration depth and its influence on benthic macro-invertebrate habitat, redox sequences, or other biological and chemical issues. Alternatively, the model could examine the volume surface water-groundwater mixing along an entire reach to track bulk temperature or concentrations. The combined use of the flume and CFD model in the with-structures and without-structures analyses allowed for testing of a scaling relationship for hyporheic exchange. We documented that restoration structures enhance the reach-scale and local-scale vertical discharge exchange rates yet decrease hyporheic flow path depth. Furthermore, the CFD model estimates of local-scale vertical velocity (e.g., measured at every 5 mm grid cell) linearly scaled to approximate the reach-scale observed vertical discharge within 5% accuracy.

[23] A critical feature of the CFD model was its ability to simulate surface water hydraulics and associated streambed pressures, allowing it to predict how hyporheic exchange varied with changes in discharge and surface roughness around the restoration structures in the pool-riffle system. At any set discharge, the installation of in-stream restoration structures had a pronounced impact on hyporheic exchange by reducing the maximum hyporheic flow path depth of penetration Dp and by increasing the bulk rate of vertical discharge Qz between surface water and groundwater. Essentially, the exchange with structures was shallower but more intense due to the hydraulics around the structures, and Dp was regulated by two counter-acting forces of flume discharge, Q, and streambed total hydraulic head, HT. Low regions of streambed HT formed over the riffle with and without the restoration structure in place (Figure 9), which caused an upwelling trend and effectively pulled Dp upward. Counter acting to this upward influence was the volume of downwelling water upstream of the riffle, which pushed Dp downward. This volume of downwelling water was driven by relatively high water surface depths upstream of the riffle. Downstream of the riffle the pressure was controlled by water depth and near bed velocities. As discharge increased, there were increases in water surface depth upstream of the riffle (due to flow blocking and backwater effects this was more pronounced behind structures than for riffles without structures), drop in total hydraulic head, HT, across the riffle (see Figures 9a and 9b), and a monotonic increase in Qz.

[24] As predicted by a straight forward application of Darcy's Law, the Qz with structures was bigger for all Q due to the larger drop in HT. Similarly, as Q increased, the depth of penetration along the entire wavelength typically increased, where the depth is noted by the line dividing groundwater and hyporheic water (see dark gray line in Figure 7). Two features to note along this dividing line are the maximum depth, Dp, and the minimum depth, where the concave upward lines form a cone. In Figure 7, this penetration cone is located at start-of-pool for simulations without-structures and at mid-pool for simulations with structures. For any Q, Dp and the depth to the penetration cone are significantly larger without structures than with-structures. This is explained by the pronounced low streambed HT downstream of the structure causing strong upwelling at mid-pool, which vertically lifts and longitudinally shifts the penetration cone compared with the without-structure case (Figure 7). By contrast, the area below the pool in the without-structure simulation is characterized by pronounced downwelling. For conditions without the hydraulic jump, as Q increased the balance of forces upward (low streambed HT) and downward (downwelling volumes) favored a steady deepening of Dp. However, when Q transitioned from 6 to 10 L/s with-structures, a hydraulic jump formed and the streambed low HT dropped nearly threefold, over-powering the increase in downwelling and causing Dp to move upward. When the hydraulic jump intensified at Q of 15 L/s, the two forces were again in counter action; however, the 50% increase in Q overpowered the drop in streambed HT and flow depths and Dp deepened. At Q of 10 and 15 L/s, the hydraulic with-structures cause an increased complexity in hyporheic flow paths compared with flow paths without structures (Figure 7, compare bottom and top row). In the pool-riffle bedform without water surface undulations and hydraulic jumps, the variation in streambed velocity and pressure induces areas of upstream-directed upwelling from the pool toward the upstream riffle. The hydraulic jump with-structures caused a flow-path reversal in this region, as well as generated other nested areas of flow-path complexity similar to those observed by Endreny et al. [2011] when hydraulic jumps formed around drop structures.

[25]Cardenas and Wilson [2007] examined the relationships between hyporheic dynamic variables and surface flow hydraulic conditions through a series of steady state groundwater simulations. They suggested a power relationship between dimensionless hyporheic flux (q* scaled with hydraulic conductivity) and Reynolds number, as q*=aReb, and a Michaelis-Menten growth relation between dimensionless penetration depth (Dp/λ) and Reynolds number (Dp/λ=Redc+Red). It is informative to compare our normalized data with these conceptual models. We compared our flume data with these conceptual models by scaling the Qz with (K × As), where K is the sediment hydraulic conductivity, and scaling Dp with λ. In both with-structure and without-structure conditions, the scaled Qz and Re followed the power relationship, even for flows with hydraulic jumps, with a goodness of fit greater than 0.99 (Figure 10a). However, the scaled Dp only fit the Michaelis-Menten trend for flows without the hydraulic jump (Figure 10b). Given the control hydraulic jumps have on streambed HT and on Dp, it is likely a predictive equation for Dp will depend on the emergence of hydraulic jumps and depend on transitions in the Froude number from subcritical to supercritical flows.

[26] Various hyporheic exchange models have been applied in previous studies to simulate the reach-scale hyporheic processes, including the vertical discharge model [Mutz et al., 2007], effective diffusion model [O'Connor and Harvey, 2008], and the pumping model [Elliott and Brooks, 1997b]. Each of these analytical models have been proven effective in analyzing important hyporheic exchange metrics, with the pumping model simulating local-scale flux direction and magnitude but limited to sinusoidal bedforms without structures. As demonstrated in this research, the vertical discharge model is extremely useful in fitting observed solute data and providing a physical exchange term to verify a CFD model. The verified CFD model, in turn, can then examine local-scale and reach-scale exchange and its response to changes in model drivers (e.g., discharge).

[27] The effective diffusion model has been applied in reach-scale flume studies [Marion et al., 2002; Packman et al., 2004] to describe the initial rapid solute decay rate in surface water, which can be simulated in the CFD model but is not equivalent to a velocity term. The effective diffusion model uses Fick's second law to describe the hyporheic exchange processes and it combines all hyporheic exchange into the single diffusion term (De). O'Connor and Harvey [2008] used dimensional analysis to regress an analytical form of the effective diffusion model to the shear Reynolds Number and Peclet Number, and thereby advanced the opportunities for using effective diffusion experiments to verify reach-scale CFD simulations. However, a constraint of using the model to represent the restoration structures' impacts is its reliance on the initial rate of decreasing water solute concentration, Cs, not using the equilibrium concentration, Ceq. We and others have found in recirculating flume experiments, it is very difficult to maintain a monotonic decrease in Cs immediately after solute injection, even when automatically administered over the course of entire circulation cycle (A. Packman, personal communication, 2011). The Cs data rarely decays at a constant rate due to the absence of a completely mixed surface water initial condition, and the decay includes cycling increases in Cs as the center of mass of solute recirculates. The initial variation in Cs decay may represent experimental error, yet despite the sensitivity of the De term to the section of the Cs decay data, there is no objective method for selecting the most representative section of the decay curve. The vertical discharge model, however, uses an exponential function to approximate the final solute concentration, which emphasizes the most stable section of decay data and deemphasizes the initial unstable section of decay data. Another advantage of the vertical discharge model was that all the terms were fixed except Qz.

[28] We modified the Mutz et al. [2007] vertical discharge model to analyze the volume of hyporheic pore water. According to our dye tracing observations, the volume of pore water involved in the solute mixing process gradually increased as a function of time until the hyporheic zone was fully filled with the solute. Thus, we introduced a time-dependent variable Vp in equation (1) in place of the constant Vp (which was presented as Vpmax in our model) used in Mutz et al. [2007]. This modification allowed us to define the volume of actively mixed pore water and better capture the dynamics of the hyporheic exchange. This method was an alternative approach to the Ren and Packman [2004] model, which used a time-dependent variable M(t) to describe the varying solute penetration depth. Additional modifications to the vertical discharge model are possible by enhancing its time rate analysis. The model could incorporate a new term Qz(solute)/Qz, defined as the fraction of total vertical discharge that is solute passing from pores into surface water. Theoretically, Qz(solute)/Qz should be equivalent to A(t)/Amax, where A(t) denotes the streambed area with solute penetrating into the hyporheic zone and Amax is the total streambed surface area. In practice, we derived the A(t)/Amax fraction through CFD simulations on various bedforms and flow conditions. The results showed that A(t)/Amax was proportional to Vp(t)Vpmax and therefore the mass conservation equation could be written as,

dCs(t)dtVs=Qzk1Vp(t)VpmaxCp(t)−QzCs(t)(8)

and

Vp(t)=k2t(9)

where k1 and k2 are constants. The analytical solutions of equations (1), (8), and (9) contained more than 15 terms and were not considered efficient for practical applications. To achieve computational efficiency, we simplified the model algorithm by equating A(t)Amax=Vp(t)Vpmax in equation (2) to approximate the Qz(Solute)/Qz term. Even though Vp(t) canceled in the combination of equations (1) and (2), the resulting equations represent the involved pore water rather than assuming all pore water is mixed, which provides a foundation for future improvement of the vertical discharge model.

[29] While the vertical discharge model provided a useful tool for analyzing our flume experiments, there are model limitations to fitting Qz due to the underlying assumption that the solute in the involved pore water was completely mixed. Given the solute is likely moving through the porous medium as an advection-dispersion-reaction process [Cerling et al., 1990; Elliott and Brooks, 1997b], it is reasonable to consider the solute incompletely mixed in the subsurface. Follow-on work should explore other vertical exchange model factors that influence hyporheic exchange rate and maximum penetration depth. Factors of interest include (1) flow depth, (2) the bedform geometry, e.g., the number of pool-riffle sequences, pool-riffle wavelength, and the amplitude/wavelength ratio, (3) the placement and geometry of the restoration structures, including footer depth and gaps, and (4) the sediment grain size for mobile bed simulations. Variations in these factors would potentially change the pattern and magnitude of the hyporheic flow and consequently influence the degree of mixing and the reach-scale hyporheic responses.

5. Conclusions

[30] In this paper, we conducted flume experiments with solute injection and decay and dye tracing to investigate how river restoration structures impacted hyporheic response. Hyporheic response was quantified with reach-scale vertical discharge rate and local-scale maximum hyporheic penetration depth. These variables were either observed or derived from the flume experimental data and were used to verify the local-scale and reach-scale CFD simulations. The validated CFD model was then applied as a virtual flume to investigate the impacts of structures on different flow discharge conditions. The major findings of this paper are as follows:

[31] 1. The river restoration structures changed streambed pressure and created relatively high vertical flux yet shallow depths of hyporheic flow compared with the pool-riffle system without-structures.

[32] 2. The restoration structures generated backwater effects that increased downwelling forces; however, these forces were nearly counter-balanced by low bed pressures downstream of the structure. These forces explain why hyporheic penetration depth increased less rapidly with increasing flow rate for runs with structures than without structures. At a threshold discharge noted by the emergence of hydraulic jumps, the balance of these forces shifted and low bed pressures caused a shallower penetration depth despite increasing discharge.

[33] 3. The surface-subsurface coupled CFD model predicted the observed bulk solute concentrations and dye tracks by tracking flows at a 5 mm grid scale, demonstrating the power of this virtual flume model for scalable analysis of hyporheic exchange.

Acknowledgments

[34] The authors would like to acknowledge Chuck Kroll of SUNY ESF and Laura Lautz of Syracuse University for informative suggestions. Funding support was provided by the NSF through EAR 09-11547, the NY Water Resources Institute through 2011NY162B, and the USDA Forest Service Northeastern Research Station. We acknowledge Bangshuai Han, Mike Fay, Bernard Reschke, and Paul Szemkow of SUNY ESF for technical support. The constructive comments from the editor and four reviewers are appreciated.