Inﬂuence of aquifer and streambed heterogeneity on the distribution of groundwater discharge

. The spatial distribution of groundwater ﬂuxes through a streambed can be highly variable, most often resulting from a heterogeneous distribution of aquifer and streambed permeabilities along the ﬂow pathways. Using a groundwater ﬂow and heat transport model, we deﬁned four scenarios of aquifer and streambed permeability distributions to simulate and assess the impact of subsurface heterogeneity on the distribution of groundwater ﬂuxes through the streambed: (a) a homogeneous low- K streambed within a heterogeneous aquifer; (b) a heterogeneous streambed within a homogeneous aquifer; (c) a well connected heterogeneous low- K streambed within a heterogeneous aquifer; and (d) a poorly connected heterogeneous low- K streambed within a heterogeneous aquifer. The simulation results were compared with a base case scenario, in which the streambed had the same properties as the aquifer, and with observed data. The results indicated that the aquifer has a stronger inﬂuence on the distribution of groundwater ﬂuxes through the streambed than the streambed itself. However, a homogeneous low- K streambed, a case often implemented in regional-scale groundwater ﬂow models, resulted in a strong homogenization of ﬂuxes, which may have important implications for the estimation of peak mass ﬂows. The ﬂux distributions simulated with heterogeneous low- K streambeds were similar to the ﬂux distributions of the base case sce-Correspondence


Introduction
Groundwater fluxes at the interface between aquifers and streams can show strong variations in space and time at different scales (e.g., Ellis et al., 2007;Krause et al., 2007).This is important since the magnitude of groundwater discharge across the streambed influences the exchange with and the size of the hyporheic zone (Boano et al., 2008;Cardenas and Wilson, 2007) which plays a critical role for the functioning of stream ecosystems (Brunke and Gonser, 1997).For example, the exchange of water between aquifers and streams has important implications for the hydrochemistry of the streambed sediments (Malcolm et al., 2003), thus influencing biogeochemical nutrient cycling and habitat quality.A heterogeneous distribution of groundwater fluxes and hyporheic exchange flows leads to a patchy distribution of biogeochemical gradients and interstitial fauna (Boulton et al., 1998;Malcolm et al., 2004).
Spatial heterogeneities of groundwater fluxes through the streambed also impact the fate and transport of contaminants between aquifers and streams (e.g., Conant et al., 2004;Kalbus et al., 2007;Chapman et al., 2007).Schmidt et al. E. Kalbus et al.: Subsurface heterogeneity and groundwater discharge (2008) showed that the highly variable groundwater fluxes observed at a small stream resulted in a significant tailing of contaminant mass flow rates compared to the theoretical homogeneous case.
It is commonly assumed that the groundwater flux across streambeds is controlled by the heterogeneity of the connected aquifer (e.g., Wondzell and Swanson, 1996;Wroblicky et al., 1998;Storey et al., 2003;Conant, 2004).The properties of the streambed sediments may further contribute to the heterogeneous distribution of groundwater fluxes (Conant, 2004;Ryan andBoufadel, 2006, 2007).Also, geomorphologic features at different spatial scales were shown to cause variabilities of water exchange across the groundwatersurface water interface (Kasahara and Wondzell, 2003;Cardenas, 2008).Infiltrating stream water caused by streambed irregularities further leads to a very complex spatial exchange pattern (Savant et al., 1987;Salehin et al., 2004;Gooseff et al., 2006).
Our focus is on the influence of heterogeneous distributions of hydraulic conductivity (K) in the aquifer and the streambed sediments on the spatial distribution of groundwater fluxes across the streambed.In a few recent modelling studies, subsurface heterogeneity was included to simulate stream-aquifer interactions.Chen and Chen (2003) considered anisotropic and layered aquifers as well as streambeds with different hydraulic conductivities in their simulations of stream-aquifer interactions, but did not include within-layer heterogeneity.Bruen and Osman (2004) studied the effect of spatial variabilities of aquifer K on stream-aquifer seepage flow, but did not consider a separate analysis of the influence of streambed properties.Cardenas et al. (2004) simulated the impact of heterogeneous streambed deposits on hyporheic zone geometry, fluxes, and residence time distributions, but did not include the groundwater component.Fleckenstein et al. (2006) investigated the effect of aquifer heterogeneity on the distribution of seepage on an intermediate (10 2 m) scale.
The objective of our study was to investigate the potential influence of the heterogeneity of both the aquifer and the streambed sediments on the spatial distribution of fluxes through the streambed on the metre-scale.In numerical simulations we used different combinations of aquifer and streambed heterogeneity to evaluate which of these hydrological units has a stronger influence on the flux distribution.Focussing on spatial variations at fixed boundary conditions, we performed steady-state simulations to look at the flux distribution at a certain point in time.This study is a theoretical investigation of flow processes between aquifers and streams.However, we based the numerical model parameters on measured field data to obtain results in realistic orders of magnitude.

