The Role of Management of Stream–Riparian Zones on Subsurface–Surface Flow Components

Abstract

A managed riparian lowland in a glacial landscape (Holtum catchment, Denmark) was studied to quantify the relative importance of subsurface and surface flow to the recipient stream. The hydrogeological characterization combined geoelectrical methods, lithological logs, and piezometric heads with monthly flow measurements of springs, a ditch, and a drain, to determine seasonality and thereby infer flow paths. In addition, groundwater discharge through the streambed was estimated using temperature and water-stable isotopes as tracers. The lowland received large groundwater inputs with minimal seasonal variations from adjacent upland aquifers. This resulted in significant amounts of groundwater-fed surface flow to the stream, via man-made preferential flow paths comprising ditches, drainage systems, and a pond, and via two natural springs. Roughly, two thirds of the stream gain was due to surface flow to the stream, mainly via anthropogenic alterations. In contrast, direct groundwater discharge through the streambed accounted for only 4% of the stream flow gain, although bank seepage (not measured) to the straightened and deepened stream potentially accounted for an additional 17%. Comparison to analogous natural flow systems in the catchment substantiate the impact of anthropogenic alterations of riparian lowlands for the subsurface and surface flow components to their streams.

1. Introduction

Non-point source contamination of streams with nutrients, applied for agricultural production, stresses natural environmental systems. In Denmark, despite a nearly 50% reduction in nitrate leaching from the root zone [1], additional reductions are required to meet the goal of the EU Water Framework Directive. Riparian lowlands may serve as natural buffers, which may attenuate the movement of nutrients to streams [2]. Management of riparian lowlands are therefore often recommended as an integral part of watershed management practices in agricultural areas to reduce nutrient loading to streams [3,4,5]. Riparian lowlands used for agricultural purposes are frequently drained, but how much such drainage increases or decreases nutrient retention is uncertain. To effectively use riparian lowlands as buffer zones, or to implement new mitigation measures in the riparian lowlands, it is imperative to first understand the flow paths. This includes how discharge of groundwater from regional aquifers is distributed among e.g., direct groundwater discharge through the streambed, through drains to the stream, or from groundwater-fed seeps (focused and/or diffusive) leading to overland flow, which is then routed to the stream.

Hill [6] proposed different conceptual flow distribution models, where the hydrological controls were mainly determined by the size of the connected regional aquifer system and the thickness of the riparian aquifer. Several studies have since then, qualitatively confirmed these models. For example, studies [7,8] have demonstrated that flow across the lowland occurred as subsurface flow and groundwater-fed surface flow. Furthermore, Shabaga and Hill [8] demonstrated by dye tracer studies that surface flow took place as a combination of water discharging through macro-pores (soil pipes) and diffusively. Preferential flow through soil pipes in peat deposits has been shown to transmit large fluxes of water and with the possibility of causing springs resulting in overland flow [9].

Measuring the individual flow paths is not a trivial task and thus not often done. Brüsch and Nilsson [10] attempted to measure total surface runoff (i.e., the sum of all overland flow) from a riparian zone and found an average surface runoff of 0.3 L s⁻¹, compared with stream flow of 24–38 L s⁻¹. Furthermore, diffuse discharge to the surface only accounted for 3% of the runoff, so most of this occurred through focused seeps. Shabaga and Hill [8] measured surface flow directly to a river at one site. Johansen et al. [11] measured flow in three springs and a network of ditches. In none of these studies, the direct seepage through the streambed was measured. It was therefore not possible to quantify the relative role of subsurface flow to overland flow. Langhoff et al. [12] attempted to measure all flow paths. Direct groundwater discharge to the streams was measured by use of seepage meters at eight sites. Bank seepage and overland flow was estimated or measured as well. At three of the sites, more than 50% of the discharge was composed of direct groundwater discharge through the streambed; at one site, surface flow contributed all the stream gain. From their study, it is evident that surface flow is a non-negligible flow path with a potential to represent a significant part of the overall stream–gain–water balance. Frederiksen et al. [13] performed a similar, but less intense, study to evaluate methods for estimating groundwater discharge at different scales. Here they compared the relative contribution of groundwater discharge measured by seepage meters, tile drainage, and overland flow to that of stream flow gain from differential gauging and hydrograph separation. Their study catchment could be divided into two separate units; an upstream reach, which was tile drainage dominated and a lower reach, which was groundwater dominated. In the latter, groundwater discharge was >60% of stream flow gain.

The objectives and overall methodology of this study were to: (1) investigate the connection between the groundwater system and the recipient stream in a managed riparian lowland, based on frequent monitoring of major water flux components of the stream–gain–water balance throughout a year, (2) quantify the relative importance of subsurface and surface flow to the stream, and (3) discuss these results in relation to similar studies in the same catchment on of how management (or mismanagement) of a riparian lowland may change flow paths. The hydrogeology of the riparian lowland and adjacent upland was characterized extensively using geoelectrical methods, hand-drillings, water-stable isotopes, piezometric head observations, along with monthly monitoring of flow in several springs, a ditch, and through the streambed.

2. Study Site

As study site, a riparian lowland located in Holtum catchment (central Jutland, western Denmark) was chosen (Figure 1). From north to south through Holtum catchment runs the Main Stationary Line (MSL) of the late Weichselian glaciation [14]. West of the MSL outwash plain sediments dominate, while clayey tills are frequent east of the MSL. Land use in the catchment is mainly agriculture (49%), followed by forest (31%), urban (15%), and wetland and surface-water areas (5%). The annual precipitation in the catchment is 984 ± 55 mm year⁻¹ (mean ± 1 standard deviation), and actual evapotranspiration is 510 ± 45 mm year⁻¹ [15].

The riparian lowland is located along a section of Holtum stream running in the valley of an old sub-glacial stream trench, and hence the riparian aquifer comprises sandy outwash deposits. Hillslopes north and south of the lowland area comprise clayey tills. Beneath these Quaternary deposits sits a regional Miocene micaceous sand aquifer. The stream is a lowland second-order gaining stream with a width of approx. three meters at the field site. The stream is perennial and has been straightened and likely deepened at the field site, a very frequent management practice along streams in agricultural landscapes. The height of the stream bank is about 1–1.3 m. The catchment boundary lies about 1 km southwest and 0.4 km northeast of the study site (distances perpendicular to the stream). Land use southwest of the riparian lowland is agriculture all the way to the catchment boundary. A forest plantation, mainly with pine trees, is located to the northeast of the riparian lowland and here extends to the catchment boundary.

