The Significance of Vertical and Lateral Groundwater–Surface Water Exchange Fluxes in Riverbeds and Riverbanks: Comparing 1D Analytical Flux Estimates with 3D Groundwater Modelling

Abstract

Riverbed temperature profiles are frequently used to estimate vertical river–aquifer exchange fluxes. Often in this approach, strictly vertical flow is assumed. However, riverbeds are heterogeneous structures often characterised by complex flow fields, possibly violating this assumption. We characterise the meter-scale variability of river–aquifer interaction at two sections of the Aa River, Belgium, and compare vertical flux estimates obtained with a 1D analytical solution to the heat transport equation with fluxes simulated with a 3D groundwater model (MODFLOW) using spatially distributed fields of riverbed hydraulic conductivity. Based on 115 point-in-time riverbed temperature profiles, vertical flux estimates that are obtained with the 1D solution are found to be higher near the banks than in the center of the river. The total exchange flux estimated with the 3D groundwater model is around twice as high as the estimate based on the 1D solution, while vertical flux estimates from both methods are within a 10% margin. This is due to an important contribution of non-vertical flows, especially through the riverbanks. Quasi-vertical flow is only found near the center of the river. This quantitative underestimation should be considered when interpreting exchange fluxes based on 1D solutions. More research is necessary to assess conditions for which using a 1D analytical approach is justified to more accurately characterise river–aquifer exchange fluxes.

1. Introduction

The interaction between rivers and groundwater plays an important role in a wide range of contemporary challenges, including the provision of drinking water, the characterisation and management of environmental flow regimes, the conservation or restoration of riverine ecosystem health and functions [1,2], and the attenuation of contaminants [3]. It is therefore crucial to understand and quantify exchange processes between rivers and groundwater for the best possible protection of our water resources.

One of the parameters describing groundwater–surface water interactions in terms of water flow is the exchange flux (Darcy flux) between rivers, riverbeds, and connected aquifers. Nowadays, a wide range of methods is available for the quantification of these fluxes (see [4,5]). These methods can broadly be subdivided into hydraulic, chemistry-based, and tracer approaches [6] and are applicable at different spatial and temporal scales. At reach and smaller scales, river–aquifer exchange is often estimated from information on the hydraulic gradient and hydraulic conductivity [7,8,9] as input to calculations making use of Darcy’s Law [10]. As an alternative, the application of environmental tracers (including heat, stable isotopes, major ions, or dissolved oxygen) has become well established. Here, the idea is to derive information on river–aquifer interaction by linking the quantity or change of a tracer to the movement of water [11,12]. For example, Cook [13] estimated groundwater discharge to rivers using individual ion concentrations, electrical conductivity, stable isotopes, and the dissolved gases helium, chlorofluorocarbons, and radon, while Cox et al. [14] applied a combination of tracers including heat, chloride, and electrical conductivity for a similar purpose.

The application of heat as a tracer has received high interest in the last two decades [15,16,17]. Heat as a tracer is a fast and cost-effective method [4] that can be used to obtain a high-resolution spatial coverage of point-in-space exchange flux estimates [18]. One of the advantages of heat as a tracer over hydraulic approaches is that subsurface thermal properties vary over a much narrower range than hydraulic properties, such as riverbed hydraulic conductivity, potentially reducing the uncertainty on the exchange flux estimates [11]. All heat-tracing techniques are based on the principle that the movement of water through a porous medium influences the temperature distribution in that medium [19,20]. The temperature distribution in a water-saturated riverbed is the result of the processes of advection caused by the movement of water through the pores and heat conduction through and between solid and liquid phases [16]. Exchange fluxes can then be estimated, e.g., by fitting a steady-state analytical solution to measured vertical riverbed temperature profiles [18,21,22,23], time-series analysis of riverbed temperature data [24,25,26], or by coupled 3D groundwater flow and heat transport modeling [27,28,29,30]. Various 1D steady-state [31,32] and transient-in-time solutions [33,34,35,36,37] to determine vertical exchange fluxes have found widespread application and have been implemented into various software packages [38,39,40,41].

One of the main limitations of the 1D methods with respect to estimating representative exchange fluxes is the assumption of strictly vertical flow [28,42] in a homogeneous riverbed. The importance of both vertical and horizontal flow processes in riverbeds and their impact on flux estimates is the subject of ongoing research [43,44,45]. Although the exact contributions of vertical and non-vertical flow components to exchange fluxes are site-specific, previous research has shown that fully vertical flow is most likely to occur only beneath the center of a river at rather shallow depths. Near its banks and with increasing depth, horizontal and lateral flow components increase [46,47]. It also seems that for very small vertical velocities the 1D analytical solutions tend to overestimate vertical fluxes [48].

The contribution of vertical and non-vertical flow components to exchange fluxes is closely related to riverbed heterogeneity [49]. Even at the meter-scale, riverbeds can be highly heterogeneous structures [50,51], characterised by complex 3D flow fields [52]. This heterogeneity often results in significant spatial variability in temperature, flow, and discharge patterns in riverbeds [53,54,55]. As such, a proper characterisation of riverbed heterogeneity can prove beneficial when describing river–aquifer exchange flux patterns at the meter-scale.

Few studies have used temperature measurements and hydraulic data to uncover reach-scale spatial patterns of groundwater–surface water exchange fluxes in heterogeneous riverbeds. For example, Conant [9] investigated the spatial distribution of groundwater discharge along a 60 m long section of the Pine River in Ontario, Canada, and found an empirical relationship between Darcy fluxes derived from hydraulic head and conductivity measurements and riverbed temperatures. Schmidt et al. [18] then used these river and riverbed temperatures mapped by Conant [9] to calculate river–aquifer exchange fluxes with a steady-state analytical solution and demonstrate their spatial variability. In another study, Lautz and Ribaudo [56] combined a time-series of vertical temperature profiles with point-in-time riverbed temperatures to show the spatial variation of exchange fluxes over a 30-m-long reach of Ninemile Creek, NY, USA. Munz, et al. [57] used information on streambed temperatures, water levels, and hydraulic conductivities to determine spatial and temporal variations in flow patterns around geomorphological structures along a reach of the Selke River, Germany.