Background
Along a 220 m reach of a small stream in Germany, streambed temperatures were mapped with high resolution by Schmidt et al. (2006).The stream is a man-made stream which partially penetrates a Quaternary alluvial aquifer.It is about 3 m wide and has an average water depth of 0.6 m.The mean annual discharge is 0.2 m 3 s −1 at a gradient of 0.0008 m m −1 .The streambed consists of a 0.6 m layer of crushed rock.The interstices of the coarse crushed rock grains are filled with allochthonous, sandy, alluvial material.The connected aquifer is unconfined with a mean saturated thickness of about 8 m and consists of sandy gravel.Further information about the study site can be found in Schmidt et al. (2006Schmidt et al. ( , 2008) ) and Kalbus et al. (2007Kalbus et al. ( , 2008a)).
The streambed temperatures were mapped with a multilevel temperature probe at depths between 0.1 and 0.5 m below the streambed surface.Based on the observed temperature profiles, Schmidt et al. (2006) estimated groundwater fluxes through the streambed by applying a one-dimensional analytical solution of the heat-advection-diffusion-equation.From both the temperature observations and the flux calculations, considerable spatial heterogeneity of the groundwater discharge was observed, ranging from no discharge up to a flux of 455 L m −2 d −1 with a reach-average flux of 58.2 L m −2 d −1 .
The observed spatial heterogeneity was assumed by Kalbus et al. (2008a,b) to result from the permeability distribution of the connected aquifer.Even though observed streambed temperatures are temporally highly variable (e.g., Westhoff et al., 2007), the temperature distribution observed at a certain point in time is a consequence of the spatial distribution of subsurface permeabilities and the head and temperature gradient between groundwater and stream at the time of observation.As long as the temperature observations are recorded at a sufficient depth below the streambed surface, they are not influenced by diurnal temperature oscillations in the surface water and the system can be considered to be at a quasi-steady state for the short duration of observation (Schmidt et al., 2007).Focussing on the spatial variabilities, Kalbus et al. (2008a,b) simulated groundwater flow and heat transport through the streambed at the stream reach investigated by Schmidt et al. (2006).They included stochastically generated fields of K to represent the aquifer properties.The variance of ln(K), one parameter for the generation of the Kfields, represents the heterogeneity of the aquifer permeability.After developing a relation between this parameter and the variance of observed temperatures, it was calibrated until the simulation results matched the observed distribution of temperatures and groundwater fluxes through the streambed.From 50 realizations of K-fields used for the simulations, 10 were selected which reproduced best the field observations.Kalbus et al. (2008a,b) assumed in their simulations that the streambed had the same properties as the aquifer and thus they did not parametrize the streambed elements in the model differently from the aquifer elements.However, it is often presumed that streambed sediments are characterized by lower permeabilities due to clogging effects resulting from the deposition of fine-grained sediment and organic matter (e.g., Sophocleous et al., 1995;Su et al., 2004), siltation around macrophytes (e.g., Wharton et al., 2006), or bacterial growth and biofilms (e.g., Boulton et al., 1998;Pusch et al., 1998).These low-K layers could effect the distribution of fluxes across the streambed.Moreover, a heterogeneous streambed with a parameter distribution independent of the aquifer could lead to altered discharge patterns.Therefore, we performed subsequent simulations with different combinations of aquifer and streambed heterogeneity to identify the roles of the aquifer and the streambed in the generation of heterogeneous flux distributions.