Figure 1A,B displays the observation points, land use, and drainage features, and defines the naming used henceforth in this paper. South of the stream, the riparian lowland consists of an 85 m wide wetland functioning as a cow pasture. The wetland is bordered by a steep hillslope to the south rising 12 m. The hillslope is cut by a small gorge in the center. Several old drainage wells and two springs (S1 and S2, Figure 1) are present in the wetland. One of the drainage wells (DW) is constantly overflowing. Springs and the drainage well contribute to diffusive overland flow. The wetland’s drainage system has not been well maintained and it is likely that drains connecting to the drainage wells are not properly functioning. A submerged drainpipe in the stream in the eastern part of the wetland exists, also indicated in Figure 1A,B. An elevated small pond is present in the middle of the wetland. The banks of the pond are to the north (towards the stream) and northeast constructed from earthen levees. The pond is drained through an overflow pipe (Pond Outlet, PO) also contributing to overland flow. As a result, a small rivulet (RI) runs from the wetland into the stream. To the west, a nearly overgrown ditch (DI) is deep enough to reach sand beneath the foot of the clayey till hillslope; this forms a spring at the southern end of DI. The water captured by DI discharges directly into the stream. Several similar ditches exist downstream of the field site (between Q2 and Q3, Figure 1C), with discharge mainly consisting of groundwater. S1, S2, DW, PO, RI and DI are collectively termed “springs” from here on as they represent a visible and measurable flow. Furthermore, all “springs” one way or another represent water emerging to the surface in a, at least apparently, spring-like fashion. Along the foot of the hillslope, seeps have been observed maintaining wet surface areas. The wetland area is therefore mostly very wet, with areas fully saturated throughout the year, and consists of the main part of saturated peat and mud.

In 2015, a transect consisting of ten piezometers with screens at varying depths was established [16]. Some of these piezometers are shown in Figure 1A (along the ditch; S1, C5, C10, C20, C30).

On the north side of the stream, land use consisted of an approx. 35 m wide buffer strip of grass functioning as a cow pasture and a 145 m wide field for crop production extending from the buffer strip to the northern valley hillslope. Field crops have alternated between corn and oats. A field drain (DR), leading directly into the stream (not submerged and therefore easily measurable), drains a part of the agricultural area which is prone to waterlogging due to a depression in the surface topography.

3. Materials and Methods

3.1. Hydrogeological Characterization

Electrical Resistivity Tomography (ERT) surveys were done on both sides of the stream similar to surveys carried out in comparable settings [17,18]. A SYSCAL PRO resistivity meter with a Wenner configuration was used. The ERT1 and ERT2 profiles were measured in November 2015; the ERT3 profile was measured in March 2016 and the ERT4 profile in March 2017. All profiles (locations shown in Figure 1A) were measured using 96 electrodes with 2 m electrode spacing; equaling a length of each profile of 192 m and a maximum penetration depth of ~35 m in the center of the profiles. The location and elevation (meters above sea level; MASL) of each electrode were measured with a Differential Global Positioning System (DGPS) using a Trimble^® R8s GNSS receiver with a Trimble TSC3 controller (Trimble Inc., Sunnyvale, CA, USA). The inversion of the measured data was done with Res2DInv by Geotomo software (Res2DInv ver. 3.57; Geotomo Software, Penang, Malaysia), using a least squared inversion technique to produce the tomographic profiles.

Common offset Ground-Penetrating Radar (GPR) surveys were conducted across four transects in the wetland south of the stream (Figure 1B). The GPR survey was carried out to identify the thickness and spatial geometry of the peat in the wetland. GPR has been frequently used to identify the peat-layer thickness in wetlands [19,20]. A PulseEKKO^® PRO SmartTow GPR system (Sensors & Softwarer Inc., Mississauga, Canada) was used with a 250 MHz transducer, with a separation of 0.4 m between the transmitter and receiver. A step size of 0.05 m was used and a stacking of 256 pulses, collected at each position along the GPR transects, with the receiver set with a sample time window of 400 ns. No topographical corrections were made since the ground surface is relatively flat in the surveyed area. The profile length was measured using an odometer. The collected data was post-processed in Ekko_View Deluxe (Sensors & Software Inc., Mississauga, Canada). A high pass (DEWOW) filter was applied to reduce low frequency noise and The Spherical Exponential Calibrated Compensation (SEC2) gain was used to compensate for spherical spreading loses and dissipation of the radar signal.

Eighteen shallow boreholes were hand-drilled, using an Eijkelkamp auger (Eijkelkamp Soil & Water, Geisbeek, the Netherlands), in the wetland to determine the peat-layer thickness and to ground truth depth of reflectors and obtain lithological logs from sediments retracted from the boreholes.

3.2. Hydrology

In addition to the existing piezometers, 14 new piezometers were installed in a northeast southwest transect perpendicular to the stream. Galvanized steel pipes (Ø2.54 cm) equipped with a 9 cm long screen, were hammered into the subsurface using a pneumatic hammer. Universal Transverse Mercator (UTM) coordinates of the piezometers including elevation of top of well and ground surface were measured with a DGPS.

Naming of piezometers, presented in Figure 1, were done in chronological order of installation. Installation of piezometers TH1d, TH2d, TH3, and TH4 took place in November 2015. Piezometers TH3 and TH4 had to be removed in March 2016 to allow the farmer to access the field. In September 2016, TH3 was re-established as TH3.2 and removed again in March 2017.

The transect was extended in March 2016 with the addition of TH1s and TH2s, next to the existing TH1d and TH2d. On the same occasion, TH6 was also installed approx. 1.5 m from the stream. On the opposite side of the stream, in the wetland, Installation of TH8, approx. 1 m from the stream, took place in January 2017. Two piezometers TH7d and TH7s, forming a nest, were installed in October 2016 and January 2017 respectively. Installation of TH10 near the hillslope, and another nest (consisting of TH9s and TH9d) took place in March 2017. Hence, four piezometer nests were installed at the field site, two on each side of the stream; nested piezometers are referred to as TH1, TH2, TH7, and TH9.

Nine of the piezometers were slug-tested using the falling head method, five of which was slug-tested at multiple screen depths. The slug-tests were conducted three consecutive times and analyzed using the Hvorslev method for unconfined aquifers.

Daily values of precipitation for 2017 were obtained from Danish Meteorological Institute’s Nørre Snede station (UTM-E: 524034, UTM-N: 6202165) approx. 4 km east of the field site.

3.3. Discharge Measurements in Stream and Springs

Stream discharge and flows from the springs were measured once every month and the monthly measurement was repeated twice. It would have been better if we had continuous stream discharge data, but, on the other hand, it was not possible for us to measure continuous flow from the springs. Therefore, the two types of measurements go in tandem and provide a snapshot of the monthly stream flow gain and spring flows. This also means that our water balance is the average of these monthly stream flow gains and spring flows.

Stream discharge was measured using an OTT Acoustic Digital Current meter (ADC) (OTT Hydromet, Kempten, Germany). The ADC uses a 6 MHz acoustic sensor for point flow velocity measurements, and corrects the measurement for water depth (water above the sensor) and temperature. Stream discharge was measured using the mean-section method at three locations; immediately upstream of the field site (Q1), right after the ditch draining the wetland (Q2) and downstream of the field site (Q3) (Figure 1C). The distances along the stream are 250 m and 750 m between Q1 and Q2, and Q2 and Q3, respectively. Measurements were performed monthly (except in July 2017) from January 2017 to January 2018.