By using spatially distributed fields of riverbed hydraulic conductivity [58,59], physically based, mechanistic groundwater models allow for the quantification of complex reach scale 2D/3D flow processes. The novelty of our study lies in the combination of high-resolution riverbed hydraulic conductivity measurements with 3D groundwater modelling to simulate complex flow patterns in the riverbed at the meter-scale. These high-resolution flows are compared with exchange fluxes estimated with a 1D heat tracer method. We hypothesise that this comparison can help us (i) better understand the complex spatial flows in riverbeds; (ii) highlight the impact of riverbed heterogeneity on the spatial variability of exchange fluxes; and (iii) advance our understanding of heat tracer-based method applicability in actual field settings.

Here, we study the spatial distribution of river–aquifer exchange fluxes at the meter-scale for two sections of the Aa River, Belgium, which are obtained with two different methodological approaches—(i) the heat tracer method and (ii) numerical groundwater flow modelling, respectively. Groundwater–surface water interaction along the Aa River and its streambed characteristics have been studied previously, using heat as a tracer [21,60,61], in-stream piezometers [24], hydraulic conductivity measurements [51,62], electrical resistivity tomography, and induced polarisation [63] in addition to regional [64] and local reach groundwater modelling [58,65]. In this study, we first estimate point-in-time vertical exchange fluxes from 115 riverbed temperature profiles using a simple 1D analytical solution to the heat transport equation. We then combine the extensive dataset of vertical flux estimates with variogram analysis and ordinary kriging to identify the spatial variability and anisotropy of the flux and generate spatially distributed maps of the river–aquifer exchange. We compare these spatially distributed flux estimates to estimates obtained from a numerical 3D groundwater model (set up in MODFLOW) that essentially quantifies fluxes based on Darcy’s law. This model makes use of previously obtained high-resolution data on riverbed hydraulic conductivity [51] to simulate riverbed heterogeneity.

2. Study Area and Field Measurements

Fieldwork was carried out along two sections of the Aa River, a typical lowland river located in northeast Flanders, Belgium (Figure 1). It is a tributary of the Kleine Nete River and has a drainage area of 237 km². In the study area, it has an average discharge of about 1.9 m³/s, an inclination of 0.48‰, and a manning coefficient of 0.06 [60,66]. The area is characterised by a temperate maritime climate with mild winters. The mean annual rainfall is 830 mm, with the monthly average rainfall fluctuating between 55 mm in January and 90 mm in July. The mean annual evapotranspiration is 670 mm [67]. The study area has a mean annual temperature of 10.5 °C, with respective values of 17.5 °C for summer and 3.6 °C for winter [68]. The dominant land use in the catchment is agriculture. The study area is located at a 1425-m-long canalised river reach, in which water levels are controlled by an upstream and a downstream weir (designated as Weir 3 and 4, respectively; Figure 1). These weirs help to regulate river discharge, resulting in relatively constant water levels throughout the year, which is favourable for the agricultural areas along the river course.

Two sections were chosen in the study area, based on earlier research of Anibas et al. [21,60]. This allowed us to compare the impact of temperature and hydraulic conductivity data on flux estimates along the river course for partially distinct (hydro)geological settings. The downstream section is located 200 m upstream of Weir 4, while the upstream section is located 185 m downstream of Weir 3. Both sites are approximately 900 m apart (Figure 1).

The geology of the local aquifer connected to the Aa River in the study area consists of sandy aquifer deposits approximately 80 m in thickness, limited by the underlying Boom aquitard, a thick clay deposit. The study area is covered by two stratigraphical units, the Kasterlee Formation and the underlying Diest Formation. According to Louwye et al. [69], both formations can be distinguished by their grain size and mineral content. The Diest Formation is a coarse, loosely compacted upper Miocene sand deposit with high glauconite content [70]. The Kasterlee Formation has low glauconite content but is highly compacted. Geological maps [71] suggest a boundary between the two formations within the study area so that the Kasterlee Formation is the top layer present below the Quaternary cover at the upstream section and the Diest Formation at the downstream section. However, the exact location of this boundary in the study area is unknown.

The riverbed is mainly composed of fine sand with a highly variable fraction of organic matter present on its top. This organic layer is more pronounced near the banks due to the presence of vegetation [51]. For additional information with respect to the field site, refer to, e.g., Anibas et al. [21], Benoit et al. [63], and Ghysels et al. [51].

The downstream section consists of a 25-m-long river stretch with an average width of 15 m. Riverbed bathymetry is relatively even, and the average water depth is shallow and around 0.6 m (Figure 2). The upstream section consists of a 14-m-long river stretch with a width of approximately 12 m. Here, the bathymetry is more heterogeneous, with an average water depth of 1.0 m. In the center of the river, the water depth is higher (up to 1.5 m; Figure 2). The differences in bathymetry are possibly related to the influence of the upstream weir which has a noticeable effect on river morphology. There, peak discharge of up to 10 m³/s is observed at times, which may result in scouring of the riverbed in the center of the river. The organic matter layer is present on top of the sandy riverbed sediments near the banks, while this layer is absent in the center of the river. At the downstream section, some organic matter is also present near the banks, albeit less pronounced than at the upstream section.

Ghysels et al. [51] performed an extensive field campaign focusing on the meter-scale spatial variability of riverbed hydraulic conductivity ($K$) at the two studied sections (Figure 2). In total, 220 rising-head slug tests and 45 falling-head standpipe tests were conducted to estimate horizontal and vertical riverbed $K$, respectively. Both horizontal and vertical $K$ range over several orders of magnitude and show significant meter-scale spatial variability. For both sections, horizontal $K$ varies between 0.1 m/day and 23.2 m/day, with an average of 6.4 m/day (n = 220). Vertical $K$ at the upstream section varies between 3.0 × 10⁻³ m/day and 2.0 m/day, with an average of 0.4 m/day (n = 37). At the downstream section, vertical $K$ varies between 0.1 m/day and 4.0 m/day, with an average of 1.6 m/day (n = 8). Elongated zones of high horizontal $K$ along the river course were observed at both sections, indicating a link between riverbed structures, depositional environments, and flow regimes. The spatial variability of vertical $K$ seems to be governed mainly by the presence and thickness of the low permeability organic layer at the top of the riverbed. At the upstream section, this organic layer is well-developed, resulting in significantly lower vertical $K$ compared to the downstream section, where less organic matter is present. The lower vertical $K$ observed in the center of the river at the upstream section is related to the scouring of the top of the riverbed during peak discharges due to the lowering of the upstream weir. The spatial variability of the riverbed $K$ datasets was quantified by modelling omnidirectional and directional (along and across the river) variograms. Both display a clear directional anisotropy with variogram ranges along the river being twice as large as those across the river course.

