Next Article in Journal
A Comprehensive Evaluation Model of Ammonia Pollution Trends in a Groundwater Source Area along a River in Residential Areas
Next Article in Special Issue
Exploring Local Riverbank Sediment Controls on the Occurrence of Preferential Groundwater Discharge Points
Previous Article in Journal
Gendered Perspectives on Climate Change Adaptation: A Quest for Social Sustainability in Badlagaree Village, Bangladesh
Previous Article in Special Issue
Comparative Analysis of Runoff and Evaporation Assessment Methods to Evaluate Wetland–Groundwater Interaction in Mediterranean Evaporitic-Karst Aquatic Ecosystem
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Water
Volume 13
Issue 14
10.3390/w13141923
water-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Numerical Representation of Groundwater-Surface Water Exchange and the Effect on Streamflow Contribution Estimates

by
Sachin Karan
1,*,
Martin Jacobsen
2,3,
Jolanta Kazmierczak
1,
José A. Reyna-Gutiérrez
4,
Thomas Breum
5 and
Peter Engesgaard
3
1
Department of Geochemistry, Geological Survey of Denmark and Greenland (GEUS), 1350 Copenhagen, Denmark
2
Rambøll Denmark, 2300 Copenhagen, Denmark
3
Department of Geoscience and Natural Resource Management, University of Copenhagen, 1350 Copenhagen, Denmark
4
Ørsted, 2820 Gentofte, Denmark
5
GEO, 2800 Kgs. Lyngby, Denmark
*
Author to whom correspondence should be addressed.
Water 2021, 13(14), 1923; https://doi.org/10.3390/w13141923
Submission received: 10 June 2021 / Revised: 1 July 2021 / Accepted: 8 July 2021 / Published: 12 July 2021
(This article belongs to the Special Issue Advances in the Study and Understanding of Groundwater Discharge to Surface Water)
Browse Figures
Versions Notes

Abstract

:
The effects of streams and drainage representation in 3D numerical catchment scale models on estimated streamflow contribution were investigated. MODFLOW-USG was used to represent complex geology and a stream network with two different conceptualizations—one with equal cell discretization in the entire model domain and another with refined cell discretization along stream reaches. Both models were calibrated against a large data set including hydraulic heads and streamflow measurements. Though the optimized hydraulic parameters and statistical performance of both model conceptualizations were comparable, their estimated streamflow contribution differed substantially. In the conceptualization with equal cell discretization, the drainage contribution to the streamflow was 13% compared to 41% in the conceptualization with refined cell discretization. The increase in drainage contribution to streamflow was attributed to the increase in drainage area in proximity to the stream reaches arising from the refined discretization. e.g., the cell refinement along stream reaches reduced the area occupied by stream cells allowing for increased drain area adjacent to the stream reaches. As such, an increase in drainage area equivalent to 7% yielded a 146% increase in drainage contribution to streamflow. In-stream field measurements of groundwater-surface water exchange fluxes that were qualitatively compared to calculated fluxes from the models indicated that estimates from the refined model discretization were more representative. Hence, the results of this study accentuate the importance of being able to represent stream and drain flow contribution correctly, that is, to achieve representative exchange fluxes that are crucial in simulating groundwater–surface water exchange of both flow and solute transport in catchment scale modeling. To that end, the in-stream measurements of exchange fluxes showed the potential to serve as a proxy to numerically estimate drainage contribution that is not readily available at the catchment scale.

1. Introduction