Methodology
Based on the study by Kalbus et al. (2008a,b), we used their model set-up and the 10 K-field realizations selected in their study as the base case for subsequent simulations.To evaluate the effect of streambed characteristics, we added different hydraulic conductivity scenarios for the streambed sediments to the model.The results were compared with the base case model results and the observed distribution of groundwater fluxes obtained by Schmidt et al. (2006) from mapped streambed temperatures.

Model set-up
A two-dimensional groundwater flow and heat transport model using the model code HEATFLOW (Molson et al., 1992;Molson and Frind, 2005) was set up according to the model used by Kalbus et al. (2008a,b).The conceptual model (Fig. 1) represents a vertical longitudinal profile along the streambed and within the underlying aquifer, and corresponds to the length of the investigated stream section and the saturated thickness of the aquifer (220 m×8 m).The upper 0.60 m hydrostratigraphic layer represents the streambed sediments, which corresponds to the thickness of the crushed rock layer at the study site.The model grid consists of 220×65 elements with a layer thickness varying from 0.20 m at the bottom to 0.05 m at the top.Since we look at spatial variations at a certain point in time, the system is assumed to be at steady state.The bottom and top boundaries are constant head boundaries, left and right boundaries are no-flow boundaries, leading to vertical flow through the system.Although the assumption of vertical groundwater discharge seems rigid for complex stream-aquifer systems it is commonly made for the interpretation of groundwater fluxes through the streambed (e.g., Cardenas and Wilson, 2007;Keery et al., 2007).The constant head values were chosen such that for each simulation the mean groundwater flux through the model equalled the mean flux through the streambed calculated from the observed temperature profiles (q z mean=58.2L m −2 d −1 ).The temperature boundary conditions correspond to the mean stream water temperature during the mapping programme (18.4 • C) at the top boundary and the constant deep groundwater temperature (10.9 • C) at the bottom boundary.No energy flux is assumed across the left and right boundaries, because in conditions of vertical flow the lateral heat transport by conduction is negligible.The thermal transport properties were taken from the literature (thermal conductivity of the saturated sediments=2 J s ).The thermal conductivity of saturated sediments varies only little between different sediment types and can therefore be reliably estimated from literature data (Stonestrom and Blasch , 2003).A porosity of 0.25 was estimated from field data.
A heterogeneous distribution of hydraulic conductivity was achieved by including stochastically generated fields of hydraulic conductivity in the simulations.With the code FGEN (Robin et al., 1993), the K-fields were generated from the mean and variance of ln(K) and the correlation lengths in each direction (Table 1).These data were obtained from field observations of K, except the variance which was calibrated with the observed temperature distribution by Kalbus et al. (2008a,b).Ten realizations of the K distribution were used for the simulation of each of the scenarios explained below.The heterogeneous K-fields of the aquifer were identical in all scenarios and the base case.