Measurements of the hydraulic conductivity of the river banks ($K_{bank}$) were performed at the downstream section by Baya Veliz [62]. $K_{bank}$ was measured using L-shaped standpipes [72]. Estimated $K_{bank}$ varies between 0.5 m/day and 16.2 m/day, with an average of 5.2 m/day (n = 8 at each bank). Mean $K_{bank}$ is more than 50% lower than $K$ of the surrounding aquifer material.

3. Methodology

3.1. Temperature Measurements

Vertical riverbed temperature profiles were collected with a ‘Temperature lance’ (T-lance). The T-lance consists of a metal rod with a T-shaped handle at the top and a pointed tip at the bottom. Near the bottom end, a slotted hole holds a single thermistor (Davis Instruments Model 7817; Hayward, CA, USA). The electric resistance is measured with an electric multimeter and converted to temperatures by using a calibration function [21].

Using the T-lance, 115 temperature profiles were collected in roaming surveys at the two measurement sections, consisting of 610 individual temperature measurements. At the upstream section, 64 temperature profiles (Figure 2) were obtained during 14–16 July 2015, consisting of 376 temperature measurements. At the downstream section, 51 temperature profiles (Figure 2) were logged on 17 August, 26 August, and 27 August 2016, consisting of a total of 234 temperature measurements. At the upstream section, temperatures were measured at depths of 0.0 (i.e., the interface between river and riverbed), 0.1, 0.2, 0.3, 0.5, and 0.7 m within the riverbed, in addition to the maximum depth accessible with the ‘T-lance’. At the downstream section, temperatures were measured at depths of 0.0, 0.1, 0.2, and 0.5 m, in addition to the maximum accessible depth.

3.2. Vertical Exchange Flux Estimates from Riverbed Temperature Profiles

The measured temperature profiles were used to calculate vertical river–aquifer exchange fluxes considering a 1D heat and fluid transport problem. The applied model assumes a thermal and hydraulic steady state, for which the 1D, vertical, anisothermal transport of heat through a homogeneous, porous medium (the riverbed) can be described as

$$k\frac{{\partial }^{2}T}{\partial z^2}-q_z{c}_w{\rho }_w\frac{\partial T}{\partial z}=0$$

(1)

where $k$ is the thermal conductivity of the soil–water matrix in J/(s·m·K), $T$ (K) is the temperature at point $z$ at time $t$ in the sediment, $c_w$ is the specific heat capacity of the fluid in J/(kg·K), ${\rho }_w$is the density of the fluid in kg/m³, and $q_z$ is the vertical component of the Darcy velocity in m/s [21,73,74]. A positive $q_z$ corresponds to water flowing from the river towards the aquifer (groundwater recharge), while a negative $q_z$ indicates the reverse flow direction. In our study, the vertical Darcy velocity $q_z$ is expressed in terms of mm/day.

Bredehoeft and Papadopulos [31] derived an analytical solution to Equation (1), with boundary conditions $T_z=T_0$ at $z=0$ and T_z = T_L at $z=L$, where $T_z$ is the temperature at any depth $z$, and $T_0$ is the temperature at $z=0$ (i.e., the uppermost temperature measurement), $T_L$ is the measurement at $z=L$ (i.e., the lowermost temperature measurement), and $L$ is the vertical distance over which temperatures are being observed. After solving Equation (1) with the specified boundary conditions, the following solution is obtained:

$$\frac{T_z-T_0}{T_L-T_0}=f(\beta ,z/L)=\frac{e^{\beta (z/L)}-1}{e^{\beta }-1}$$

(2)

with

$$\beta =\frac{q_z{c}_w{\rho }_{w}L}{k}$$

(3)

where $\beta$ is a dimensionless parameter that is estimated by minimising the objective function $F\left(\beta \right)$,

$$F\left(\beta \right)={\displaystyle \sum }_{z=0}^{z=L}{\left(\frac{T_z-T_0}{T_L-T_0}-\frac{e^{\beta (z/L)}-1}{e^{\beta }-1}\right)}^2$$

(4)

The vertical Darcy velocity $q_z$ can then be determined by rearranging Equation (3),

$$q_z=\frac{k\beta }{c_w{\rho }_{w}L}$$

(5)

To solve this problem, a thermal steady state is assumed. Seasonal influences on the temperature measurements can be approximated by a sine function with a wavelength of one astronomic year. According to Anibas et al. [61], a steady-state assumption is valid close to the function’s maximum (summer) and minimum (winter). Our measurements were performed in summer making a steady-state assumption appropriate. Diurnal temperature variations can also have an influence on vertical flux estimates [75]. However, the high-frequent diurnal temperature signal has a limited penetration depth in the riverbed. It is therefore advisable to exclude measurements close to the surface water-groundwater interface from further analyses when using a steady-state solution [76]. The diurnal temperature changes at the Aa River are in the range of 2–5 °C in summer. Anibas et al. [21] detected a significant diurnal influence only in the uppermost 0.1 m of the riverbed. The temperature at a depth of 0.1 m is therefore used as an upper boundary condition. Vertical flux estimates will thus not contain the portion of upwelling hyporheic flux in the first 0.1 m of the riverbed that can be caused by small-scale flow over roughness features on top of the riverbed. For an in-depth discussion of the contribution of hyporheic versus groundwater flow to total flux see, e.g., Bhaskar et al. [77].