In catchment scale hydrological modeling, the numerical discretization needed to cover large areas, often in scales > 100 km2, instigates challenges of representing field measurements conducted on the meter scale. For groundwater (GW)–surface water (SW) exchanges, in-stream and streambed measurements provide valuable information for testing large-scale models that quantify both the magnitude and direction of exchange fluxes and the variability within the meter scale along stream reaches [1,2,3,4]. Such measurements are often not utilized in catchment scale modeling [5], for good reasons. Measured variability in GW-SW exchanges along stream reaches is difficult to represent in models where the spatial discretization greatly exceeds the scale of the measured variability. Hence, numerical studies using field measurements on the meter scale are commonly restricted to sub-catchment areas with finer model discretization [6,7,8,9]. Although successful in representing field measurements on the meter scale, sub-catchment models are inadequate for analyzing the impacts of management decisions on, for example, reducing nutrient loading from agricultural land use, which necessarily must be taken on the regional catchment scale.
Artificial drainage such as tile drains and drain ditches act as headwaters of streams and are often sources of nutrient loading from agricultural usage to downstream catchments [10,11,12]. Tile drainage to sustain crop production in poorly drained soils is a large recipient of precipitation. For instance, 20-years monitoring data from the Danish Pesticide Leaching Assessment Program show that the annual mean ratio of drainage to precipitation varies from 13 to 39% at three different agricultural sites [13]. In a 4 km2 agricultural headwater catchment, King et al. [14] showed that 21% of the annual precipitation was captured by drainage, and using regression analysis of measured drainage and precipitation they found that annual mean drainage accounted for 41% of the streamflow. In a numerical analysis, Hansen et al. [15] calibrated a small area hydrological model (1 km2) with detailed information of the drain network against observed drainage, streamflow, and hydraulic heads. By testing the calibrated model in a setup without the drain network, they found that groundwater discharge to the stream increased from 28% to 45%. Steiness et al. [16] used multiple field measuring techniques in a lowland riparian area and estimated that 67% of measured streamflow was groundwater-fed surface runoff prompted by anthropogenic changes related to ditches and drainage systems in the wetland next to the stream. The impact of artificial drainage on the hydrological processes in the near-stream environment is evident, and consequently, the impact of artificial drainage on nutrient loading to SW is often substantial [17,18]. Thus, there is a need to improve the representation and quantification of drainage contribution to SW on the catchment scale [11,12,14,18].
A few studies have investigated the efficiency in how drains are conceptualized in 3D numerical models (e.g., [19]), and to our knowledge, no studies have investigated the impact on drainage contribution to streams arising from how the model is discretized along stream reaches. Though some studies argue that discrete representation of drainage is best represented directly in 3D numerical models, they also conclude that such implementation in large-scale catchment models would yield inefficient due to disproportional long computational times and lack of data [15,19]. Further, the incorporation of discrete drain networks would require drainage measurements to assess the validity, which is impossible to achieve at the catchment scale [18]. Thus, assessing the representativity of simulated drainage contribution to streamflow in catchment scale models without measurements of drainage remains a challenge. A common approach in 3D catchment modeling is to incorporate drains as a head-dependent boundary in the entire model domain accounting for tile drains, drain ditches, and other smaller tributaries not discretely defined (e.g., [20,21]). However, as small streams of 10 m or less in width make up the majority of streams within many catchments and receive the most water and dissolved nutrients [22], a coarse numerical discretization applied in the catchment scale exceeding the stream width may impair a correct simulation of GW-SW exchange. Consequently, using such a model for understanding non-point pollutant transport dynamics in relation to GW-SW exchange could be incorrect due to the innate deficiencies.
The objective of this study was to implement MODFLOW-USG [23] for a headwater catchment—from which a wealth of geological and hydrological data exists and where numerous investigations of GW-SW exchange have been carried in the past decade—to identify the numerical controls on flow paths to streams. Specifically, we aim to use MODFLOW-USG capabilities, allowing inclusion of complex geology and refined cell discretization in the vicinity of streams in a headwater catchment to investigate: (i) the effect of numerical discretization on the flow balance, and (ii) GW-SW exchange estimates along stream reaches. The MODFLOW-USG results are compared with common modeling approaches using regular discretization in the entire modeling domain.

2. Materials and Methods

2.1. Study Area and Setting

The study area is situated in western Denmark and comprises a 108 km2 headwater catchment including two lakes and 52 km streams. Numerous field investigations and local scale modeling studies of GW-SW exchanges have been conducted within the catchment in the past decade (e.g., [3,4,6,16,24,25,26]), and a wealth of borehole data with geological information and measured water levels exist in the Danish National well database, JUPITER [27]. Based on water level observations from JUPITER and measured water levels from the many field investigations, the groundwater catchment was established from average head observations (Figure 1). The catchment consists of marine Miocene silt and sand layers overlain by Quaternary clay, sand, and gravel deposited during three glacial-interglacial periods [28,29]. The main stationary line of Weichselian glaciation runs through the catchment [29], dividing the landscape into flat outwash plain at the western part and moraine hills in the south-eastern area. Streamflow is from east to west. The main stream, Holtum Stream, is a second-order 15-km-long perennial stream with increasing width near the headwaters from around 0.4 m to around 8 m at the catchment outlet. From the hydraulic head contours, it is evident that Holtum Stream is a gaining stream with the largest gradients towards the stream occurring at the downstream reaches.

2.2. Geological Setting

Geological model was built using GeoScene 3D and ArcGIS. Lateral extension and top elevation of the layers were interpreted in GeoScene 3D based on the data from 3664 boreholes from the JUPITER database, 20 electrical resistivity tomography (ERT) profiles, two SkyTem data sets, Danish soil map, and Digital Elevation Model (Figure 2). The highest borehole density was in the western part of the groundwater catchment. Seventy-three percent of the boreholes had a depth <20 m and 0.01% had a depth >100 m with the deepest borehole (334 m) located south-east of Lake Hampen. ERT profiles were measured in the area of Lake Ejstrup and Lake Hampen. SkyTem data sets covered the eastern and central parts of the catchments. The top elevation of the layers was interpolated to a 50 × 50 m grid in ArcGIS following the shape of geological structures and discontinuities in the layer extensions. The interpolated elevations were adjusted to ensure no mutual crossing of the layers.

2.3. Numerical Modeling

2.3.1. Conceptual Model

Simulations were performed with MODFLOW-USG [23] and two steady-state model conceptualizations were used; a horizontally regular grid (REG) with a discretization of 100 × 100 m and a horizontally irregular grid (IRREG) with a base size of 10 m along with all stream segments. As such, the IRREG had a numerical discretization of 6.25 × 6.25 m along the stream cells increasing to 100 × 100 m away from the stream. Flow in the streams was simulated using the SFR2 package [30], where stream width and elevation were specified according to data from Karan et al. [31]. The streambed thickness is assumed to be 1 m. A general head boundary (GHB) was used to define the lakes and time-averaged measured lake stages were specified as head stage. Drains were assigned in the remaining modeling area in a depth of 1 m to account for drainage from agricultural land use (around 54%), riparian lowlands, as well as ditches and smaller tributaries not represented separately due to lacking information. This is common practice in many regional scale models (e.g., [20,21]), and drainage is not activated unless the simulated water level exceeds the specified drain elevation. Unsaturated zone processes were not in focus in the current study, and recharge was based on average streamflow measurements from the gauging station near the stream outlet (Figure 1). The recharge flux was calculated from the average streamflow divided by the catchment area. It is thus assumed that the average streamflow represents steady-state flow conditions and that all recharge is captured by the stream (and drains to streams). Hence, the outer boundary of the model domain was specified as no-flow.