Scenarios
The case simulated by Kalbus et al. (2008a,b) was taken as the base case for comparison with the four streambed scenarios described in the following list.For the base case, the streambed was assumed to be part of the aquifer and have exactly the same permeability characteristics.The K-fields were generated for the entire model domain.Cases A and B were chosen to isolate the influence of the aquifer and the streambed, respectively.Cases C and D were defined to include effects of clogging together with different concepts www.hydrol-earth-syst-sci.net/13/69/2009/ Hydrol.Earth Syst.Sci., 13, 69-77, 2009 of streambeds: one considering the streambed as part of the aquifer (Case C) and one considering the streambed as a separate unit (Case D).
Case A: The aquifer was assumed heterogeneous as in the base case, the streambed was assumed homogeneous with K two orders of magnitude less than the mean aquifer K.This scenario was selected to demonstrate the effect of a low-K stream boundary condition using a conductance term as it is often used in regional-scale groundwater flow models.The rather small value of the streambed K was selected to represent a clogged layer and cause clear effects.
Case B: To investigate the potential of the streambed sediment layer alone to cause a heterogeneous distribution of groundwater discharge, the aquifer was assumed homogeneous with the same mean K as in the base case and the streambed was assumed heterogeneous with the same statistical properties as the aquifer in the base case.
Case C represents a naturally developed streambed which basically consists of the same material as the underlying aquifer, but is assumed to have experienced some clogging.The aquifer and streambed were both assumed heterogeneous (using the same variance and correlation lengths as in the base case), while streambed clogging was simulated by dividing the K value of each streambed element by 100.The streambed thus has the same degree of heterogeneity as the aquifer, but the mean K is two orders of magnitude less.
Case D: The streambed properties were assumed to be independent of the aquifer, which may occur for instance in streambeds with high sediment turnover rates or in man-made streams.The connectivity between aquifer and streambed is lower than in Case C, which was achieved by generating new K-fields for the streambed layers only.As in Case C, the mean streambed K was chosen two orders of magnitude less than the mean aquifer K.The other statistical parameters for the K-field generation (variance, correlation lengths) were adopted from the aquifer statistics to enable a direct comparison with Case C.
The aquifer and streambed properties used for the generation of K-fields for the simulations of the base case and the four scenarios are summarized in Table 1.