“Spring” discharge in the wetland (Figure 1A,B) was measured using a cutthroat flume (Baski Inc.; Denver, CO USA). Using a spirit level to keep the base of the flume horizontal in both the longitudinal (direction of flow) and transverse directions. The banks on either side of the flume were dammed to make sure water would not flow around the flume. After installation, the flume was left for up to half an hour before measurements were taken, to allow flow to equilibrate and to ensure free flow conditions had been reached. Experimental work [21] have established a relationship between the head $\left(h_a\right)$, upstream of the flume throat, and discharge rate $\left(Q\right)$, given different sizes of cutthroat flumes under free flow conditions:

$$Q_{freeflow}=C_f{h}_a^{n_1}$$

(1)

where $Q_{freeflow}$ is the flow rate in $L^3{t}^{-1}$; $C_f$ is the free flow coefficient, which is a function of flume length, given by a flume length coefficient, and throat width of the flume $\left(C_f=K_l·W^{1.025}\right)$; and $n_1$ is the free flow exponent.

3.4. Direct Seepage through the Streambed: 1D Vertical Temperature Profiling

The mean near-surface groundwater temperature for Danish conditions and within the catchment is relatively constant at ~8 °C [17,22,23]. This feature allows the use of natural temperature variations as tracer to determine the flux between groundwater and stream. Vertical groundwater fluxes were hence calculated from measurements of Vertical Temperature Profiles (VTPs) and thermal conductivity of the streambed.

The measurements were carried out in September 2016, where the contrast between groundwater and stream water was relatively high after the summer period and flow and heat transport conditions were assumed to be at steady state. The same method has previously been used to estimate stream- and groundwater exchange in Holtum stream [16,17,22]. Notice that we here only have one estimate (September) of the streambed seepage. We quickly found that streambed seepage was much smaller than flows from the springs, so we did not repeat this monthly. This is quite a time-consuming process.

Three VTPs, one near each of the banks and one in the center of the stream, were measured for 11 transects along a 90 m stream section (Figure 1B). Each VTP was recorded after a 10-min equilibration time [24] at 0, 0.025, 0.05, 0.075, 0.1, 0.15, 0.2, 0.3, 0.4, and 0.5 m depth below the streambed surface. A steady-state analytical solution [25] for a 1-dimensional conduction–convection system was used to estimate the vertical groundwater fluxes.

$$T\left(z\right)=T_s+\left(T_g-T_s\right)\frac{exp\left(\frac{N_{p{e}_z}}{L}-1\right)}{exp\left(N_{pe}-1\right)}$$

(2)

where $T\left(z\right)$ is the temperature (°C) of the streambed sediment at depth $z\left(m\right)$; $T_s$ is the temperature (°C) of the stream water at the streambed (measured at depth 0); $T_g$ is the groundwater temperature (°C) at depth L $\left(m\right)$ and $N_{pe}$ is the Peclet number expressing the ratio caused by convection and conduction. For this, study $T_g$ and L are defined as 8.3 °C and 4 m.

$$N_{pe}=\frac{q_z·{\rho }_f·C_f·L}{{\mathsf{\kappa }}_e}$$

(3)

where $q_z$ is the vertical Darcy flux $\left(\mathrm{m}{\mathrm{s}}^{-1}\right)$; ${\rho }_f·C_f$ is the volumetric heat capacity of the water $\left(\mathrm{J}{\mathrm{m}}^{-3}{}^{°}{\mathrm{C}}^{-1}\right)$ and ${\mathsf{\kappa }}_e$ is the thermal conductivity of the sediment $\left(\mathrm{J}{\mathrm{m}}^{-1}{\mathrm{s}}^{-1}{}^{°}{\mathrm{C}}^{-1}\right)$. The accuracy of the temperature measurements have been reported to be ±0.2 °C [16,24].

The thermal conductivity (${\mathsf{\kappa }}_e$) of the top 3–10 cm of the streambed sediment was measured in the field using a KD2 pro with the SH-1 sensor (Decagon Devices, Pullman, WA, USA). The instrument has an accuracy of ±10% at thermal conductivities between 0.002 and 2 $\mathrm{W}{\mathrm{m}}^{-1}{}^{°}{\mathrm{C}}^{-1}$. The sensor is equipped with a heater needle spaced 6 mm apart from a temperature sensor needle. A read time of 2 min was used, during which a heat pulse was generated for the first half of the read time and then data (temperature) was recorded to compute thermal properties for the full read time. Prior to initiation of the heat pulse, the instrument allows 30 s for temperature equilibration of the sensor. Measurements were conducted next to each VTP measurement, yielding the bulk thermal conductivity and a measurement error for each location. Previous studies have used the same instrument to measure in situ ${\mathsf{\kappa }}_e$ [24] and ${\mathsf{\kappa }}_e$ of cores collected in Holtum catchment [26]. Both studies found that the thermal conductivity can assume a wide range of values. In general, the thermal conductivity of a porous material depends on the sediment mineral composition, grain size, pore space, and water content hereof (in this case fully saturated). Thermal properties of selected materials are available in the literature [27].

3.5. Water Sampling and Analysis

Water samples were collected for water-stable isotopes (²H and ¹⁸O) and water quality, an- and cations, alkalinity, Fe²⁺ and H₂S. In this paper, only the water-stable isotopes will be reported. Samples from springs in the wetland (see Figure 1A,B) were collected monthly between January 2017 and January 2018 (excluding July 2017). Grab samples of stream water and the drain on the crop side were collected for the same period.

Groundwater samples from the piezometers were collected on several occasion. During installation of piezometers screened deeper than 2 m below ground level (MBGL) (TH1d, TH1s, TH2d, TH2s, TH7d, TH7d and TH9d) multilevel water sampling was conducted by taking samples at several depths. The piezometers were further sampled 2–5 times between January 2015 and December 2017.

Groundwater samples were retrieved using a peristaltic pump. Following three purges, water was pumped through a flow cell equipped with electrodes for dissolved oxygen (DO; WTW FDO 925 IDS electrode), pH (Hach PHC101 electrode), Electrical Conductivity (EC) and temperature (both recorded using an Hach CDC401 electrode) connected to WTW Oxi 3310 IDS and HACH HQ30 flexi instruments. Sampling commenced when values were stable. Samples were taken before water entered the flow cell, using a syringe and a three-way valve.

All samples of stable hydrogen and oxygen water isotopes were filtered with a 0.20 µm Minisart cellulose acetate syringe filter by Sartorius (Göttingen, Germany), into a 1.5 mL glass vial in the field and then stored in a refrigerator. The samples were analyzed at the Geological Survey of Denmark and Greenland (GEUS) on a PICARRO L2120-i (PICARO Inc., Santa Clara, CA, USA) using cavity ring-down spectroscopy (CRDS). All values for isotopes are expressed in per-mill (‰) with the δ-notation indicating the deviation from the VSMOW (Vienna Standard Mean Ocean Water) standard:

$$\mathsf{\delta }\text{-}\mathrm{value}(‰)=\left(\frac{R_{Sample}-R_{VSMOW}}{R_{VSMOW}}\right)\times 1000‰$$

(4)

where R is the ratio of the heavy to the light isotope.