2.3.2. Model Calibration

The parameter estimation code PEST [32] was used for calibration of the models, and these were calibrated against 419 hydraulic head data representing seven aquifer units. The hydraulic head data were a mix of single measurements and the average time series of head observations. The uncertainty of the hydraulic head data was calculated based on procedures in Henriksen et al. [33] and yielded a standard deviation of around 3 m. As the clay units of Quaternary and Miocene origin are not present throughout the model domain, the sand aquifers are not clearly separated. For simplicity, the hydraulic head data were treated equally, i.e., the weighting of the data within each of the six sandy geological layers representing aquifers was the same. The average streamflow was also used in the calibration and the weight was specified to assure that the contribution to the total objective function was similar for both hydraulic head and streamflow data. The average streamflow is based on observations from 2007 to 2016. The streamflow measurements have an estimated uncertainty of 10%. The average streamflow was compared to the sum of simulated streamflow and drainage. In MODFLOW, the simulated drainage (using the drain package) is treated as water leaving the model domain, and therefore the drain contribution must be added to the simulated streamflow when comparing to observed flow.

3. Results

3.1. Geology

The catchment comprises Miocene and Quaternary gravel, sand, silt, and clay layers (Figure 2). The bottom of the geological model is at the top of a fine grained, low permeable Vejle Fjord Formation. The Miocene sequence consists of alternating silt and sand layers locally overlain by Maarbaek clay formation. The total thickness of the Miocene deposits increases towards the east and reaches up to 266 m. Miocene formations were partially eroded during glaciations and a buried valley structure reaching −120 m below sea level was created at the eastern boundary of the catchment. The buried valley is mostly filled with sand and gravel covered by a thick clay layer. Miocene and the buried valley deposits are overlain by Quaternary sand and gravel with clay lenses. The thickness of the Quaternary sediments outside of the buried valley structure is up to 100 m. The uppermost Quaternary layer in the major part of the catchment consists of sand and gravel deposited at an outwash plain. In the south-eastern part of the catchment, east from the main stationary line, the sedimentary sequence is covered with moraine clay forming a hilly landscape.

3.2. Model Calibration

Each of the interpreted geological units was assigned a horizontal hydraulic conductivity and an anisotropy ratio for the horizontal and vertical hydraulic conductivity (Kh/Kv). Since the aquifers were not clearly separated (Figure 3), the parameters for the Quaternary and Miocene sand and clay layers were tied. This left four geological units to be included in the calibration. The stream network was divided to represent one streambed vertical hydraulic conductivity (Kv streambed) for Holtum Stream and one for all its tributaries. PEST was initially used to find the most sensitive parameters plotted as relative composite sensitivities in Figure 4. A sensitive parameter is here defined as having a sensitivity ≥ 3% of the maximum sensitivity. During the sensitivity analyses, different initial values were tested, which showed that with high conductance values for drains and lake GHB (generally > 1 1/d), the parameter became insensitive. Due to correlation between the two lake GHB conductances, these were also tied during the calibration. Likewise, the anisotropy ratio was fixed to 10, because of the correlation to Kh. For both model conceptualizations, the Kh Quaternary clays were insensitive and therefore not included in the calibration and set to 0.1 m/d similar to values previously used to simulate Quaternary clays in the area [6]. A lower Kh for the Quaternary clays was tested and although the model converged successfully, the Kh value of 0.1 m/d was maintained due to longer model run times with lower values, which was deemed acceptable as the Kh Quaternary clays were insensitive.
During the calibrations, it was evident that for the IRREG model, Kv for the streambed for the tributaries could not be uniquely optimized, resulting in large confidence intervals. Hence, for simplicity, in both model conceptualizations, Kv for the streambed was optimized for a single parameter representing tributaries and Holtum. However, as final the calibration results were deemed satisfactory based on calculated root mean square error (RMSE), coefficient of determination (R2), and mean absolute error (MAE), the selected approach was maintained.
For both the REG and IREG models, root mean square error (RMSE) for hydraulic heads was well below the estimated uncertainty of 3 m and coefficient of determination (R2) close to one (Figure 5). The two models represent equally well the mean hydraulic heads over the range from 50 to 117 m. This was also supported by the mean absolute error (MAE) of 1.80 m and 1.77 m for the REG and IRREG models, respectively. For both models, the simulated streamflow was within 0.1% and thus well within the uncertainty of 10% for the flow observation.
The optimized parameters for the models yield comparable values for the REG and IRREG models (Table 1). The calibrated hydraulic conductivities for geological layers and the streambed were well constrained with relatively narrow confidence intervals. For the lake GHB’s, the calibrated conductance was only well constrained in the IRREG model, while the GHB- and drain conductances in the REG model were not. Several calibrations were performed with different initial values for both conceptualizations and weights adjusted to secure equal contribution of the observation groups to the total objective function. Still, it was not possible to obtain a well constrained parameter value for the drain conductance in the REG model. Since the drain conductance was insensitive with values >1 1/d and the drain conductance in the IRREG model calibration kept reaching the upper bound of the specified parameter range (which was also tested), it was decided to fix the parameter in the model. Therefore, based on the final REG model calibration of the drain conductance, a fixed drain conductance of 0.05 1/d was used in the calibration of the IRREG model (Table 1).