The solution after Bredehoeft and Papadopulos [31] has been applied previously in studies of groundwater–surface water exchange by, e.g., Anibas et al. [21,61], while Schmidt et al. [18] used a similar solution to estimate fluxes from point-in-time measurements. Here, it is used to fit the 115 temperature profiles measured at both sections of the Aa River following Arriaga and Leap [78]. All physical properties and boundary conditions used to solve Equation (5) are shown in

Table 1. Physical properties and boundary conditions that are used to solve the analytical solution of Bredehoeft and Papadopulos [32].

Parameter	Value	Unit
Thermal conductivity k	1.8	J/(s·m·K)
$\mathrm{Density}\mathrm{of}\mathrm{water}{\rho }_w$	1000	kg/m3
$\mathrm{Specific}\mathrm{heat}\mathrm{capacity}\mathrm{of}\mathrm{water}{c}_w$	4183	J/(kg·K)
$\mathrm{Upper}\mathrm{boundary}\mathrm{condition}{T_0}^{\mathrm{a}}$	Variable	°C
$\mathrm{Lower}\mathrm{boundary}\mathrm{condition}{T_L}^{\mathrm{b}}$	11.2	°C

^a Located at the interface between surface water and riverbed (depth of 0.1 m). ^b Located at a depth of 5 m.

. We assumed a thermal conductivity value of 1.8 J/(s·m·K), taking into account the organic content of the riverbed sediments [21]. Thermal properties are not expected to have a large influence on flux estimates at our study site due to their narrow range.

3.3. Variogram Analysis and Spatial Interpolation

An exploratory variogram analysis is carried out to study meter-scale spatial variability and directional anisotropy in the exchange flux dataset. A variogram $\gamma \left(h\right)$ can be used to quantify the spatial variability in a dataset and is defined as [79]

$$\gamma \left(h\right)=\frac{1}{2N\left(h\right)}{\displaystyle \sum }_{i=1}^{N\left(h\right)}{\left(x_i-x_i^{\prime }\right)}^2$$

(6)

where $N\left(h\right)$ is the number of data pairs, while $x_i$ and $x_i^{\prime }$ represent one data pair separated by the lag distance $h$. Geostatistical interpolation techniques, such as kriging, require a variogram model to reproduce the spatial variability in the dataset [80].

Both omnidirectional and directional (across and along the river course) variograms were computed using the Stanford Geostatistical Modelling Software (SGeMS, Stanford University, Stanford, CA, USA) [81]. The experimental variograms were manually fitted by variogram models that consist of the sum of a nugget and a spherical model