In June 2017, a campaign was undertaken to investigate the spatial patterns of the water-stable isotopes of groundwater directly underneath the peat cover of the wetland. Nineteen point samples were obtained using stainless steel mini-piezometers with an inner diameter of 0.5 cm and a 5 cm screen length. The mini-piezometer was pushed vertically down to ~1 m below the surface by hand or with the use of a hammer, and withdrawn again after sampling. Sampling followed the same procedure as described above. Likewise, a mini-piezometer was later (November 2018) used to sample water 1 m beneath the streambed.

4. Results

4.1. Hydrogeological Characterization

Figure 2 shows the geology as interpreted by comparing resistivity cross sections to four borehole logs available from the Danish National Well Database (Jupiter), in which boreholes are denoted with a unique borehole number (DGU no.). The logs are from DGU Nos. 96.2187, 96.1367 and 96.1026, the locations of which are shown in Figure 1A. DGU No. 96.1076 is located about 425 m north of the end of ERT1.

A zone near the surface with resistivity values above 600 Ωm is visible to the southwest (from 0 to about 38 m along the profile). This is most likely due to the steep increase in slope, which facilitates the presence of a deeper unsaturated zone. Just below the unsaturated zone is a ~4 m thick clay lens with low resistivity values (30–70 Ωm), which could be part of a more extensive clay till layer seen in borehole 96.2187 from 67.5 to 64 MASL. A saturated gravel or sand layer (100–400 Ωm) was found underneath, again in accordance with 96.2187 with gravelly outwash deposits from 64 to 54.5 MASL. The upland aquifer depth to the south (southwest), considering only the unconfined top unit of permeable sediments, is 15.5 m thick.

The geology described in borehole 96.1367 consists of a peat layer of 1 m thickness followed by 14 m sand. Data from ERT4 is in reasonable agreement with the lithology of the borehole. Resistivities in the top of the profile are below 30 Ωm representing peat. Below, resistivities are in the range of 70–400 Ωm corresponding to saturated gravelly or coarse sand. The peat has been formed in a small depression in the top-center of the profile. The peat layer has a length of about 40 m and a thickness ranging from less than 0.5 m up to approx. 3 m. A gravelly sand layer with varying thickness of ~5–8 m is found below the peat followed by a sand layer (resistivity values ~70 Ωm) extending to the bottom of the profile. At the southern edges and deepest part of profile ERT4, a clay layer is likely present, since resistivity values are <50 Ωm. Borehole 96.2187 shows a micaceous clay layer at a depth of 54.5–26 MASL and it may be that this clay layer has a trough-like bathymetry extending beneath the whole stream valley (but not captured with ERT).

On the crop side of the stream (ERT1), the sandy aquifer is unconfined with resistivity values ranging from 70 to 400 Ωm. The aquifer is covered by approx. 0.5–1 m sandy soil with high resistivity values (>400 Ωm), where the soil is unsaturated. Smaller areas near the surface along the profile have resistivity values ranging from ~90–190 Ωm indicating saturated soils all the way to the surface. These are areas prone to water logging and are drained (DR).

At the far northern end of profile ERT1, resistivity values are high associated with increase in thickness of the unsaturated zone—or indicative of high resistivity sediments such as gravel. The geology of borehole 96.1026 consists of sand and gravel down to 60.3 MASL, then from 60.3 to 53 MASL a layer of sandy clay is present and below this is a clay layer (53–39 MASL). This is also indicated by the resistivity values found at the foot of the hillslope for ERT1, with sand and gravelly deposits on top of a moraine till/clay. The aquifer depth to the north (northeast), down to the micaceous clay, is about 13 m thick. This is comparable to the thickness found to the south of the stream.

The hydraulic conductivities (K) range from ~0.1 to 50 m day⁻¹ with a mean of 14 m day⁻¹. These K-values (Figure 2) are all representative for sandy aquifer materials and are similar to those found for a sandy riparian aquifer in the same catchment [17]. High K-values are found at depths of 2–6 MBGL for TH7, where also the ERT4 profile indicated the presence of gravelly or coarse sand. A similar agreement is with TH2 and TH1, both with higher hydraulic conductivities from 2–6 MBGL and 5–12 MBGL, respectively, in areas with coarse sands as indicated by ERT1. In addition, 11 out of the 18 hand-drilled boreholes were described as having coarse sand and some with pieces of gravel-sized particles.

A joint interpretation of ERT4, the GPR survey, and shallow borehole logs is shown in Figure 3, which displays the resulting thickness of the peat that overlies the sandy aquifer in the wetland. Peat thicknesses varies between 0.1 and 1.5 m. Springs are in elongate bands with peat thickness >0.5 m which runs parallel to the stream.

4.2. Hydrology

Figure 4 displays the seasonal fluctuations in precipitation, groundwater hydraulic heads, and stream stage. The stream stage and hydraulic heads followed the fluctuations of precipitation (Figure 4A,B). June had 114 mm of rainfall, which resulted in increased hydraulic heads in all piezometers within the same month. The same was seen for October with 136 mm of rainfall.

In all piezometers and throughout the monitoring period, hydraulic heads were higher than the stream stage (Figure 4). The corresponding average gradients between the most northern (F100) and southernmost (TH10) piezometers (Figure 4C) and the stream were 0.01 and 0.02, respectively (positive towards to the stream).

The two nested piezometers TH1d and TH1s, on the crop side of the stream and screened at 16 and 6 MBGL, respectively, showed a significant downward gradient. On average, the head observed in TH1s was 0.7 m higher than that of TH1d, and TH1s was the only piezometer with an artesian head at the field site. The nested piezometers TH2s and TH2d, located closer to the stream and screened at 5.5 and 12 MBGL, respectively, showed vertical head difference from 0.05 to 0.07 m over the monitoring period, yielding a persistent upward hydraulic gradient of about 0.01.

To the other side of the stream, the TH7d screened at nine MBGL always had higher hydraulic heads than the shallower TH7s piezometer (4 MBGL), varying from 0.03 to 0.08 m giving upward hydraulic gradients of between 0.006 and 0.016. TH9d and TH9s that are also located in the wetland, but further away from the stream, screened at 5.5 and 1.5 MBGL, respectively. Here, the hydraulic heads were generally higher for the shallow TH9s, yielding downward gradients ranging from 0.001 to 0.004.

The difference between maximum and minimum heads ranges from 0.16 to 0.37 m with an average of 0.21 m. The largest fluctuations were observed in F100 with −0.21 to 0.16 m of fluctuations from the observed average value of 59.91 MASL. Piezometer TH10 also had slightly higher fluctuations, when compared to piezometers further away from the edges of the riparian lowland, with −0.06 to 0.14 m of fluctuations from the average head of 59.39 MASL.

In general, the temporal fluctuations throughout the monitoring period on both sides of the stream were identical (Figure 4). Although the magnitude of the fluctuations in observed heads were slightly larger in the piezometers on the crop side, ranging from 0.16 m to 0.24 m, not including F100 furthest from the stream. Fluctuations were more dampened in the observed hydraulic head in the wetland ranging from 0.16 to 0.18 m, not including TH10 at the hillslope.

4.3. Stream and Spring Discharge

4.3.1. Stream Discharge