3.3. Stream Flow Contribution and Groundwater-Surface Water Exchange

The flow balance relative to the total outflow showed that the main difference between the two conceptualizations was related to the simulated discharge to streams and drainage to streams (Table 2). As the total estimated streamflow was represented by the sum of simulated streamflow and drainage, the drain contribution to the total streamflow was equivalent to 13% and 31% for the REG- and IRREG model, respectively.
Simulated groundwater fluxes to Holtum stream were extracted from the model results and compared to estimated groundwater discharge rates from previous field measurements (Figure 6). For the REG model, the simulated groundwater discharge ranged from 0.3 to 3.5 m/d with an average of 1.41 m/d, while in the IRREG model, the groundwater discharge ranged from 0.01 to 2.6 m/d with an average of 0.9 m/d. Hence, the groundwater discharge to the stream in the REG model was substantially larger than in the IRREG model. The average of GW-SW exchange fluxes from field measurements was 0.41 m/d and ranged from 0.06 to 1.73 m/d mainly around 3 km of the most downstream part of the stream. Overall, both models showed decreasing groundwater discharge moving downstream. At the upstream part, from around 1700 m to 6500 m, where the largest groundwater discharge was calculated, there was a pattern of areas being flooded in both model conceptualizations. The sudden drops in calculated groundwater discharge in both models, for instance at around 5000 m, coincided with stream cells defined in areas represented by clay.

4. Discussions

4.1. Parameter Calibration

The calibrated Kh for the Quaternary sands was higher compared to the Miocene sands with a ratio of approximately 2:1. Kidmose et al. [6] found a similar ratio in a model investigating the GW-SW exchange at Lake Hampen, which is situated in the north-eastern part of the catchment (Figure 1). Kidmose et al. [6], however, attained slightly higher Kh values for the Quaternary sand of around 70 m/d. The calibrated Kh of Quaternary sand is in good agreement with measured Kh of surface-near aquifers in the riparian zones. Karan et al. [7] reported values ranging from 0.2 to 19.6 m/d, Karan et al. [25] values ranging from 16 to 69 m/d, and Steiness et al. [16] values ranging 0.2 to 50 m/d. Further, the calibrated Kh for Quaternary sand was well constrained with relatively narrow confidence intervals (Table 1). In addition, different initial values were tested for the parameters in the calibrations to ensure that no local minima were encountered, and the optimizations yielded similar values for the calibrated parameters (not shown). Thus, the calibrated Kh for Quaternary sand is deemed representative. Though no measured values of Kh for the Miocene sand exist for the area, the ratio correspondence to previous studies [6] and the well constrained parameter value from the calibrations (Table 1) indicate Kh to be representative for Miocene sands. Likewise, the calibrated Kh of Miocene clays were also well constrained and comparable to previously calibrated values for the Miocene clay in the area [6]. Further, the calibrated Kh of Miocene clays supports the assumed 0.1 m/d for Quaternary clays.
The calibrated Kv for the streambed was similar in the two model conceptualizations with well constrained confidence intervals (Table 1). Compared to previous studies [2,4] of the streambed hydraulic conductivities in Holtum stream, the calibrated values were comparable although in the lower range. In a detailed investigation of the horizontal- and vertical streambed hydraulic conductivities, Sebok et al. [4] estimated Kv streambed in the range of 0.01 to 8.4 m/d. Karan et al. [2] used measured horizontal streambed hydraulic conductivities in a numerical estimation of the anisotropy ratio. Based on these ratios, the estimated Kv for the streambed would range from 0.14 to 13.6 m/d. As shown by Sebok et al. [4], Kv can be subject to spatial and temporal variations (scouring and sedimentation). Here, the entire stream network with Holtum stream encompassing 15 km was represented with a single Kv. Therefore, a calibrated Kv value in the lower range of previously estimated values is likely representative in the applied conceptualizations. A more sophisticated approach of dividing the stream into multiple segments each having its own Kv and varying the weighting schemes of the observations to contribute differently to the total objective function could potentially lead to a better optimization with well constrained Kv for the individual stream reaches. However, as the final calibration results were deemed satisfactory based on calculated RMSE, R2, and MAE (Figure 5), the selected approach was maintained. In future work, the spatially distributed Kv will be further investigated by dividing the stream network into multiple reaches to enable calibration of Kv streambed at the reach level rather than in the entire lengths of the stream.
The conductances for drain and GHB were not well constrained in the REG model and although several initial values were tested, it was not possible to reach unique estimates of the parameters. For the IRREG model, where the drain conductance was fixed, the calibrated lake GHB conductance was well constrained (Table 1), and thus, there seemed to be a correlation between the two parameters. Nevertheless, the calibrated values for drain conductance were comparable to the National regional hydrological model, also covering the study area [35], where drain conductance in sandy areas was estimated to 0.04 1/d.

4.2. Streamflow Contribution of Flow Balance Components