Results and discussion
The groundwater fluxes simulated in the base case scenario are highly variable with a standard deviation of σ (q)=63.7 L m −2 d −1 matching well the observed variation (σ (q)=65.5 L m −2 d −1 ).Homogeneous low-K streambeds (Case A) significantly dampen the groundwater fluxes compared to the base case scenario and result in a relatively uniform spatial pattern of fluxes with a small standard deviation of σ (q)=6.4L m −2 d −1 (Fig. 2a).The range of fluxes is much smaller than in the base case (Fig. 3).Homogeneous low-K streambeds thus serve as homogenizing layers which reduce the influence of the aquifer texture.It is highly unlikely that this case occurs in reality, since all naturally developed streambeds as well as artificially constructed streambeds develop some degree of heterogeneity resulting from groundwater fluxes, sediment turnover, hyporheic fluxes, or activities of the interstitial and benthic fauna.
Nevertheless, homogeneous low-K streambeds are often implemented in regional-scale groundwater flow models (e.g., McDonald and Harbaugh, 1988) and in the analysis of stream flow depletion through pumping (Chen et al., 2008) where the stream-aquifer interaction is governed by a conductance term representing the resistance of the streambed (Rushton, 2007).
This approach may be sufficient for evaluating average water budgets on a regional scale, but prediction accuracies may be low because it would not provide a range of possible fluxes.For a detailed analysis of flow and transport processes, a homogeneous representation of the streambed may not be appropriate.For instance, in cases of contaminated groundwater discharging to a stream, maximum contaminant mass flow rates may be underestimated since areas of high groundwater discharge contribute more mass flow than low-discharge areas.Schmidt et al. (2008) also showed that a heterogeneous distribution of groundwater discharge strongly influences the time scales of contaminant release from a contaminated streambed.Hence, for small-scale investigations of stream-aquifer interactions, a representation of the streambed in flow models using a boundary condition with a uniform conductance term is not recommended.The streambed conductance should rather be resolved on a small scale to cover the range of high-and low-permeability zones and thus the range of high and low groundwater fluxes in the streambed.
In Case B, a heterogeneous streambed on top of a homogeneous aquifer leads to a wider distribution of fluxes than in Case A (Fig. 2b) with a standard deviation of σ (q)=22.5 L m −2 d −1 , but the range is still much smaller than in the base case (Fig. 3).This case is also highly unlikely to occur in reality, since all aquifers show some degree of heterogeneity.Nevertheless, it shows that the streambed alone does not cause the observed distribution of fluxes, at least not for the considered ranges and patterns.We performed some simulations increasing the variance of ln(K) of the streambed K-fields to see how large it would have to be to cause a flux distribution similar to the base case (data not shown), but we did not get close to the base case flux distribution within a reasonable range of variances.Gelhar (1993) gave a range of variances from 0.16 to 4.41 for alluvial aquifers.In our simulations, even with a variance of 10, which is a highly unreasonable value, the simulated range of fluxes was still too small (σ (q)=33.6L m −2 d −1 ).The passage through the streambed, which is much shorter compared to the passage through the aquifer, seems insufficient to cause highly diverse flow paths.Larger structures are necessary to direct the flow into highly permeable zones resulting in higher flow velocities.A heterogeneous streambed with a mean K two orders of magnitude less than the mean K of the heterogeneous aquifer (Case C) shows a similar pattern of fluxes to the base case (Fig. 2c).The high-and low-discharge zones are at the same locations and the range of fluxes is similar to the range of the base case (Fig. 3).The maximum fluxes are even higher than those of the base case (σ (q)=89.0L m −2 d −1 ).This is a result of the larger gradient which had to be www.hydrol-earth-syst-sci.net/13/69/2009/ Hydrol.Earth Syst.Sci., 13, 69-77, 2009 implemented in the models to achieve the reach-average flux of 58.2 L m −2 d −1 (average hydraulic gradient=0.011;base case: 0.002).Within high-permeability zones, this higher gradient leads to increased fluxes compared to the base case with a lower gradient.When reaching the streambed, the short passage through the less permeable streambed does not have much influence on the flow velocities in these zones since the permeability is still higher than in the neighbouring low-discharge zones.
In case of an independent heterogeneity of the streambed (Case D), the pattern is still similar to that of a related heterogeneity as in Case C, but the locations of high-and low-discharge zones have been slightly displaced, some peaks have disappeared, while other peaks have developed (Fig. 2d).The range of fluxes is almost identical with the range of the base case (Fig. 3) with a standard deviation of σ (q)=74.8L m −2 d −1 .Again, the higher gradient (average hydraulic gradient=0.014)leads to increased flow velocities through the high-permeability zones of the aquifer.As opposed to Case C, however, groundwater flow from high-K zones within the aquifer may now intersect low-permeability zones in the streambed and will thus be diverted to neighbouring zones with higher permeabilities.This attenuates some of the peak flows observed in Case C and creates new peaks at other locations.
Comparing the mean (solid line) and median (dashed line) in Fig. 3, it becomes apparent that greater spatial heterogeneity mainly leads to an increase in the proportion of high fluxes.Because we assumed vertical flow through the model domain, the fluxes cannot become less than zero, but well connected high-permeability zones can lead to very high fluxes which are concentrated in small areas.This is even more evident from Fig. 4, which shows the relative contribution of the streambed area to the cumulative flux.In Cases A and B, the band representing the range between maximum and minimum fluxes of all K-field realizations is narrow and almost straight with a slope of 1:1.In these cases, a certain proportion of streambed area thus contributes a similar proportion of cumulative flux.For instance, 20% of the streambed area contributes 22% (Case A) to 30-33% (Case B) of the cumulative flux.In Cases C and D, a much smaller proportion of streambed area contributes a larger proportion of cumulative flux.For instance, in Case C, 20% of the streambed area contributes 50-74% of the cumulative flux along the modelled reach.The band is much wider in Cases C and D, indicating considerable variation between the different K-field realizations.In Case C, the variations are even larger than in Case D, which is a result of a better connectivity between aquifer and streambed.We used the transport connectivity indicator CT 1 by Knudby and Carrera (2005) to analyse the relation between connectivity and discharge variations.CT 1 was defined by Knudby and Carrera (2005) as the ratio between the average arrival time t AVE of a solute travelling through the model domain and the time t 5 at which 5% of the solute has arrived.Better connectivity leads to shorter early arrival times.Here we use the reciprocal, 1/CT 1 =t 5 /t AVE , and thus better connectivity results in values of 1/CT 1 approaching zero.
t AVE was obtained from the domain length in flow direction and the mean flux across the domain.t 5 was determined through particle tracking.5000 particles were released at the bottom of the model domain and t 5 was determined from the breakthrough curves.
Figure 5 shows the standard deviation of the flux across the streambed σ (q) in relation to the connectivity indicator 1/CT 1 .It is apparent that a higher connectivity results in a higher variation of discharge rates.In Case C, smaller values of 1/CT 1 are found than in Case D, which confirms that the larger variations in fluxes observed in Fig. 4 (Case C) are a result of better connectivities between aquifer and streambed.