$$\gamma \left(h\right)=\{\begin{array}{c}0 if h=0\\ C_0+ C_1\left(\frac{3h}{2a}-\frac{1}{2}{\left(\frac{h}{a}\right)}^3\right) if 0<h \le a\\ C_0+ C_1 if h>a\end{array}$$

(7)

where $C_0$ is the nugget effect, $C_1$ is the sill minus $C_0$, and $a$ is the range. Based on the resulting variogram parameters, spatially distributed fields of vertical flux estimates were interpolated using ordinary kriging.

3.4. Groundwater Flow Modelling

In this study, we could take advantage of an extensive dataset of riverbed $K$ collected by Ghysels et al. [51], which is briefly described in Section 2, based on which spatially distributed fields of riverbed $K$ were generated.

A 3D groundwater flow model (MODFLOW) of the entire stretch of the Aa River shown in Figure 1 had been set up and modified in previous studies [58,64,65]. This larger reach-scale model is 1500 m by 1800 m in size and has a horizontal resolution of 2.5 m by 2.5 m. The model domain is bounded by constant-head boundaries at all sides, based on the results of a regional model for the Kleine Nete River basin [64]. The bottom boundary is a no-flow boundary due to the presence of the Boom Aquitard. On top of the model, a recharge boundary is applied, with spatially distributed recharge calculated with the WetSpass model [82]. The aquifer over the model domain is subdivided into two layers—(i) Quaternary and Pleistocene cover deposits together with the Kasterlee Formation and (ii) the Diest Formation combined with the Berchem Formation. Model thickness is approximately 90 m.

Based on this larger reach-scale model, two fine-scale models are developed here, focusing on the upstream and downstream sections (Figure 1). Both models are 200 m by 200 m in size and have a horizontal resolution of 0.5 m by 0.5 m. A constant-head type boundary is assigned to all side boundaries, based on simulated heads of the larger-scale reach model. The same recharge as used in the larger-scale model is used. The two model layers of the larger-scale model are subdivided into 11 model layers in total for the two fine-scale models. This increased vertical resolution is needed for an accurate simulation of the river–aquifer interaction. Both fine-scale models are steady-state models, are set up using FloPy [83], and are run with MODFLOW-2005 [84].

To account for the influence of possible non-vertical flow components on the exchange flux, the approach presented by Ghysels et al. [58] is followed, according to which both vertical flow through the riverbed and horizontal flow through the riverbed and river banks can be simulated. The river water itself is modelled as a constant-head boundary with the hydraulic head equal to the river stage. To simulate river–aquifer exchange fluxes in detail at both upstream and downstream sections, the riverbed is divided into four separate layers of 0.25 m thickness each. The effect of a clogging riverbank layer is modelled with the Horizontal Flow Barrier (HFB) package. A thickness of the riverbank clogging layer of 0.5 m is assumed. This thickness was deduced from previously conducted standpipe tests [62] and is equivalent to the length of the horizontal section of the L-shaped standpipe. An average $K_{bank}$ value of 5.2 m/day is assigned to all river bank cells, obtained from standpipe measurements [62].

Spatially distributed fields of horizontal ($K_h$) and vertical ($K_v$) riverbed hydraulic conductivity are used as input for the groundwater models and resulting river–aquifer exchange fluxes and flow paths through the riverbed and riverbanks are analysed. Note that for the cells near the riverbanks, the simulated exchange fluxes at the top of the riverbed are a combination of vertical fluxes through the riverbed and horizontal fluxes through the riverbanks. The spatially distributed field for $K_h$ is based on the $K_h$ dataset described in Section 2 [51]. $K_h$-fields are generated for both upstream and downstream sections, using ordinary kriging in SGeMS on a 3D grid with a discretisation of 0.5 m in the horizontal direction and 0.25 m in the vertical direction. At locations outside these sections, average $K_h$ of 6.4 m/day (upstream section) and 5.5 m/day (downstream section) are assigned. $K_v$ is assumed to be 1/10 of $K_h$ [85].

4. Results and Discussion

4.1. Vertical Exchange Flux Estimates Based on 1D Analytical Solution

4.1.1. Temperature Measurements with T-Lance

The average recorded temperature at a depth of 0.1 m (upper boundary) in the riverbed is 16.7 °C for the downstream section and 16.8 °C for the upstream section. Temperatures at this depth show a strong spatial variability (Figure 3 and Figure 4a,b), ranging from 11.7 °C to 18.4 °C and from 12.2 °C to 19.3 °C, respectively. At the deepest measurement depth, the average measured temperature is 12.8 °C at the downstream section and 12.6 °C at the upstream section. The temperature gradient between 0.1 m and the deepest measurement ranges between 0.9 °C/m and 9.4 °C/m, with an average of 4.8 °C/m.

All recorded temperature profiles show a decrease in temperature with depth, which is consistent with higher temperatures of surface water compared to groundwater during summer (Figure 3). Temperature profiles for both sections are relatively similar and all profiles are convex upward, indicating an upward flux towards the river at the time of measurement. Stronger curvatures correspond to higher discharges. The Root Mean Square Error (RMSE) of observed versus simulated temperatures for the upstream and downstream sections are 0.216 °C and 0.253 °C, respectively.

4.1.2. Point-In-Space 1D-Flux Estimates

All estimated vertical fluxes are negative, indicating gaining conditions in the river. The average flux is −72.5 mm/day for the downstream section and −69.4 mm/day for the upstream section. The flux estimates range from −32.0 mm/day to −266.9 mm/day for the downstream section and from −28.9 mm/day to −164.6 mm/day for the upstream section. In general, estimates at both sections are in the same order of magnitude. The highest fluxes at the downstream section are observed for the temperature profiles with the lowest temperatures near the riverbed top (Figure 4). Average vertical flux estimates relate well to previous studies for a different section of the Aa River (−65 mm/day in Anibas et al. [21]) and the Slootbeek (−92 mm/day in Anibas et al. [24]), a small tributary near the upstream section, where flux estimates were based on temperature–time-series analysis.

Point maps of measured temperatures at a depth of 0.1 m (Figure 4a,b) and of estimated vertical flux (Figure 4c,d) show that both display a clear spatial variability, even on the meter-scale. The observed spatial patterns are different for both sections. For the downstream section, the temperature near the banks is low, while it is higher in the center of the river. For the upstream section, the pattern is more complex. Low-temperature zones are present to the northeast (NE) and southwest (SW) of the section, near the riverbanks. However, in the northwest (NW), temperatures near the bank are clearly higher. Similar trends can be identified when looking at the vertical exchange flux (Figure 4c,d). Areas showing low temperatures at the upper boundary correspond to zones of high vertical flux. In the upstream section, two zones of higher flux are identified in the NE and SW. In the downstream section, the center of the river is characterised by low fluxes, while high fluxes are estimated near the banks.

4.1.3. Spatial Variability and Anisotropy

To quantify the spatial variability and directional anisotropy of the flux, omnidirectional and directional (along and across river flow) variograms were calculated and modelled for both the downstream and upstream section (Section 3.3 and Figure 5a–d). An overview of the obtained variogram parameters is provided in

Table 2. Summary of parameters of the variogram models of estimated exchange flux for the upstream and downstream sections.

Section	Direction	Nugget (mm²/day²)	Sill (mm²/day²)	Range (m)
Downstream	Omni	0	1900	8
	Along	0	1600	8
	Across	0	2200	8
Upstream	Omni	0	1100	8
	Along	0	1100	11
	Across	0	1100	6.5

. The omnidirectional variograms for both sections are fitted with spherical variogram models, however, nugget and range still vary (Table 2). The directional variograms along and across the river for the downstream section (Figure 5b) are modelled with zonal anisotropy, i.e., the range in both directions is the same (and equal to the range of the omnidirectional model), but the magnitude of the sill is different for each direction. For the upstream section, the directional variograms (Figure 5d) have the same nugget and sill as the omnidirectional variogram but have a different range, i.e., the anisotropy is geometric.

The ranges of the variograms presented here (Table 2) are in the same order of magnitude as those for horizontal and vertical riverbed hydraulic conductivity calculated by Ghysels et al. [51] for the same sections. Variograms of chargeability and normalised chargeability, which are measured at the downstream section [63], show similar ranges, while electrical resistivity shows a larger range. Variogram ranges along the river are larger than those across the river for both the exchange fluxes and riverbed hydraulic conductivity, which suggests a similar spatial variability and directionality in these riverbed parameters.

Based on the modelled directional variograms, the point-in-space flux estimates are interpolated using ordinary kriging on a grid with a resolution of 0.05 m by 0.05 m, resulting in 2D maps of exchange flux for both the upstream and downstream sections (Figure 4e,f). These maps clearly show distinct flux patterns for both sections. At the downstream section, a clear trend of higher fluxes near the banks and low fluxes towards the center of the river can be identified. Storey et al. [86] and Anibas et al. [21,61] made similar observations, with fluxes being about twice as high at the sides of their investigated river reaches than in the center. For the upstream section, elongated zones of high and low flux can be identified. Near the east bank, fluxes are higher all over as compared to the stream center, while near the west bank, fluxes are high in the downstream but low in the upstream part. In general, colder areas show increased exfiltration as compared to areas with higher streambed top temperatures (compare with Figure 4a,b).

4.1.4. Correlation to Riverbed Hydraulic Conductivity

Studies have shown that riverbed $K$ is often higher in the center of a river than near its banks [50,87]. This is also the case for the Aa River, as demonstrated by $K$ measurements for our site in Ghysels et al. [51]. This effect is especially pronounced at the upstream section, where a low-permeable organic layer at the top of the riverbed is well-developed near the banks. Based solely on these $K$ measurements, one could expect a higher vertical flux in the center of the river as well. However, this does not match with the temperature-based flux estimates discussed in the previous sections.

Point-in-space 1D flux estimates (dots in Figure 4) were compared to horizontal and vertical hydraulic conductivity values obtained by Ghysels et al. [51] for the same spots (Figure 6). When looking at the center of the river, flux estimates for both sections do not vary greatly and therefore no correlation can be observed (Figure 6a–c). There is no correlation between exchange flux and vertical $K$ at the upstream section (Figure 6c). For the downstream section, insufficient data were available to explore the correlation.

For the points near the banks, however, weak linear correlations are found. Vertical and horizontal $K$ at the upstream section and $K_{bank}$ at the downstream show a positive linear correlation with increasing magnitude of the flux, with R² ranging from 0.16 to 0.29 (Figure 6d–f). Larger $K_{bank}$ and riverbed $K$ near the banks result in larger fluxes from the groundwater towards the river.

Our results indicate that, in this study, riverbed $K$ is not the main parameter influencing the magnitude of the exchange fluxes. On the contrary, we hypothesise that the properties of the banks and the adjacent riverbed play a much more pronounced role. Ghysels et al. [58] used groundwater modelling to show that lateral flows through river banks could be a major contributor to total river–aquifer exchange. These observations point to the bank areas as important factors in more accurate quantification of river–aquifer exchange fluxes.

4.2. Fluxes from Groundwater Flow Modelling

Unlike the analytical solution used to quantify vertical fluxes in Section 4.1, the 3D groundwater model for each section allowed us to specifically study lateral flows through the riverbed and banks and estimate detailed exchange fluxes in both sections. For cells near the riverbanks, simulated exchange fluxes at the top of the riverbed are a combination of the vertical fluxes through the riverbed and the horizontal fluxes through the riverbanks.

4.2.1. Downstream Section

Simulated total exchange fluxes vary between −34.8 mm/day and −1973.8 mm/day, with an average of −133.7 mm/day. Vertical fluxes through the riverbed vary between −34.8 mm/day and −524.3 mm/day, with an average of −89.7 mm/day. Lateral riverbank fluxes vary between −463.3 mm/day and −1582.2 mm/day, with an average of −872.0 mm/day, and are significantly higher than vertical riverbed fluxes. The spatial distribution of simulated total exchange fluxes is shown in Figure 7. Vertical fluxes near the banks are significantly higher than those in the center of the river, with fluxes of more than −100 mm/day only observed near the bank areas. In the center of the river, patterns are visible that correspond to the patterns in the riverbed $K$fields of Ghysels et al. [51]. This indicates that in the center of the river the exchange flux is correlated with riverbed $K$.

In Figure 8, flow vectors for a cross section through the river are shown, indicating groundwater flow direction and magnitude in the uppermost part of the aquifer and in the riverbed. Furthermore, simulated riverbed $K$, riverbank $K,$ and aquifer $K$ are visualised. Note that the largest flows are from the aquifer through the riverbanks to the river. Moreover, lateral flow near the banks results in higher vertical fluxes in these areas, explaining the high flux estimates near the banks in Figure 7. In the center of the river, the flow through the riverbed is mainly vertical, and its magnitude is limited compared to the flows near the banks.

4.2.2. Upstream Section

At the upstream section, a well-developed organic layer on top of the riverbed is present near the banks. The $K_h$ measurements of Ghysels et al. [51] were performed at depths of 0.25 m or greater in the riverbed. Therefore, the effect of this low $K$ organic layer is not adequately represented in this dataset and in the simulated riverbed $K$ fields. However, measurements of $K_v$ are present in this section and clearly show that $K_v$ is lower near the banks than in the center of the river. The high $K_v$ in the center can be explained by the scouring of the organic layer during peak flows when the upstream weir is lowered. Based on the $K_v$ measurements of Ghysels et al. [51], the left bank zone has a $K_{v,avg}$ of 0.07 m/day and the right bank zone a $K_{v,avg}$ of 0.14 m/day. These averages are used to change the $K$ of the uppermost riverbed layer at these bank zones. In the center of the river, the $K$ estimates, which are based on the $K_h$ measurement, are kept.

The resulting total exchange fluxes vary between −6.3 mm/day and −2339.2 mm/day with an average of −156.9 mm/day. Vertical riverbed fluxes vary between −6.3 mm/day and −207.5 mm/day, with an average of −68.3 mm/day. Again, lateral bank fluxes are significantly higher than vertical riverbed flux, with the former varying between −498.6 mm/day and −2284.3 mm/day, with an average of −1249.0 mm/day. The spatial distribution of the simulated fluxes is shown in Figure 9. Similar to the downstream section, the riverbank cells are characterised by very high fluxes. However, the effect of the organic layer with low permeability is clearly visible, as relatively low vertical fluxes occur near the bank areas. In the center of the river, the fluxes show a pattern similar to the riverbed $K$ fields of Ghysels et al. [51], and the effect of the elongated high $K$ zone can be clearly distinguished. Note that very low fluxes are simulated in the northernmost part near the center and towards the east bank. Due to the relatively low number of measurements in this area, the quality of the riverbed $K$ field, and thus the simulated exchange flux, is less reliable.

Similar to the downstream section, the largest flows occur from the aquifer through the riverbanks to the river (Figure 10). However, contrary to the downstream section, high vertical flows are present in the center of the river, through the elongated high $K$ zone. Near the banks, vertical exchange fluxes are relatively low due to the presence of the low permeable organic layer at the top of the riverbed. Note the strong horizontal flow components in the riverbed near the west bank. Groundwater flows laterally underneath the low $K$ organic layer and discharges in the center of the river. As such, Figure 10 indicates the importance of non-vertical flow components in this section, near the banks and in large parts of the riverbed.

It should be noted that both sections are modelled using a 3D forward modelling approach assuming steady-state conditions. The quality and reliability of the 3D modelling results could potentially be improved using an inverse modelling approach and calibrating the model on additional data such as high-resolution hydraulic heads when available.

4.3. Comparison of Flux Estimates from Heat Transport and Groundwater Flow Modelling

Summary statistics for the vertical exchange fluxes at the measurement locations (Figure 2) estimated with both the 1D analytical solution and the MODFLOW models are shown in

Table 3. Summary statistics of estimated vertical exchange fluxes based on heat transport modelling and groundwater flow modelling at measurement locations of temperature profiles for both the upstream and downstream sections. Abbreviations: Avg. = Average; Min. = Minimum; Max. = Maximum; St. Dev. = Standard Deviation.

Section	Method	Avg. (mm/day)	Median (mm/day)	Min. (mm/day)	Max. (mm/day)	St. Dev. (mm/day)
Downstream	Heat Transport	−72.5	−57.9	−32.0	−266.9	41.4
	MODFLOW	−68.8	−64.9	−36.7	−206.5	31.2
Upstream	Heat Transport	−69.4	−66.7	−28.9	−164.6	30.9
	MODFLOW	−76.7	−58.3	−11.1	−205.5	54.1

. Boxplots of vertical exchange fluxes are shown in Figure 11. Average values and ranges are similar for both sections. However, the point-estimates for both methods do not always correspond well, as can be seen from the scatterplots in Figure 12. This is especially the case for the upstream section, where the difference between the two methods is relatively large ($R=-0.15$; Figure 12b). For the downstream section, the agreement between both methods is somewhat better ($R=0.16$; Figure 12a).

Based on the spatially distributed maps of exchange flux for both the upstream and downstream sections and both methods (Figure 4e,f, Figure 7, and Figure 9), the total exchange flux averaged over the respective area of these sections can be calculated. For the downstream section, the total flux based on the kriging map of 1D flux estimates using the heat transport model is −74.4 mm/day (Figure 4e) while the flux simulated with the groundwater model is −133.7 mm/day (Figure 7). Expressed in volumetric fluxes, this corresponds to −40.9 m³/day and −73.5 m³/day, respectively (Figure 13). For the fluxes simulated with the groundwater model, vertical fluxes through the riverbed contribute −49.3 m³/day, while lateral flows through the riverbanks account for −24.2 m³/day or roughly one-third of the total flux (Figure 13). For the upstream section, the respective fluxes are −74.4 mm/day and −156.9 mm/day (Figure 4f and Figure 7). Expressed in volumetric fluxes, this corresponds to −21.8 m³/day and −45.9 m³/day, respectively (Figure 13). For the fluxes simulated with the groundwater model, vertical fluxes through the riverbed contribute −20.0 m³/day, while lateral flows through the riverbanks account for −25.9 m³/day or 56% of the total flux (Figure 13).

When looking at the total exchange flux, the groundwater model estimates are about twice as high for both sections as those obtained using the 1D heat transport equation. Schneidewind et al. [88] quantified non-vertical flow components for a section of the Slootbeek, a small tributary to the Aa River. Their results also indicated that non-vertical flow components are in the same order of magnitude, but often larger than vertical flow components and that lateral fluxes through the riverbanks have a strong contribution to total river–aquifer exchange. Hence, estimates solely based on 1D heat transport can potentially lead to an underestimation of exchange fluxes. For our study site, comparing solely vertical flux estimates, we obtain similar results with both methods although spatial patterns are significantly different. As the contribution of lateral fluxes through the riverbanks is added, total exchange flux increases significantly for both sections (Figure 13).

Studies have shown that neglecting or misrepresenting the degree of riverbed heterogeneity [42], anisotropy [49], and the form of the flow field [43,45] will often result in significant errors in vertical flux estimates, depending on the overall flow conditions. For the Aa River, we only encountered fully vertical flow at the center of the river, while closer to the riverbanks non-vertical flow components increased in importance. A similar observation was made by Shanafield et al. [47] in their numerical modelling study of a synthetic channel.

4.4. Limitations

4.4.1. 1D Analytical Heat Tracer Method

The application of 1D heat tracing techniques is easy, fast, and relatively inexpensive. Roaming surveys result in good spatial coverage of riverbed temperature profiles in a relatively short time. One-dimensional techniques commonly assume that the vertical temperature distribution in the riverbed is a result of vertical advective heat transport by the flow of water and by heat conduction through the sediments. However, advective heat transport, due to lateral water flow through the riverbanks and in the riverbed, also affects the temperature distribution in the riverbed. For example, the lateral inflow of relatively cold groundwater in summer results in a decrease in temperature near the top of the riverbed, changing the shape of the vertical riverbed temperature profiles. In a strictly vertical 1D approach, this results in an overestimation of fluxes. This effect might explain the larger vertical flux estimates near the banks in this and other studies [21]. There is an ongoing scientific debate about whether neglecting non-vertical flow components leads to errors in flux estimates [44,46,48]. Irvine et al. [49] found that errors in flux estimates increase as the horizontal flow component increases relative to the vertical component. Roshan et al. [48] stated that, under horizontal flow-through conditions in a riverbed, errors of the analytical 1D method are strong for low vertical velocities. In that case, the analytical 1D method significantly overestimates the vertical velocity. For higher vertical velocities, good estimates are obtained, with higher errors at the gaining side of the stream and with increasing depth in the streambed.

Lateral thermal diffusion might contribute to the attenuation of temperature signals, resulting in wrongly estimated fluxes at both the high and low end of the range. This might also result in the alteration of temperature in the areas surrounding these locations, and thus, the observed temperatures are not solely the result of vertical flow [9].

The assumption of constant thermal parameters over an entire stream section can potentially lead to errors in local vertical flux estimates [89,90]. For example, at two Danish study sites, Sebok and Müller [91] found that thermal conductivity can be spatially variable and that the presence of organic matter and the spatial variability of thermal conductivity can influence the spatial distribution of flux estimates. Aside from conceptual and model structure limitations, uncertainty in vertical flux estimates can also arise from limitations regarding temperature sensor accuracy and placement [89] and measurement resolution [92]. The time of day when the temperature measurements are taken can also have an effect on flux estimates, as discussed by Irvine et al. [75]. To overcome some of these limitations, active heat tracing techniques are being developed that allow for the injection of artificial heat pulses into the riverbed and the quantification of the resulting thermal plume in multiple directions [52,93,94] However, these techniques are more expensive and much more complex in data analysis.

Vertical flux estimates are also potentially impacted by the use of the Bredehoeft and Papadopulos solution [31] itself. Jensen and Engesgaard [22] compared the Bredehoeft and Papadopulos solution to a transient solution and found mean flux estimates from both solution types to be in good agreement at their study site over the course of their study period. Schornberg et al. [42] showed in simulations that the error on the flux estimate due to the assumption of thermal steady-state conditions at the riverbed top (upper boundary) can become significant, especially in low flow zones due to the influence of the diurnal temperature signal. They showed that thermal steady-state condition is a valid assumption essentially only for upwelling conditions, high-temperature gradients between surface water and groundwater (often encountered only in summer or winter), and flux estimates of >100 mm/day. To reduce the impact of the diurnal temperature signal on the estimation error, it could thus be necessary to remove temperature measurements from the upper 5–10 cm of the riverbed.

4.4.2. 3D Numerical Groundwater Flow Model

The advantage of numerical groundwater flow models is that they are physically based, mechanistic simulations, with which river–aquifer exchange fluxes can be easily simulated in three dimensions. Moreover, they allow for the analysis of lateral exchange fluxes through the banks. However, as such, they are limited by their need for spatially distributed fields of riverbed properties as input for more accurate flux estimates [59]. Often, high-resolution riverbed $K$ data are not available because measuring riverbed conductivity is labor-intensive and time-consuming. Even if high-resolution riverbed $K$ data are available, interpolation techniques are needed to generate 3D fields of riverbed $K$. However, these techniques often cannot accurately represent the complexity and connectivity of riverbed structures [95]. Combining point-in-space riverbed $K$ measurements with methods that result in continuous fields (e.g., geoelectrical methods [63]) can improve the characterisation of complex riverbeds.

The reliability of groundwater flow models is strongly influenced by different sources of uncertainty, including the uncertainty on input data, model parameters, and model structure [96]. Because these uncertainties also influence 1D estimates, we would like to again highlight that the quality of any model output is strongly dependent on the quality of the input data and the conceptual model assumptions. As such, a thorough analysis of all sources of uncertainty can increase the reliability of the model predictions.

Riverbank conductivity ($K_{bank}$) is an input to the model that has an important influence on the modelled exchange fluxes in our study. Our estimate of $K_{bank}$ is based on a limited number of measurements (n = 16) at one of the studied sections. Future analysis would contribute well by focusing on the effect of spatially variable $K_{bank}$ on river–aquifer exchange fluxes to increase the reliability of the model results.

Flux estimates can potentially be improved through model calibration on measured hydraulic head or temperature data (inverse modelling). However, to fully appreciate the use of spatially distributed high-resolution riverbed $K$ data, these data need to be available at a similar spatial resolution, which is often not feasible.

5. Conclusions

In this paper, we estimated river–aquifer exchange fluxes both with an analytical solution to the 1D heat transport equation and by applying 3D groundwater flow modelling for two sections of the Aa River, Belgium. All vertical flux estimates based on the 1D analytical solution are negative, indicating upward flow towards the river, i.e., the river is gaining over the observed period. Vertical flux estimates range from −28.9 mm/day to −266.9 mm/day with an average of −71.0 mm/day. The spatial variability of these fluxes on the meter-scale was quantified using variograms. Both studied sections show significant spatial variability with variogram ranges in the same order of magnitude as those of riverbed $K$ and induced polarisation parameters for the same river sections obtained in previous studies [51,63]. Spatially distributed fields of flux were generated using ordinary kriging. No clear correlation between riverbed $K$ and vertical flux estimates was found, except for some minor correlations with riverbed $K$ near the banks and riverbank $K$. Here, we mostly observed higher fluxes near the banks, while other studies (e.g., [50,87]) observed the opposite for riverbed $K$. These results indicate that, in our study, riverbed $K$ is not the main parameter influencing the magnitude of river–aquifer exchange fluxes but rather the properties of the banks and the adjacent riverbed.

Additionally, exchange fluxes were simulated using 3D groundwater flow models (MODFLOW) for both sections, following a modified approach put forward by Ghysels et al. [58]. Fields of spatially distributed riverbed $K$ based on an extensive dataset of $K$ measurements were used to model the riverbed. Simulated fluxes are high near the banks due to the contribution of lateral flows through the riverbanks. Fluxes near the banks are about one order of magnitude larger than vertical fluxes in the center of the river, where patterns of simulated flux reflect patterns in riverbed $K$. Lateral fluxes through the riverbanks account for 33% and 56% of total river–aquifer exchange fluxes for the downstream and upstream sections, respectively. Total simulated fluxes based on groundwater modelling are up to twice as high as those based on the 1D analytical solution, while both methods estimate vertical fluxes in the same order of magnitude. However, co-located point estimates of both methods do not show a clear correlation.

Although both methods used in this study have their advantages and limitations, a robust comparison of flux estimates from 3D groundwater flow modelling to those obtained from applying the 1D heat transport equation is possible, which leads us to conclude that for both river sections, significant non-vertical flow components are present. Non-vertical flow increases considerably towards the riverbanks, and quasi-vertical flow is only found near the center of the river. Our study highlights that care should be taken when inferring on total riverbed exchange fluxes solely based on 1D solutions because these do not provide estimates of horizontal or lateral flows and will likely, therefore, underestimate the total exchange flows.

Riverbed heterogeneity and anisotropy further add to the complexity of groundwater–surface water exchange patterns. For detailed simulations of horizontal and vertical flows in the riverbed, we took advantage of an extensive dataset of riverbed $K$ collected by Ghysels et al. [51], with which we generated detailed spatially distributed fields of riverbed K as input to our 3D groundwater model. However, we are aware that obtaining such detailed information on the spatial distribution of riverbed K in the field is time-consuming and goes beyond the scope of most studies.

More research is necessary to measure non-vertical flows in riverbeds and riverbanks to improve our understanding of riverbed processes and the validity of flux estimates for a variety of field settings. An integrated 3D heat and groundwater flow approach is advised to assess under which conditions a 1D analytical approach is justified or when a more elaborate 3D approach is needed.