The calibration of the two model conceptualizations yielded similar results in terms of calibration statistics with slightly better performance of the IRREG model (Figure 5). However, the resulting contribution to streamflow differed between the two conceptualizations where the drain contribution was substantially larger in the IRREG model (Table 2). The drainage contribution to streamflow was 31% in the IRREG model compared to 13% in the REG model. The reason for the large difference in the drain flow contribution is due to the numerical discretization along the streams. Due to the finer discretization in the IRREG model, the drainage area increased by almost 8 km2 as it was possible to specify a drain coverage closer to the stream cells (compare a cell area of 6.25 × 6.25 m for the IRREG model to a cell area of 100 × 100 m in the REG model). An example of the different drainage and stream areas is given in Figure 7. The additional drainage area in the IRREG model was equivalent to a 7% increase compared to the REG model but lead to a 146% increase in the drainage contribution to streamflow in the IRREG model. These results clearly showed the effect of refining the discretization along stream segments and accentuate the importance of how streams and drains in stream valleys are represented in numerical models. Substantial drainage contribution to the total streamflow has been shown in other studies. The drainage contribution to streamflow was 41% in the study by King et al. [14] or as high as 67% in the field investigations by Steiness et al. [16]. It is well known that riparian lowlands, although representing a small fraction of a catchment area, their ability to remove e.g., nitrate, is very high [36,37]. The most abrupt changes in topography and hydraulic gradients are often at the hill slope of stream valleys bordering riparian lowlands and their adjacent upland areas. Thus, it can be the ratio of catchment to lowland area that defines the split between seepage and drain. Therefore, with too large discretization, the likelihood of misrepresenting the stream and drainage areas is eminent in catchment scale models, which will affect the split between how much of the regional groundwater input reaches the stream by drains or seepage through the streambed. In conjunction, modeling of non-point nutrient dynamics with such models would be prone to misinterpretation given the potential erroneous flow balance contribution to streamflow.

4.3. GW-SW Exchange Fluxes

At the downstream part of Holtum stream, the groundwater fluxes to the stream estimated in the field were comparable with calculated groundwater fluxed from the two conceptualizations (Figure 6). Moving upstream from around 7500 m to the headwaters, the calculated fluxes—up to 3.5 m/d and 2.6 m/d for the REG and IRREG model—are somewhat higher than the measured fluxes (Figure 6) and the previously largest simulated fluxes of around 1.6 m/d [2]. As these relatively large calculated groundwater fluxes coincided with areas being flooded, it is concluded that the GW-SW exchange is not correctly represented in these areas within both models. This was supported by the interpreted hydraulic head contour map, Figure 1. Here, there were no strong gradients towards the stream that would indicate gaining reaches in the upstream part of Holtum stream. As such, the drainage contribution to the streamflow could be larger than estimated in these conceptualizations (Table 2). Thus, with a larger contribution from drainage in the upper part of Holtum stream, estimated groundwater fluxes could also be correspondingly lower and similar to measured groundwater flux ranges. Hansen et al. [15] showed that including drainage in a prior undrained area reduced groundwater fluxes substantially. Compared to the REG model, the calculated range and average exchange flux in the IRREG model yielded values in closer correspondence with the measured values. The lowest GW-SW exchange fluxes from the measurements, 0.06 m/d, were also captured in the IRREG model with values as low as 0.01 m/d, while the lowest calculated flux in the REG model was 0.3 m/d. The calculated groundwater fluxes in stream reaches with measurements showed a general pattern of being larger (Figure 6). However, it is noted that the measured groundwater fluxes do not represent steady-state but represent a snapshot in time of transient fluxes [2,4,24].
GW-SW exchange fluxes from especially the upstream area could likely be better represented if a more rigorous calibration was performed by using regularization such as pilot points [38] for Kh in the different geological units, dividing the stream reaches into several segments representing different Kv streambed values, and/or by varying the drainage conductance throughout the catchment. Although this was not part of the scope for the current study, the worth of incorporating measured GW-SW exchange fluxes is clear. With large differences in especially upstream calculations of GW-SW exchange fluxes, measured exchange fluxes could be used to constrain the models and thus aid in achieving representative flow balance contributions to streamflow. A common approach for assessing groundwater–stream exchange in catchment models is evaluating the simulations against observed hydraulic heads and stream flows [5,39]. Schilling et al. [5] argue that incorporating unconventional data such as GW-SW exchange fluxes in addition to hydraulic heads and streamflow in systematic model calibration could enhance the model predictability in terms of flow and mass balance.

5. Conclusions

The results of this study showed that the discretization used to specify the streams and drainage in their proximity had a substantial impact on simulated flow balance contribution to streamflow. As these differences in flow balance contribution consequently lead to differences in exchange fluxes that are crucial in simulating GW-SW exchange of flow and nutrients dynamics, it is paramount that future 3D numerical modeling can accurately estimate flow balance contribution to streamflow. A key challenge is to gather measurements of the drainage contribution on the catchment scale [18]. Here, we showed that in-stream measurements of GW-SW exchange fluxes are readily applicable for comparison with model calculated GW-SW exchange fluxes. Therefore, with the flexibility in calibration tools (e.g., PEST) to incorporate GW-SW exchange fluxes as calibration targets [5] and numerical codes allowing for flexible implementation of complex geology and discretization (e.g., MODFLOW-USG), more research is needed to investigate if and how GW-SW exchange fluxes can also be used to serve as a proxy to estimate the drainage contribution to streamflow on the catchment scale.

Author Contributions