Conclusions
Previous simulations of groundwater flow and heat transport through a streambed have revealed that strong spatial variations in groundwater discharge to a stream are caused by a heterogeneous distribution of aquifer hydraulic conductivity.The influence of the streambed on the distribution of fluxes was investigated in subsequent simulations with different scenarios of aquifer and streambed hydraulic conductivity.The aquifer was found to have a stronger influence on the spatial distribution of fluxes than the streambed.However, the implementation of a homogeneous low-K streambed within a heterogeneous aquifer caused a significant homogenization of the fluxes.This behaviour should be considered when using the concept of streambed conductance in regional-scale groundwater models.A heterogeneous distribution of hydraulic conductivity only in the streambed was not sufficient to cause strong flux variations.Simulation results with heterogeneous low-K streambeds were similar to the results from the model without a distinction between aquifer and streambed properties.Thus, if streambed clogging, which leads to a reduced permeability of the streambed sediments compared to the aquifer, has to be considered in a model, it might be appropriate to implement a heterogeneous distribution of streambed hydraulic conductivity even at large scales to avoid an underestimation of peak flows.These results also confirm the applicability of the methodology proposed by Kalbus et al. (2008a,b) to use measured streambed temperatures for calibration of aquifer properties even without distinguishing between the aquifer and streambed.
Observed distributions of groundwater fluxes through the streambed may often be a result of both aquifer and streambed heterogeneity, with the aquifer having a stronger influence.Numerical model predictions of groundwater flow and solute transport may thus significantly benefit from heterogeneous distributions of aquifer and streambed properties.Since mass fluxes of dissolved compounds across streambeds are governed by the flux of water, the consideration of heterogeneous flux distributions is essential for the prediction of contaminant transport and for biogeochemical modelling at the groundwater-surface water interface.

Fig. 2 .
Fig. 2. Observed (top left; afterSchmidt et al., 2006) and simulated (base case and Cases A-D) results showing temperature (colour maps) and flux distributions (white curves) in the streambed (represented by the upper grey zone in Fig.1).Temperature data are shown at streambed depths between 0.1 and 0.5 m corresponding to the observations.Simulated results are shown from one example out of ten K-field realizations (the same realization is shown in all scenarios).Vertical exaggeration is approx.100×.

Fig. 3 .
Fig. 3.Box plots of the groundwater discharge through the streambed showing 95th and 5th percentile (dots), 90th and 10th percentile (error bars), 75th and 25th percentile (box), arithmetic mean (solid line), and median (dashed line).Observed data are complete data of the mapping programme (n=140), simulated data are the complete data set from all 10 realizations (n=2200) for each case.

EFig. 4 .
Fig. 4. Distribution of groundwater fluxes through the streambed in relation to the streambed area.Bands show the full range between maximum and minimum values of observations and modelling results, respectively.

Fig. 5 .
Fig. 5. Standard deviation of groundwater flux through the streambed, σ (q), in relation to the connectivity indicator 1/CT 1 for Case C and Case D.

Table 1 .
Aquifer and streambed properties of all simulation cases.K=hydraulic conductivity, σ 2 =variance of ln(K), λ x and λ z =correlation lengths in the xand z-directions.