Measured stream discharges (Q_s) at the three locations are shown in Figure 5. Mainly for August, October (in 2017), and January (2018), the higher discharges can be attributed to higher intensity of precipitation (Figure 4A)—stream discharges were highest in August for all three locations. The average Q_s, excluding measurements affected by precipitation, for Q1, Q2, and Q3 were 0.239 m³ s⁻¹, 0.257 m³ s⁻¹, and 0.421 m³ s⁻¹, respectively.

The differences between Q1 and Q2 monthly are barely noticeable in Figure 5. The stream section between the two locations is, however, gaining for six out of the 11 months. There was no difference in discharge in February. Discharge measurement associated with rainfall event, was also when the observed largest gaining (April and October 2017) and losing conditions (August 2017 and Jan 2018) occurred in the stream between Q1 and Q2. The significant decrease in discharge between Q1 and Q2 (losing conditions, −0.22 m³ s⁻¹) seen in August, is certainly because Q2 was measured last and after a rainfall had stopped—whereas Q1 and Q3 was measured during rainfall. Excluding the month where discharge can be coupled with rainfall, the average gain between Q1 and Q2 yields 5.1 × 10⁻⁵ m³ s⁻¹ m⁻¹ with a standard deviation of 7.6 × 10⁻⁵ m³ s⁻¹ m⁻¹. Hence, it is not strictly possible to determine whether the stream section between Q1 and Q2 is gaining or loosing. For later comparison, stream gain per meter between Q1 and Q2 would correspond to a total gain of 4.58 × 10⁻³ m³ s⁻¹ along a 90 m stream reach, corresponding to where the VTP measurements was conducted.

Between Q1 and Q3, the stream was gaining and had an average gain of 1.9 × 10⁻⁴ m³ s⁻¹ m⁻¹ (standard deviation (SD) 5.4 × 10⁻⁵ m³ s⁻¹ m⁻¹) which is considerable higher than between Q1 and Q2. Seasonal variations are more visible at Q3 with lower discharge during late spring to late summer (May through September) and with higher discharge in the winter.

4.3.2. Wetland Springs Discharge

Measured discharges at the six springs are shown in Figure 6. Discharge ranged from 6.1 × 10⁻⁵ m³ s⁻¹ at S2 in December 2017 to 1.2 × 10⁻³ m³ s⁻¹ at S1 in February 2017. S1 had the highest average discharge of 9.4 × 10⁻⁴ m³ s⁻¹ followed by the rivulet (RI) with an average discharge of 6.7 × 10⁻⁴ m³ s⁻¹.

The seasonality of the discharge at RI, S1 and PO was especially similar (Figure 6). For all three locations, discharge was highest in late winter to early spring (February–April). Discharge decreased until September, with the lowest discharges. In early fall discharge increased again, until the last measurements in January 2018. The similarity in the seasonality of PO and RI is expected since the PO drains the pond into the bulrush-covered central part of the wetland. Immediately on the opposite side of the bulrush (towards the stream), the rivulet RI is formed and discharges into the stream. RI therefore likely receives water from the pond. However, discharge measured at RI is larger than PO (Figure 6), indicating that RI receives water not only from PO. This may include water from other springs upstream of RI, being S1, DW, and/or S2 or from diffusive groundwater upwelling in the central part of the wetland.

The ditch (DI) and drainage well (DW) display a sudden increase in discharge in August and continuing into the fall. Possibly, part of the water from DW runs as surface flow to DI and thus added to the groundwater emerging at the south end of DI. S2 showed an overall decrease in discharge over the monitoring period, as well as the lowest average discharge of the springs (1.2 × 10⁻⁴ m³ s⁻¹). Drain discharge from the crop field north of the stream showed an average discharge of 3.17 × 10⁻⁴ m³ s⁻¹ corresponding to a contribution of about 7% of the stream gain of 4.58 × 10⁻³ m³ s⁻¹ for a 90 m stream section.

If all springs in the wetland were assumed to discharge to the stream individually, the average combined spring discharge would amount to 2.7 × 10⁻³ m³ s⁻¹, corresponding to 59% of the average stream gain. In September 2016, when discharge through the streambed was measured (see below), the total discharge from all springs amounted to 2.5 × 10⁻³ m³ s⁻¹, which is very similar to the average spring discharge measured throughout the full monitoring period. DI and RI are the only two springs visibly discharging directly into the stream and may receive water from one or more of the other springs. Therefore, these springs two (DI and RI) can be considered separately. Their combined discharge only accounts for about one third of the combined discharge of DW, PO, S1, and S2. Furthermore, the combined discharge from DI and RI equals 0.9 × 10⁻³ m³ s⁻¹ accounting for 20% of the gain in the stream.

4.4. Streambed and Bank Fluxes

For each VTP location, the vertical groundwater Darcy flux (q) was estimated using the in situ, locally determined thermal conductivity, which ranged from 0.44 to 2.70 $\mathrm{W}{\mathrm{m}}^{-1}{}^{°}{\mathrm{C}}^{-1}$ (avg. 1.88 $\mathrm{W}{\mathrm{m}}^{-1}{}^{°}{\mathrm{C}}^{-1}$) (Figure 7A). Estimated fluxes through the streambed are presented in Figure 7B, with VTP1 located furthest downstream and VTP11 furthest upstream (locations shown in Figure 1B).

Estimated q values always indicated upwards flow to the stream. The mean estimated q was 7.02 × 10⁻⁷ m s⁻¹. The maximum q of 1.6 × 10⁻⁶ m s⁻¹ was measured at VTP2, i.e., at the north bank and the center of the stream. Here, VTP2 had the highest observed ${\mathsf{\kappa }}_e$, reflected in the highest average q across the stream of 1.57 × 10⁻⁶. The lowest flux of 1.35 × 10⁻⁷ m s⁻¹ was estimated for the north bank profile of VTP7, which also had the lowest observed ${\mathsf{\kappa }}_e$. Accordingly, q varied by a factor of five.

Assuming a stream width of 3 m along the 90 m stream segment surveyed with VTP, the total direct groundwater discharge to the stream amounts to 2.08 × 10⁻⁴ m³ s⁻¹. The groundwater discharge to the stream only represents 4% of the stream gain. Adding the combined average discharge of all springs (59%) together with the groundwater discharge (assuming half comes from the south and the other half from the north), this accounts for 61% of the gain in discharge in the stream gain that comes from the southern wetland section.

For comparison, Darcy’s Law can be applied to calculate the bank seepage to the stream, using the average stream stage, and average hydraulic head and distance to the stream’s bank for piezometers TH6 (north) and TH8 (south). The distance used is the piezometer distance to the center of the stream minus 1.5 m, assuming the stream is 3 m wide. The K-value measured in TH6 (20.6 m day⁻¹) is used on the north bank, while the measured K-value in TH7 at 2 MBGL of 31 m day⁻¹ is used on the southern side. The flow area used for both banks is the 90 m stream reach times a water depth of 0.35 m. The potential bank seepage to the stream from the northern bank is then 3.85 × 10⁻⁴ m³ s⁻¹ and through the southern bank 3.97 × 10⁻⁴ m³ s⁻¹. Potentially, this would amount to nearly 17% of the total stream gain, twice the amount, from each individual bank, than the flux estimated through the stream bottom using VTP-data.

4.5. Water-Stable Isotopes

Stable oxygen and hydrogen isotope composition of precipitation for a 2-year period from a nearby station (Voulund, UTM-E: 510000, UTM-N: 6210000) has previously been reported [28]. The δ¹⁸O values of rainfall at Voulund had a weighted average of −7.69‰ and ranged from −5.5 to −10.0‰. δ²H values had a weighted average of −51.5‰ and ranged from −38 to −73‰. Stable isotope concentrations, from all piezometers and for all samples collected in the monitoring period, yielded an average δ¹⁸O value of −7.79‰, ranging from −8.27‰ to −7.20‰. δ²H had an average of −51.5‰ ranging from −47.3‰ to −54.0‰.

Figure 8 shows the δ¹⁸O value of springs, the stream and the drain (DR) on the crop side of the stream as a function of time. About 15% of the samples had δ¹⁸O concentrations outside the δ¹⁸O mean groundwater signal ± 1 standard deviation, indicated in Figure 8 by grey shading. The springs include three surface-water bodies with free-flowing water (DI, RI and PO) where the isotope signal can be affected due to evaporation or precipitation. Indeed, DI and RI had δ¹⁸O values outside the range of the groundwater; RI was more was more enriched in May 2017 and the same was the case for DI in January, February, and April 2017. Both DI and RI, and the stream samples were more depleted in January. Direct precipitation, upwelling groundwater, evaporation, and their relative magnitude, determines the δ¹⁸O and δ²H composition at these three locations in the wetland.

DI δ¹⁸O values were variable the first four months after which they became steadier until the last measurement in January 2018. DI had the two most enriched δ¹⁸O values measured (January and April) and had a δ¹⁸O range of −8.67‰ to −6.63‰, which is the largest range for all the springs. The rivulet (RI) had the second largest range (−8.12‰ to −7.43‰) with a similar pattern to that of DI. The δ¹⁸O range of the PO draining the pond was −7.98‰ to −7.63‰, which could indicate enough residence time in the pond, to buffer changes due to variations in groundwater δ¹⁸O input.

The δ¹⁸O of the stream (ST) and the drain (DR) followed the same pattern with the same peaks and troughs. Water from the drain was, for all months except for March, enriched compared to the stream. S1, S2 and DW had the lowest fluctuations with a δ¹⁸O range of −7.99‰ to −7.48‰, −7.82‰ to −7.56‰, and −8.04‰ to −7.56‰, respectively. This can also be seen from the low standard deviations (legend, Figure 8). The δ¹⁸O ranges of S1, S2, and DW were all lower than what were measured for the groundwater at the field site.

Samples collected of the shallow groundwater in June 2017 under the wetland, using mini-piezometers and analyzed for δ¹⁸O are shown in Figure 9A. Average δ¹⁸O was −7.69‰ ranging from −8.06‰ to −7.11‰ and average δ²H was −51.06‰ ranging from −48.27‰ to −52.84‰. The spatial patterns of the measured δ¹⁸O show enriched signals in the center of the wetland—whereas depleted signals, compared to the average δ¹⁸O signal reported by Müller et al. [28], were mainly found along the hillslope to the southwest and 10–15 m from the stream (except for one measurement with a δ¹⁸O of −7.11‰). Water from DW, PO, and S1 discharge to the surface of the wetland and flows into the center of the wetland, where the peat has its maximum thickness. In this area samples from six of the mini-piezometers showed slight enrichment (−7.37 to −7.66‰, yellow area in Figure 9A).

From the piezometric head observation, it was established that the groundwater table in the wetland was below the ground surface, throughout the monitoring period. This suggests that water, originating from precipitation and the springs, on top of the peat represents a perched water table resulting in a surface-water-saturated soil. In addition, a downward gradient was present between piezometers TH9s and TH9d. The enriched isotope signals in the center of the wetland could therefore indicate water enriched from exposure to evaporation percolating through the peat and recharging the groundwater system.

A cross-section of the field site (Figure 9B) showing all samples collected for δ¹⁸O in both piezometers and mini-piezometers have been divided in to three broad sub-groups and interpolated by hand. A plume of enriched water is seen near the surface along the entire cross-section, except for directly under the stream where δ¹⁸O in groundwater is depleted. Samples taken under the streambed (−8.2‰) and from TH6 (−7.87‰ SD 0.17‰) and TH8 (−7.87‰ SD 0.05‰) right next to the stream showed consistently depleted signals. Depleted isotopic groundwater also enters the wetland from the southwest and either flows to the stream as groundwater or initially as surface flow. The surface flow is enriched and the observations of an enriched signal from the mini-piezometers indicate that water re-infiltrates. This also fits with the observed downward gradient between the nested piezometers TH9s and TH9d. Piezometer TH9d, sampled every meter at the time of installation, had enriched values of δ¹⁸O in the top of the well, but samples became more depleted with depth. δ¹⁸O concentrations from TH7d, sampled with depth in October 2016, were slightly more enriched and variable with depth. Downward seepage through the peat upstream TH7 followed by the apparent upward flow due to upward gradient between TH7d and TH7s, maybe the cause of a more uniform distribution.

Samples taken from TH10 at the hillslope (March and November 2017) also showed depleted δ¹⁸O (−7.91‰ and −8.02‰, respectively). The depleted δ¹⁸O values measured along the hillslope in mini-piezometers and in TH10 are similar to those measured in the two deepest piezometers; TH1d (16 MBGL) and TH2d (12 MBGL) on the crop side. Stable isotopes were sampled five times from TH1d, November 2015, February and September of 2016, and March and November of 2017, with an average δ¹⁸O of −8.02‰ and a δ¹⁸O minimum of −8.15‰ and maximum of −7.82‰. TH2d was sampled on four occasions in November 2015, February and September of 2016 and in March 2017 with an average δ¹⁸O of −8.11‰ with maximum and minimum of −7.96‰ and −8.26‰, respectively. Provided the isotope signals measured in TH1d and TH2d are representative of well-mixed groundwater, and not event-based signals, it would appear that a groundwater source with a similar signature is entering the wetland from the southwest. Some of this water bypasses the peat and flow to the stream through the 5–8 m thick gravelly/coarse sand. Deeper groundwater may also flow upwards to the stream as indicated by the depleted δ¹⁸O signals in the deeper parts of TH1 and TH2.

Piezometer TH3 had enriched δ¹⁸O values with depth indicating a downward flux (recharge area); the same is seen with TH4. From TH1d and TH2d δ¹⁸O signals became depleted with depth and from five MBGL and seven MBGL δ¹⁸O signals became steadier. The δ¹⁸O data from the crop side therefore indicates vertical downward seepage. In addition, upward leakage from the bottom of the profile cannot be ruled out.

5. Discussion

The main objective of the study has been to understand and quantify the heterogeneity in water fluxes to riparian zones and streams in a lowland catchment and specifically address the relative importance of groundwater-fed surface flow to subsurface flow directly to the stream. We will here try to discuss our results in relation to existing conceptual models of riparian zones e.g., [6,12,29,30]. For example, Vidon and Hill [30] proposed a conceptual model, where aquifer thickness and slope are two characteristics determining the link between riparian zones and the upland. Several studies [6,12,29] focus on how flow paths in typical riparian zones are defined by the landscape type, riparian zone–upland link, riparian aquifer thickness, and permeability of a confining peat layer. However, only one of these studies includes the geomorphology of the stream (width, sinuosity) [12].

Most conceptual flow path models emphasize different parts of the link between stream–riparian zone regional aquifer, but with no guidelines to evaluate the relative contribution of subsurface and groundwater-fed surface flow. The size of the sub-catchment to discharge station Q₂ at our field site is approximately 28 km² with roughly 23 km of stream network. Precipitation and actual evapotranspiration average 984 mm year⁻¹ and 510 mm year⁻¹, respectively, giving an average recharge of 474 mm/year [15]. Average stream flow Q₂ was measured in 2017 to ~10⁷ m³ year⁻¹ corresponding to a runoff of 377 mm/year. Thus, roughly 20% of recharge is not captured by the stream network, but flows out of the sub-catchment via groundwater. If it is assumed that all the stream flow originates as base flow (direct groundwater discharge) then one can compute the average flux. Assuming the average width is ~1.5 m then the 23 km network has a streambed area of 34000 m². This gives an average streambed flux of 1.1 m day⁻¹. Clearly, actual rates in this type of catchments are orders of magnitude lower; we measured 0.06 m day⁻¹. Accordingly, a lot of water is fed into the stream valley and only a small portion is captured by the stream directly though the streambed. Flow to the stream therefore takes alternative pathways.

5.1. Hydraulic Connectivity

Features such as upland linkage, topography, and hydrogeological setting, both adjacent to and in the riparian lowland, are known to influence the flow distribution in riparian lowlands [29,30]. Landscape characteristics such as slope of riparian zone and depth of upland permeable sediments influence the hydraulic connectivity, water table fluctuation and groundwater flow direction [30]. In the study area, both to the north and south, relatively thick unconfined aquifers with a gently (north, ~2.5%) too steeply (south, ~7%) sloping topography are present. Vidon and Hill [30] have shown that such landscape and hydrogeological characteristics would give the southern riparian lowland a permanent upland aquifer linkage with a constant flow direction, a continuous large flux, and small water table variations. For the northern riparian lowland, the upland–riparian link to the upland aquifer would still be permanent, but with small to moderate fluxes and moderate water table variations. This is confirmed by our measurements. For example, the largest fluctuations in hydraulic head was found near the edges of the riparian lowland, but notably with smaller fluctuations on the south side (TH10 with hydraulic head variation of ±0.2 m) compared with the north side (F100 ±0.37 m). Nevertheless, generally, observations showed small variations in the groundwater table, indicating a permanent link to an upland aquifer. This ensures a constant groundwater input and maintains a year-round hydraulic gradient towards the stream. The largest fluctuations observed were to the north, indicating a linkage to a smaller upland aquifer. Considering that the upland aquifer thickness is slightly less and the length to the upland boundary is shorter compared to the southern side of the stream, this could explain the slightly larger fluctuations. To the south, the connectivity of the riparian lowland appears to maintain a more continuous groundwater flux giving less seasonal variation near the hillslope.

The lack of identifying a confining layer at the bottom of the riparian aquifer suggests the potential of deeper groundwater flow paths bringing deeper and more regional water to the riparian zone from further away in the catchment. The δ¹⁸O signal (Figure 9) confirms this, where measurements below the streambed were close to the meteoric average.

5.2. Heterogeneity in Flow Paths to Stream

Conceptual models of flow paths in riparian zones as a function of the link to an upland aquifer and the thickness of the riparian aquifer have been proposed [6]. This model have since been extended to include groundwater-fed surface flow [8]. A link to larger upland aquifer results in large groundwater inputs and small water table fluctuations (see above) [28]. With large volumes of water transferred to the riparian aquifer and with a limited area of discharge more water is routed to the surface as seeps and springs and reach the stream via overland flow.

The studied riparian aquifer consists of sand. On the southern side of the stream, a semi-confining peat layer is present on top of the sand. The depth of the aquifer increases from ~15 m near the hillslope to >30 m in the central part (a deeper confining layer of the riparian aquifer has not been identified from boreholes or the ERT profiles). The large thickness and the permanent link to a more regional aquifer fits well with other conceptual models [6,8] where a significant part of the regional flow to the riparian zone is routed to the stream via seeps and springs—despite the considerable thickness of the aquifer.

This was best observed on the south side of the stream. For example, the observed diffusive seepage occurring in the area above the pond. Furthermore, the presence of springs and the more extensive need for drainage (ditch, DW) is evidence of a larger groundwater to surface-water component. This is not the case on the north side and may be explained by a link to a much smaller regional aquifer. Here flow to the stream is mainly governed by groundwater, except for the part captured by the drain near the depression in the landscape.

Like in other studies there is an exchange of water between flow paths (e.g., [8,31]). TH9s and TH9d on the southern side of the stream show a small downward gradient. TH9s is screened just below the peat layer and TH9d is screened in the coarse sand deposited under the peat located approx. 2–6 MBGL (Figure 2). The more conductive coarse sand could produce a downward gradient; generating a preferential flow path to the stream under the peat. Surface flow could therefore infiltrate again (Figure 9), flow towards the stream and upwards due to the upward gradient at TH7s and TH7d (TH7s is screened in the high conductive sand and TH7d screened below in a less conductive material) and the presence of the stream. The enriched δ¹⁸O plume near TH9 and TH7 confirms this.

Part of the groundwater will pass through the aquifer to the stream from both the north and the south (Figure 9). The individual fluxes estimated by VTP across a 90 m stream reach vary with close to an order of magnitude. Groundwater discharge has been demonstrated by other local studies, [16,17,22] and elsewhere [32,33] to vary substantially in time and space both across and along a stream reach. Compared to these local studies further downstream of Holtum stream, the estimated streambed fluxes in this study is at the lower end (average of 7 × 10⁻⁷ m s⁻¹). Two local studies [16,34] reported mean fluxes of 4–5 × 10⁻⁶ m s⁻¹ approximately a factor 6 higher on average than the fluxes measured at our site. This corresponds nicely with the observations made in one local study [16] where it was found that stream flow in areas downstream of our site was less rain-event controlled. The difference between our site and the other sites in Holtum is that our site is, or rather was, managed to some degree and now appears with overgrown ditches, overflowing drainage wells, and a possible broken drainage system.

In September 2016, the overland flow component was a factor of 12 higher than discharge estimated through the streambed. This clearly indicates that groundwater-fed surface flow is the most important flow component at the site. This ratio may change over the season. We do not have the data to quantify this as discharge through the streambed was only measured once.

Langhoff et al. [12] proposed a coefficient C called “the apparent relative seepage width” to indicate the proportion of water entering the stream through the streambed relative to overland flow. C is defined as the width of the wet zone of the flood plain (b_r) divided by the effective width of the stream (Sb_s), i.e., C = b_r/(Sb_s). Here, S is the stream sinuosity (which for Holtum stream is S = 1, i.e., a straight stream) and b_s is the width of the stream. High values of C correspond to relative low contribution to stream flow gain from direct seepage through the streambed and vice versa. At our field site, the wet zone is on the south side (b_r 85 m) and with a stream width of b_s = 3 m, C can be computed as ~28. This is similar to Site 3 studied in Langhoff et al. [12], who also found that the majority of the water flowed to the stream as overland flow (actually 110%). It thus appears that the C ratio can be an efficient way of quantifying heterogeneity in flow paths, because it takes into account the stream valley morphology and wetness.

5.3. Water Balance

The gain in stream discharge between Q1 and Q2 was estimated to be on average 5.09 × 10⁻⁵ m³ s⁻¹ m⁻¹. Higher discharges were attributed to higher intensity of precipitation, for where the largest gaining and losing conditions in the stream also was prevalent—these (four) measurements was not included in the average stream gain. The stream gain had a standard deviation of 7.6 × 10⁻⁵ meaning that the stream between Q1 and Q2 can experience both losing and gaining conditions. However, only two of the discharge measurements represented losing conditions and the stream section between Q1 and Q3 can be classified as gaining. The ADC is expected to have an uncertainty of 5% of the measured discharge [35]. This limits the minimum distance between two measurement stations when using the ADC, as the accretion of a discharge should be higher or less than the uncertainty—the average accretion for our measurements was 6%, excluding months where discharge measurements was known to be affected by precipitation. It would have been better to have continuous (daily) measurements of stream discharge, as this would have reduced the uncertainty on the stream gain. However, as explained above, it was not possible to at the same time measure continuous (daily) spring flows. This means that the water balance is to be regarded as the average of monthly snapshots of stream flow gains, spring flows, and only one measurement of direct seepage through the streambed (which turned out to be minor).

The estimated gain between Q1 and Q2 was significantly less than that estimated from Q1 to Q3. This could be explained by the large portion of arable land with drain tiles to the south and the riparian lowland with several ditches and a small tributary between stations Q2 and Q3 (Figure 1C). The dynamics of the stream stage followed the fluctuations of precipitation (flashy) and that of the hydraulic heads. This might reflect the riparian zone with conductive layers; therefore, it quickly drains and returns to its state of equilibrium. The discharge measurement of August, where discharge at Q2 was measured after a rainfall event, showed a significant decrease in discharge compared to the measurements at Q1 and Q3 made during the event. The stream therefore responds quickly to rainfall events. This is supported by the findings of Poulsen et al. [16] who estimated that the contribution to stream discharge from rainfall events for the same sub-catchment as studied here could be as high as 65%.

The average estimated flux though the streambed can only account for 4% of the stream gain. Langhoff et al. [12] measured bank seepage to a stream that had been deepened and found that bank seepage could contribute of up to 25% of the stream gain. With a similar setting at the field site, where the streambed is situated at a depth below the peat layer to the south, it is likely that bank seepage plays a moderate to important role. It was calculated based on data from the two nearest piezometers that bank seepage potentially could contribute with 17% of the stream gain.

The drain located on crop field discharges directly into the stream contributing about 7% to the stream gain. The most prominent contribution to the stream gain is water discharging through the springs in the wetland to the south. The estimated contribution of measured spring discharge to the stream gain was from 20% and up to 59% (depending on whether they contribute individually or not).

With a potential of close to 66% of the stream gain coming from focused sources a contribution of 4% and 17% from direct groundwater discharge to the stream seem insignificant. With up to 87% of the stream gain accounted for the water balance is not closed. Diffusive overbank flow, unmeasured discharge from (one known) submerged drains or highly focused bank seepage could account for the missing 13%. These numbers are similar to those found in the study by Langhoff et al. [12]. Here the measured or unmeasured bank seepage and overland flow are of similar magnitude. The results emphasize the important role of groundwater-fed overland flow [8].

We cannot exclude the possibility that some of the springs are unnatural in that they could be the result of a broken drainpipe system. We know of the existence of the (overflowing) drainage well (DW) and the submerged drain outlet to the stream. This point to the potential role of an unmanaged drainage system of the hydrology of a wetland system. It has not been possible to obtain information about the drainage system, but the landscape with a steep increase in terrain to the south suggests that a drain network is present only in the wetland and does not drain the sandy/gravelly agricultural fields located further south (see e.g., Figure 2 with well DGU 92.2187). The drains thus capture groundwater that probably otherwise would have resulted in groundwater-fed overland flow. The human alteration of the system (channelization, drain pipes, ditches) has potentially increased the lateral connectivity of the stream and the terrestrial landscape [36]. As mentioned above, the fluxes through the streambed is roughly six times lower compared to other more natural (unmanaged) systems downstream of our site. Our results are also different to those reported by Frederiksen et al. [13] for a groundwater-dominated stream west of our site. Here overland flow was negligible, and groundwater was the main flow path (>60% of stream flow gain). A drainage system (broken or not, ditch) thus have re-organized the flow paths and most groundwater enters the stream by groundwater-fed surface flow than groundwater-fed stream flow. This will of course also have implications for nitrate transport and removal, which is currently being studied. With that said, some springs may be entirely natural and could be due to a thin peat cover (e.g., S1 and S2). For example, Johansen et al. [11] found that peat thickness controlled vertical seepage to a rich fen.

6. Conclusions

From the present study, it can be concluded that management for agricultural purposes have changed the natural flow paths through the riparian lowland to the stream. In the present study we have addressed the connection of a groundwater system and a riparian lowland to specifically determine the relative importance of subsurface and surface flow to a headwater stream.

The investigated riparian lowland had a permanent upland linkage ensuring continuous large groundwater inputs, small groundwater table fluctuations, and maintaining a constant direction of flow towards the stream. The permanent upland connectivity and aquifer thickness result in large groundwater inputs. Thus, not all the groundwater from the upland aquifer can be accommodated as seepage only through the streambed of the relatively small headwater stream. In our case we were able to demonstrate that the direct discharge through the streambed was approximately a factor of six lower compared to other riparian–stream systems in the same catchment (more natural riparian zones with no or little human alterations). In other words, the highly managed wetland system captures a high proportion of the regional groundwater discharging to the riparian zone (~2/3 of the regional flow bypasses the subsurface). This can explain the need for anthropogenic alterations of the area in the past—changes needed to make the land usable for agricultural production or as, in this case, a cow pasture. The anthropogenic changes include ditches, drainage wells, and the man-made pond. Most of the “springs” are therefore man-made preferential flow paths. The straightening of the stream has likely made these changes even more needed.

Such anthropogenic changes are common for many riparian lowlands situated in agricultural lowland catchments. Many of these are easily detectable from maps (stream straightening, ditches) and of course by site visits. To determine the relative importance of groundwater-fed overland flow to direct stream seepage is of course more complicated as this requires detailed field monitoring. This just highlights the importance of measuring all flow components in a water budget for a riparian–stream system. In the end, anthropogenic alterations will have consequences for the effectiveness of N-removal in the subsurface as the main part of groundwater is routed to the stream via drains, ditches and diffuse surface flow/rivulet with short residence times and little contact with subsurface reactive materials.