S.K.; conceptualization, methodology, investigation, formal analysis, visualization, writing—original draft preparation, writing—review and editing. M.J.; investigation, formal analysis, writing—review and editing. J.K.; investigation, formal analysis, visualization, writing—review and editing. J.A.R.-G.; investigation, formal analysis, writing—review and editing. T.B.; investigation, formal analysis, writing—review and editing. P.E.; investigation, formal analysis, writing—review and editing. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Constantz, J. Heat as a tracer to determine streambed water exchanges. Water Resour. Res. 2008, 44. [Google Scholar] [CrossRef]
  2. Karan, S.; Engesgaard, P.; Rasmussen, J. Dynamic streambed fluxes during rainfall-runoff events. Water Resour. Res. 2014, 50, 2293–2311. [Google Scholar] [CrossRef]
  3. Poulsen, J.R.; Sebok, E.; Duque, C.; Tetzlaff, D.; Engesgaard, P.K. Detecting groundwater discharge dynamics from point-to-catchment scale in a lowland stream: Combining hydraulic and tracer methods. Hydrol. Earth Syst. Sci. 2015, 19, 1871–1886. [Google Scholar] [CrossRef] [Green Version]
  4. Sebok, E.; Duque, C.; Engesgaard, P.; Boegh, E. Spatial variability in streambed hydraulic conductivity of contrasting stream morphologies: Channel bend and straight channel. Hydrol. Process. 2015, 29, 458–472. [Google Scholar] [CrossRef]
  5. Schilling, O.S.; Cook, P.G.; Brunner, P. Beyond Classical Observations in Hydrogeology: The Advantages of Including Exchange Flux, Temperature, Tracer Concentration, Residence Time, and Soil Moisture Observations in Groundwater Model Calibration. Rev. Geophys. 2019, 57, 146–182. [Google Scholar] [CrossRef] [Green Version]
  6. Kidmose, J.; Engesgaard, P.; Nilsson, B.; Laier, T.; Looms, M.C. Spatial Distribution of Seepage at a Flow-Through Lake: Lake Hampen, Western Denmark. Vadose Zone J. 2011, 10, 110–124. [Google Scholar] [CrossRef]
  7. Karan, S.; Engesgaard, P.; Looms, M.C.; Laier, T.; Kazmierczak, J. Groundwater flow and mixing in a wetland–stream system: Field study and numerical modeling. J. Hydrol. 2013, 488, 73–83. [Google Scholar] [CrossRef]
  8. Rahimi, M.; Essaid, H.I.; Wilson, J.T. The role of dynamic surface water-groundwater exchange on streambed denitrification in a first-order, low-relief agricultural watershed. Water Resour. Res. 2015, 51, 9514–9538. [Google Scholar] [CrossRef]
  9. Munz, M.; Oswald, S.E.; Schmidt, C. Coupled Long-Term Simulation of Reach-Scale Water and Heat Fluxes Across the River-Groundwater Interface for Retrieving Hyporheic Residence Times and Temperature Dynamics. Water Resour. Res. 2017, 53, 8900–8924. [Google Scholar] [CrossRef]
  10. Williams, M.R.; King, K.W.; Fausey, N.R. Contribution of tile drains to basin discharge and nitrogen export in a headwater agricultural watershed. Agric. Water Manag. 2015, 158, 42–50. [Google Scholar] [CrossRef]
  11. Ford, W.I.; King, K.; Williams, M.R. Upland and in-stream controls on baseflow nutrient dynamics in tile-drained agroecosystem watersheds. J. Hydrol. 2018, 556, 800–812. [Google Scholar] [CrossRef]
  12. Goeller, B.C.; Febria, C.M.; Warburton, H.J.; Hogsden, K.L.; Collins, K.E.; Devlin, H.S.; Harding, J.S.; McIntosh, A.R. Springs drive downstream nitrate export from artificially-drained agricultural headwater catchments. Sci. Total Environ. 2019, 671, 119–128. [Google Scholar] [CrossRef] [PubMed]
  13. Rosenbom, A.E.; Karan, S.; Badawi, N.; Gudmundsson, L.; Hansen, C.H.; Nielsen, C.B.; Plauborg, F.; Olsen, P. The Danish Pesticide Leaching Assessment Programme: Monitoring Results, May 1999–June 2019; Geological Survey of Denmark and Greenland: Copenhagen, Denmark, 2020.
  14. King, K.W.; Fausey, N.R.; Williams, M.R. Effect of subsurface drainage on streamflow in an agricultural headwater watershed. J. Hydrol. 2014, 519, 438–445. [Google Scholar] [CrossRef]
  15. Hansen, A.L.; Jakobsen, R.; Refsgaard, J.C.; Højberg, A.L.; Iversen, B.V.; Kjærgaard, C. Groundwater dynamics and effect of tile drainage on water flow across the redox interface in a Danish Weichsel till area. Adv. Water Resour. 2019, 123, 23–39. [Google Scholar] [CrossRef]
  16. Steiness, M.; Jessen, S.; Spitilli, M.; van’t Veen, S.G.; Højberg, A.L.; Engesgaard, P. The role of management of stream–riparian zones on subsurface–surface flow components. Water 2019, 11, 1905. [Google Scholar] [CrossRef] [Green Version]
  17. Steiness, M.; Jessen, S.; van’t Veen, S.G.; Kofod, T.; Højberg, A.L.; Engesgaard, P. Nitrogen-Loads to Streams: Importance of Bypass Flow and Nitrate Removal Processes. J. Geophys. Res. Biogeosciences 2021, 126, e2020JG006111. [Google Scholar] [CrossRef]
  18. King, K.W.; Williams, M.R.; Macrae, M.L.; Fausey, N.R.; Frankenberger, J.; Smith, D.R.; Kleinman, P.J.A.; Brown, L.C. Phosphorus Transport in Agricultural Subsurface Drainage: A Review. J. Environ. Qual. 2015, 44, 467–485. [Google Scholar] [CrossRef] [Green Version]
  19. De Schepper, G.; Therrien, R.; Refsgaard, J.C.; Hansen, A.L. Simulating coupled surface and subsurface water flow in a tile-drained agricultural catchment. J. Hydrol. 2015, 521, 374–388. [Google Scholar] [CrossRef]
  20. Sonnenborg, T.O.; Scharling, P.B.; Hinsby, K.; Rasmussen, E.S.; Engesgaard, P. Aquifer Vulnerability Assessment Based on Sequence Stratigraphic and 39 Ar Transport Modeling. Ground Water 2015, 54, 214–230. [Google Scholar] [CrossRef]
  21. Vervloet, L.S.; Binning, P.J.; Børgesen, C.D.; Højberg, A.L. Delay in catchment nitrogen load to streams following restrictions on fertilizer application. Sci. Total Environ. 2018, 627, 1154–1166. [Google Scholar] [CrossRef]
  22. Ranalli, A.J.; Macalady, D.L. The importance of the riparian zone and in-stream processes in nitrate attenuation in undisturbed and agricultural watersheds—A review of the scientific literature. J. Hydrol. 2010, 389, 406–415. [Google Scholar] [CrossRef]
  23. Panday, S.; Langevin, C.D.; Niswonger, R.G.; Ibaraki, M.; Hughes, J.D. MODFLOW–USG version 1: An unstructured grid version of MODFLOW for simulating groundwater flow and tightly coupled processes using a control volume finite-difference formulation. Tech. Methods 2013. [Google Scholar] [CrossRef]
  24. Jensen, J.K.; Engesgaard, P. Nonuniform Groundwater Discharge across a Streambed: Heat as a Tracer. Vadose Zone J. 2011, 10, 98–109. [Google Scholar] [CrossRef]
  25. Karan, S.; Kidmose, J.; Engesgaard, P.; Nilsson, B.; Frandsen, M.; Ommen, D.; Flindt, M.R.; Andersen, F.Ø.; Pedersen, O. Role of a groundwater–lake interface in controlling seepage of water and nitrate. J. Hydrol. 2014, 517, 791–802. [Google Scholar] [CrossRef]
  26. Hajati, M.C.; Frandsen, M.; Pedersen, O.; Nilsson, B.; Duque, C.; Engesgaard, P. Flow reversals in groundwater–lake interactions: A natural tracer study using δ18O. Limnologica 2018, 68, 26–35. [Google Scholar] [CrossRef]
  27. Geological Survey of Denmark and Greenland (GEUS). National Well Database (Jupiter). Available online: https://eng.geus.dk/products-services-facilities/data-and-maps/national-well-database-jupiter (accessed on 7 June 2021).
  28. Rasmussen, E.S.; Dybkjær, K.; Piasecki, S. Lithostratigraphy of the Upper Oligocene—Miocene succession of Denmark. Geol. Surv. Den. Greenl. Bull. 2010, 22, 1–92. [Google Scholar] [CrossRef] [Green Version]
  29. Houmark-Nielsen, M. Pleistocene Glaciations in Denmark: A Closer Look at Chronology, Ice Dynamics and Landforms. Dev. Quat. Sci. 2011, 15, 47–58. [Google Scholar] [CrossRef]
  30. Niswonger, R.G.; Prudic, D.E. Documentation of the Streamflow-Routing (SFR2) Package to Include Unsaturated Flow Beneath Streams—A Modification to SFR1. Tech. Methods 2005. [Google Scholar] [CrossRef] [Green Version]
  31. Karan, S.; Sebok, E.; Engesgaard, P. Air/water/sediment temperature contrasts in small streams to identify groundwater seepage locations. Hydrol. Process. 2017, 31, 1258–1270. [Google Scholar] [CrossRef]
  32. Doherty, J. Calibration and Uncertainty Analysis for Complex Environmental Models; Watermark Numerical Computing: Brisbane, Australia, 2015. [Google Scholar]
  33. Henriksen, H.J.; Troldborg, L.; Sonnenborg, T.; Højberg, A.L.; Stisen, S.; Kidmose, J.B.; Refsgaard, J.C. Hydrologisk Geovejledning—God Praksis i Hydrologisk Modellering (Hydrological Geoguide—Best Practice in Hydrological Modeling). Guide 2017/1; De Nationale Geologiske Undersøgelser for Danmark og Grønland—GEUS (Geological Survey of Denmark and Greenland): Copenhagen, Denmark, 2017; p. 130.
  34. Sebok, E.; Müller, S. The effect of sediment thermal conductivity on vertical groundwater flux estimates. Hydrol. Earth Syst. Sci. 2019, 23, 3305–3317. [Google Scholar] [CrossRef] [Green Version]
  35. Stisen, S.; Ondracek, M.; Troldborg, L.; Schneider, R.J.M.; van Til, M.J. National Vandressource Model. Modelopstilling og kalibrering af DK-model 2019. In Danmarks Og Grønlands Geologiske Undersøgelse Rapport 2019/31; Geological Survey of Denmark and Greenland: Copenhagen, Denmark, 2019. [Google Scholar]
  36. Hill, A.R. Nitrate Removal in Stream Riparian Zones. J. Environ. Qual. 1996, 25, 743–755. [Google Scholar] [CrossRef]
  37. Hill, A.R. Groundwater nitrate removal in riparian buffer zones: A review of research progress in the past 20 years. Biogeochemistry 2019, 143, 347–369. [Google Scholar] [CrossRef]
  38. Anderson, M.P.; Woessner, W.W.; Hunt, R.J. Applied Groundwater Modeling: Simulation of Flow and Advective Transport; Academic Press: Cambridge, MA, USA, 2015. [Google Scholar]
  39. Hansen, A.L.; Refsgaard, J.C.; Christensen, B.S.B.; Jensen, K.H. Importance of including small-scale tile drain discharge in the calibration of a coupled groundwater-surface water catchment model. Water Resour. Res. 2013, 49, 585–603. [Google Scholar] [CrossRef] [Green Version]
Water 13 01923 g001 550
Figure 1. Study area. The insertion depicts the study area in red situated in the western part of Denmark.
Figure 1. Study area. The insertion depicts the study area in red situated in the western part of Denmark.
Water 13 01923 g001
Water 13 01923 g002 550
Figure 2. Data sets used in the setup of the geological model of the groundwater catchment. The cross-section through the geological model is shown in Figure 3.
Figure 2. Data sets used in the setup of the geological model of the groundwater catchment. The cross-section through the geological model is shown in Figure 3.
Water 13 01923 g002
Water 13 01923 g003 550
Figure 3. Geological cross-section through the study area (for location see Figure 2). There is a buried valley cutting through Miocene deposits in the NE area of the groundwater catchment. Top elevations of the Miocene formations are at the greater depths east of the main stationary line.
Figure 3. Geological cross-section through the study area (for location see Figure 2). There is a buried valley cutting through Miocene deposits in the NE area of the groundwater catchment. Top elevations of the Miocene formations are at the greater depths east of the main stationary line.
Water 13 01923 g003
Water 13 01923 g004 550
Figure 4. Relative sensitivities in percent.
Figure 4. Relative sensitivities in percent.
Water 13 01923 g004
Water 13 01923 g005 550
Figure 5. Simulated versus observed hydraulic heads in the two model conceptualizations. The gray line corresponds to the 1:1 line.
Figure 5. Simulated versus observed hydraulic heads in the two model conceptualizations. The gray line corresponds to the 1:1 line.
Water 13 01923 g005
Water 13 01923 g006 550
Figure 6. Calculated groundwater discharge along Holtum stream from the two model conceptualizations and previously estimated values [4,7,16,31,34].
Figure 6. Calculated groundwater discharge along Holtum stream from the two model conceptualizations and previously estimated values [4,7,16,31,34].
Water 13 01923 g006
Water 13 01923 g007 550
Figure 7. Drainage area in the two conceptualizations for a part of Holtum stream. Blank cells represent the area occupied by stream cells in the IRREG model. The largest cell is 100 × 100 m. The dark green area delineates the drain polygon of the REG model. The light green area delineates the additional drainage area in the IRREG model and the area occupied by stream cells in the REG model.
Figure 7. Drainage area in the two conceptualizations for a part of Holtum stream. Blank cells represent the area occupied by stream cells in the IRREG model. The largest cell is 100 × 100 m. The dark green area delineates the drain polygon of the REG model. The light green area delineates the additional drainage area in the IRREG model and the area occupied by stream cells in the REG model.
Water 13 01923 g007
Table
Table 1. Calibrated parameter values and their associated 95% confidence intervals.
Table 1. Calibrated parameter values and their associated 95% confidence intervals.
Model ConceptualizationParameterOptimized ValueLower 95%
Confidence		Upper 95%
Confidence		Units
Horizontally
regular grid (REG)		Kh Quaternary sand321471(m/d)
Kh Miocene sand161419(m/d)
Kh Miocene clay0.0650.0320.13(m/d)
Kv streambed0.730.391.3(m/d)
Conductance drain0.0780.0001932(1/d)
Conductance general head0.0320.00860.12(1/d)
Horizontally
irregular grid (IRREG)		Kh Quaternary sand301464(m/d)
Kh Miocene sand171420(m/d)
Kh Miocene clay0.0670.0320.14(m/d)
Kv streambed0.590.311.2(m/d)
Conductance drain *0.05--(1/d)
Conductance general head0.0140.00610.030(1/d)
* Fixed value.
Table
Table 2. Flow balance distribution relative to the total outflow in percent.
Table 2. Flow balance distribution relative to the total outflow in percent.
Flow BudgetHorizontally Regular
Grid (REG)		Horizontally Irregular
Grid (IRREG)
Stream flow80.563.0
Drainage12.530.9
General head boundary7.06.1
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Karan, S.; Jacobsen, M.; Kazmierczak, J.; Reyna-Gutiérrez, J.A.; Breum, T.; Engesgaard, P. Numerical Representation of Groundwater-Surface Water Exchange and the Effect on Streamflow Contribution Estimates. Water 2021, 13, 1923. https://doi.org/10.3390/w13141923

AMA Style

Karan S, Jacobsen M, Kazmierczak J, Reyna-Gutiérrez JA, Breum T, Engesgaard P. Numerical Representation of Groundwater-Surface Water Exchange and the Effect on Streamflow Contribution Estimates. Water. 2021; 13(14):1923. https://doi.org/10.3390/w13141923

Chicago/Turabian Style

Karan, Sachin, Martin Jacobsen, Jolanta Kazmierczak, José A. Reyna-Gutiérrez, Thomas Breum, and Peter Engesgaard. 2021. "Numerical Representation of Groundwater-Surface Water Exchange and the Effect on Streamflow Contribution Estimates" Water 13, no. 14: 1923. https://doi.org/10.3390/w13141923

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop