Calibration of groundwater seepage against the spatial distribution of the stream network to assess catchment-scale hydraulic properties

The assessment of effective hydraulic properties at the catchment scale, i.e., hydraulic conductivity (K) and transmissivity (T), is particularly challenging due to the sparse availability of hydrological monitoring systems through stream gauges and boreholes. To overcome this challenge, we propose a calibration methodology which only considers information from a digital elevation model (DEM) and the spatial distribution of the stream network. The methodology is built on the assumption that the groundwater system is the main driver controlling the stream density and extension, where the perennial stream network reflects the intersection of the groundwater table with the topography. Indeed, the groundwater seepage at the surface is primarily controlled by the topography, the aquifer thickness and the dimensionless parameter $K/R$, where R is the average recharge rate. Here, we use a process-based and parsimonious 3D groundwater flow model to calibrate $K/R$ by minimizing the relative distances between the observed and the simulated stream network generated from groundwater seepage zones. By deploying the methodology in 24 selected headwater catchments located in northwestern France, we demonstrate that the method successfully predicts the stream network extent for 80 % of the cases. Results show a high sensitivity of $K/R$ to the extension of the low-order streams and limited impacts of the DEM resolution as long the DEM remains consistent with the stream network observations. By assuming an average recharge rate, we found that effective K values vary between $\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{5}}$ and $\mathrm{1.1}\times {\mathrm{10}}^{-\mathrm{4}}$ m s⁻¹, in agreement with local estimates derived from hydraulic tests and independent calibrated groundwater model. With the emergence of global remote-sensing databases compiling information on high-resolution DEM and stream networks, this approach provides new opportunities to assess hydraulic properties of unconfined aquifers in ungauged basins.

1 Introduction

Evaluating the availability of water resources and its evolution under global changes requires a quantitative assessment of water fluxes at the catchment scale (Fan et al., 2019). Such an evaluation involves the development of advanced hydrological models resolving relevant hillslope- to catchment-scale processes (Refsgaard et al., 2010; Holman et al., 2012; Wada et al., 2010) in a wide variety of high-stake areas (Elshall et al., 2020; Vergnes et al., 2020). Within the local hydrological cycle, aquifers ensure the storage of water during and after recharge periods, increasing the availability of the resources (Fan, 2015; Fan et al., 2015), and transfer this water to surface systems during rain-free periods (Winter, 1999; Sophocleous, 2002; Alley et al., 2002; Anderson et al., 2015). Quantifying groundwater fluxes remains a challenge, as the hydraulic properties of aquifers, i.e., hydraulic conductivity (K) and transmissivity (T), have classically been constrained through sparse borehole-scale characterization (Anderson et al., 2015; Carrera et al., 2005). They are classically estimated using hydraulic tests at centimeter scales for laboratory experiments up to decameter scales for well tests (Renard, 2005; Freeze and Cherry, 1979; Domenico and Schwartz, 1990). Other methods have been proposed at larger scales based on the analysis of streamflow dynamics (Troch et al., 2013; Mendoza et al., 2003; Vannier et al., 2014; Brutsaert and Nieber, 1977), earth tides (Hsieh et al., 1987; Rotzoll and El-Kadi, 2008) and borehole head dynamics (Zlotnik and Zurbuchen, 2003; Jiménez-Martínez et al., 2013) as well as from the calibration of large-scale hydrological models (Eckhardt and Ulbrich, 2003; Etter et al., 2020; Chow et al., 2016). Multi-objective calibration has been proposed to reduce uncertainties, considering complementary data like temperature (Bravo et al., 2002), groundwater ages derived from environmental tracers (Kolbe et al., 2016) or continuous geochemical monitoring (Schilling et al., 2019). In addition, recent advances in machine learning technics show promising results to evaluate hydraulic properties at the regional scale (Cromwell et al., 2021; Marçais and de Dreuzy, 2017; Reichstein et al., 2019).

To tackle the numerous challenges related to the upscaling of hydraulic properties from the local to the regional or global scales, several databases provide exhaustive compilations of measurements performed all around the world (Comunian and Renard, 2009; Achtziger-Zupančič et al., 2017; Ranjram et al., 2015; Kuang and Jiao, 2014). By compiling values obtained from calibrated groundwater models, Gleeson et al. (2014) proposed a global-scale hydraulic conductivity map (GLobal HYdrogeology MaPS, GLHYMPS), with an update by Huscroft et al. (2018), where values have been interpolated based on a high-resolution Global LIthological Map (GliM; Hartmann and Moosdorf, 2012). Besides inconsistencies and methodological biases supported in Gleeson et al. (2014), the compiled permeabilities above the regional scale (>5 km) are not suitable at the catchment scale. Therefore, estimating subsurface hydraulic properties that correctly represent observed catchment-scale processes remains a major challenge for the hydrological community (Blöschl et al., 2019). New opportunities have been identified through the increasing availabilities of surface observations (Beven et al., 2020; Gleeson et al., 2021), specifically with an application for ungauged basins.

Information on the spatial distribution of groundwater seepage appears to be a critical observation to use for the calibration of subsurface hydraulic properties in hydrologic models (Grayson and Blöschl, 2000). This approach can be applied on the assumption that the density and extent of the stream network is primarily controlled by groundwater flow. This assumption is valid in temperate and wet climates for unconfined aquifers, where surface and subsurface hydrological systems are well connected (Cuthbert et al., 2019; Fan et al., 2013) with important discharge of the aquifer directly into the streams (Haitjema and Mitchell-Bruker, 2005). Indeed, a groundwater seepage network dominantly controls the structure of a continuous-stream network, its spatial extent and its ramification (De Vries, 1994; Devauchelle et al., 2012; Strahler, 1964; Leibowitz et al., 2018; Pederson, 2001). Under steady-state conditions, the distribution of groundwater seepage is then controlled by the characteristic hillslope geometry, the recharge rate (R) and the aquifer transmissivity (T), i.e., the product of the hydraulic conductivity (K) and the saturated aquifer thickness (d_sat) (Litwin et al., 2022; Luijendijk, 2021; Bresciani et al., 2014; Haitjema and Mitchell-Bruker, 2005; Gleeson and Manning, 2008). At a given recharge rate, low transmissive aquifers display high groundwater table elevations and, consequently, dense stream networks in the upper part of the catchments. Conversely, highly transmissive aquifers will display lower groundwater tables, higher discharge rates in fewer seepage areas and, consequently, sparser stream networks confined in the lower-elevation valleys (Day, 1980; Lovill et al., 2018; Dunne, 1975; Luo et al., 2016; Dietrich and Dunne, 1993; Godsey and Kirchner, 2014; Prancevic and Kirchner, 2019).

Previous studies have focused on the comparison of an observed hydrographic network with a simulated hydrographic network computed from different methods. The organization of stream networks has been predicted directly from a digital elevation model (DEM) based on a predefined accumulation threshold value determining whether an upstream surface is capable of producing significant flow (Mardhel et al., 2021; Le Moine, 2008; Schneider et al., 2017; Luo and Stepinski, 2008; Lehner et al., 2013). Lumped parameter models, such as TOPMODEL (Beven and Kirkby, 1979), have also been extensively used to predict the spatial patterns of seepage areas (Merot et al., 2003) allowing us to constrain the subsurface hydraulic properties (Blazkova et al., 2002; Güntner et al., 2004; Franks et al., 1998). Luo et al. (2010) proposed a method leading to a spatial distribution of hydraulic conductivities by constraining a 1D groundwater model based on the Dupuit–Forcheimer (DF) assumption on drainage dissection patterns. However, these approaches are limited to considering a subsurface flow path equal and parallel to the downslope topography (Luo and Stepinski, 2008). Relying on explicit simulations of the spatial stream network with a process-based hillslope model following the DF assumption (Weiler and McDonnell, 2004), Stoll and Weiler (2010) proposed to overcome this limitation by routing the downslope subsurface flow from the groundwater table with a grid cell by grid cell approach (Wigmosta and Lettenmaier, 1999). The assessed subsurface hydraulic properties are intended to guide the calibration of hydrological models in ungauged basins, using only a DEM and a hydrographic network map. These approaches are especially relevant as rapid advances in remote sensing are improving the description of global river networks (Schneider et al., 2017; Lehner and Grill, 2013), wetlands (Tootchi et al., 2019; Rapinel et al., 2023) and soil moisture (Vergopolan et al., 2021). Lidar and high-resolution satellite imagery offer new opportunities to determine the surface characteristics of landscapes (Levizzani and Cattani, 2019; Blöschl et al., 2019) and, by extension, the hydrological parameters of local to continental ungauged catchments (Barclay et al., 2020; Dembélé et al., 2020).

Figure 1 Model workflow for the calibration of subsurface hydraulic properties from the observed stream network.

In this work, we propose a new methodology to quantify effective hydraulic properties of unconfined aquifers from topographical and stream network observations now available at high resolution. From a parsimonious 3D groundwater flow model, we aim at estimating the catchment-scale K, based on surface information only, which is the spatial distribution of the stream network. We propose a novel performance criterion to assess the similarity between the simulated seepage areas and the observed stream network, coupled with a stand-alone calibration procedure. We present the full methodology and its sensitivity to different hydrographic network observation products in 24 catchments covering various geological contexts in northwestern France. We finally discuss opportunities and perspectives to systematically characterize aquifers in ungauged catchments from surface observations.

2 Materials and methods

An overview of the method workflow is illustrated in Fig. 1. Each block refers to a specific subsection detailed below (from Sect. 2.1 to 2.5).

1. A digital elevation method (DEM) is used as the top boundary of the groundwater flow model (Sect. 2.1).

2. Three-dimensional groundwater flow is solved in the model domain, and simulated seepage areas are extracted (Sect. 2.2).

3. A selected stream network independent of the DEM is taken as the observed reference (Sect. 2.3).

4. The dimensionless ratio $K/R$ [–] is calibrated to reach the best match between the simulated seepage areas and the extent of the observed stream network (Sect. 2.4).

5. From the optimized $K/R$, the optimal hydraulic conductivity K_optim [L T⁻¹] is deduced by considering the recharge R. The optimal transmissivity T_optim [L² T⁻¹] is obtained considering the average thickness of the saturated aquifer d_sat [L] computed by the model (Sect. 2.5).

2.1 Definition of the model domain based on the analysis of the topography

We first select the digital elevation model (DEM) that will be defined as the upper boundary of the groundwater model. In this study, we use the 75 m grid resolution DEM available at the scale of France. It is generated from photogrammetric restitution and provided by BD ALTI (IGN, 2021). We also explore the impact of different DEM resolutions on the final estimations of $K/R$. We consider two higher-resolution DEMs of 5 and 25 m also provided by BD ALTI. For coarser resolutions, the 25 m DEM was downsampled with the nearest-neighbor option to larger cell sizes, i.e., 100, 200 and 300 m.

Geospatial processing is performed using the software WhiteboxTools available in Python (Lindsay, 2016), labeled WBT with the respective functions quoted in brackets in the following. First, the raw DEM is corrected by filling all depressions and by removing flat areas (WBT.FillDepressions) to ensure continuous flow between grid cells. The vector point shapefile of the outlet is moved to the location coincident with the highest flow accumulation value (WBT.D8FlowAccumulation) within a specified maximum distance taken as twice the DEM resolution (e.g., 150 m for a 75 m resolution DEM) (WBT.SnapPourPoints). A flow direction raster (WBT.D8Pointer) is used to extract the drainage basin (WBT.Watershed).

2.2 Groundwater flow model parameterization

The MODFLOW software suite (Harbaugh, 2005) is used to solve the groundwater flow equation under steady-state conditions for an unconfined aquifer using a three-dimensional finite-difference approach (Harbaugh, 2005; Niswonger et al., 2011). The hydraulic head, h [L], is calculated as a function of the hydraulic conductivity, K [L T⁻¹], and the recharge rate, R [L T⁻¹], applied on top of the water column. At the surface of the model domain, the drain package (DRN) of MODFLOW accounts for the intersection of the groundwater table with the surface and for the induced seepage. Overland flows and surface water re-infiltration are not integrated as they remain marginal in the conditions of temperate climate and low topographical gradients of the studied sites.

We use the FloPy Python package (Bakker et al., 2016) to set and handle simulations. To reduce uncertainties linked to potential flow across topographic boundaries, a buffer zone is added to the topographical catchment boundaries, increasing the modeled domain area by 10 %. The 3D model domain is discretized laterally using the regular mesh of the DEM and vertically into six layers of equal thickness. Convergence tests are performed to ensure the stability of the result independently of numerical discretization. In agreement with field observations undertaken in the region, a homogeneous thickness of the aquifer, d [L], is set to a constant value of 30 m. This thickness represents the typical depth of the interface between the shallow weathered and/or fractured zone with the underlying fresh bedrock (Dewandel et al., 2006; Roques et al., 2016; Mougin et al., 2008; Kolbe et al., 2016). The model assumes a uniform and isotropic hydraulic conductivity.

The recharge is estimated by the land EXternalized SURFace model (SURFEX version 8.1; Le Moigne et al., 2020). For more detailed information, the reader is referred to https://www.umr-cnrm.fr/surfex/ (last access: 1 September 2023). Supplied by meteorological variables, SURFEX computes the energy and water fluxes at the interfaces between soil, vegetation and the atmosphere (Noilhan and Mahfouf, 1996). The groundwater recharge of SURFEX is computed as the proportion of the water mobilized down to the aquifer after infiltrating through the soil column (Vergnes et al., 2020). SURFEX was supplied by the SAFRAN meteorological reanalysis (Vidal et al., 2010; Quintana-Seguí et al., 2008), available over the French metropolitan area at a 64 km² (8×8 km) resolution. Here, we took the steady-state recharge as the long-term average recharge rates computed over the period 1960–2019 and applied them uniformly over the surface of the 3D domain.

2.3 Spatial distribution of the observed stream network

The observed stream network is extracted from the French hydrographic network database BD TOPAGE as a vector format at the scale of 1:10 000 (IGN and OFB, 2019). The main vector file labeled “Cours d'eau” (rivers) of BD TOPAGE represents the majority of perennial sections of the stream network, i.e., filled and/or continuous-flow segments throughout the year. Note that the information classifying perennial or intermittent streams collected in the database is still under development to increase accuracy (Schneider et al., 2017). It has been rasterized at a grid resolution similar to the groundwater flow model in order to facilitate the comparison of the results (WBT.VectorLinesToRaster). Due to the uncertainty in the positioning of the stream network vector with respect to the DEM, an error of the order of 1 pixel is considered (plus or minus 75 m in this case). The influence of this error is analyzed in the results presented in Sect. 3.1 (Fig. 4). We also quantify the impacts of other DEM product resolutions by considering five other hydrographic network products from three different sources: the global-scale database HydroRIVERS (labeled case A), the French database BD TOPAGE (cases B, C and D) and local-scale inventories (cases E and F) performed within the framework of SAGE (Schéma d'Aménagement et de Gestion des Eaux). The HydroRivers product is derived from the processing of the DEM at lower resolution (approximately 500 m at the Equator), while the local inventories are completed by more detailed field observations. More information on these products can be found in Appendix A.

2.4 Calibration criteria between observed and simulated spatial patterns

For each pixel where seepage is simulated by the groundwater flow model, we trace the nearest downslope flow path to the observed stream and compute its distance D_so [L] (WBT.TraceDownslopeFlowpaths) (Fig. 2a). This function of WhiteboxTools uses the topographic structure to compute the path from cells on the surface to the catchment outlet. This procedure converts the initial discontinuous spatial pattern of seepage zones simulated by the groundwater flow model into a continuous-stream network (Fig. 2b).

Figure 2 (a) Definition of the main metrics used for calibration, with $\stackrel{\mathrm{‾}}{D_{\mathrm{so}}}$, the average distances computed from simulated stream pixels (in orange) to the nearest (downslope flow path) observed stream pixels (in blue) and $\stackrel{\mathrm{‾}}{D_{\mathrm{os}}}$ the average distances obtained conversely. (b) Three-dimensional conceptual diagram of the groundwater flow model and of a cross section through the catchment. Continuous streams are generated from pixels where the simulated groundwater table intercepts the topography. By comparison with the observed stream network, some of the simulated streams are correctly estimated (valid, in green), overestimated (excess, in red) or underestimated (missing, in blue).

The distances of the simulated stream network to the observed D_so are calculated and averaged into the criterion labeled $\stackrel{\mathrm{‾}}{D_{\mathrm{so}}}$ (WBT.DownslopeDistanceToStream). High $\stackrel{\mathrm{‾}}{D_{\mathrm{so}}}$ values are characteristic of stream networks extending far away from the observed steams. We also compute the mean distance of the observed to the simulated stream networks following a similar procedure. The distance D_os [L] from each observed stream network pixel to the simulated stream is computed along the steepest downslope path. In the following, we consider its average $\stackrel{\mathrm{‾}}{D_{\mathrm{os}}}$ obtained over all pixels of the observed streams. High $\stackrel{\mathrm{‾}}{D_{\mathrm{os}}}$ values are characteristic of an underdeveloped stream network. The minimum absolute difference between $\stackrel{\mathrm{‾}}{D_{\mathrm{so}}}$ and $\stackrel{\mathrm{‾}}{D_{\mathrm{os}}}$ (Eq. 1), labeled J, is used as the calibration criterion expressing the closest match of the observed and simulated streams or, in other words, the most relevant combination of missing and excess streams (Fig. 2):

$$\begin{array}{}\text{(1)}& J=\left|\stackrel{\mathrm{‾}}{D_{\mathrm{so}}}-\stackrel{\mathrm{‾}}{D_{\mathrm{os}}}\right|.\end{array}$$

$\stackrel{\mathrm{‾}}{D_{\mathrm{so}}}$ and $\stackrel{\mathrm{‾}}{D_{\mathrm{os}}}$ intersect when the calibration criterion J is met. This criterion based on both $\overline{D_{so}}$ and $\stackrel{\mathrm{‾}}{D_{\mathrm{os}}}$ achieves the best equilibrium between over- and underestimations.

At this point, we define the distance D_optim [L] as the average of $\stackrel{\mathrm{‾}}{D_{\mathrm{so}}}$ and $\stackrel{\mathrm{‾}}{D_{\mathrm{os}}}$ (Eq. 2):

$$\begin{array}{}\text{(2)}& D_{\mathrm{optim}}={\displaystyle \frac{\stackrel{\mathrm{‾}}{D_{\mathrm{so}}}+\stackrel{\mathrm{‾}}{D_{\mathrm{os}}}}{\mathrm{2}}}.\end{array}$$

The smaller the value of D_optim, the better the match of the simulated seepage pattern and the observed stream network. D_optim will thus be used as an indicator of the calibration performance. In order to compare cases with different DEM resolution DEM_res [L], D_optim is normalized by the DEM resolution:

$$\begin{array}{}\text{(3)}& r_{\mathrm{optim}}={\displaystyle \frac{D_{\mathrm{optim}}}{{\mathrm{DEM}}_{\mathrm{res}}}}.\end{array}$$

r_optim [–] should remain small to ensure the consistency of the observed and simulated stream networks. It will practically be limited to 2 considering that the mismatch cannot exceed the resolution of 2 pixels:

$$\begin{array}{}\text{(4)}& r_{\mathrm{optim}}\le \mathrm{2}.\end{array}$$

2.5 Estimating the optimal hydraulic conductivities

The model parameter $K/R$ ratio is calibrated by minimizing the objective function defined by Eq. (1), for a given aquifer thickness (d). Optimization is performed by a dichotomy approach (Burden and Faires, 1985). The convergence criterion is reached when $K/R$ varies by less than 1 %. In order to ensure that K estimates are representative of catchment-scale processes driving the spatial distribution of the stream network, independently of the aquifer thickness set in the model, we computed the equivalent normalized transmissivity, $T/R$, by multiplying $K/R$ by the average saturated aquifer thickness (d_sat) computed by the model at the catchment scale (Fig. 2b). In our modeling approach, K and d are input parameters of the model, while T is an output including the computed d_sat. Finally, optimal transmissivity T_optim and hydraulic conductivity K_optim are evaluated assuming the applied average groundwater recharge rate, R, and with known aquifer thickness.

Figure 3 (a) Location of the pilot catchments (Armorican Massif in Brittany and Normandy, northwestern part of France). (b) Zoom on sites along with a simplified map of the main lithological units (1:1 000 000 scale) in sub-panels A, B, C and D.

Table 1 Main landscape characteristics, model input parameters and modeling results, including hydraulic conductivities and transmissivities for the 24 catchments studied.

2.6 Testing the methodology on selected pilot catchments

The approach is deployed in 24 selected catchments located in Brittany and Normandy (France) (Fig. 3), where an oceanic and temperate climate prevails. The average catchment area ranges from 12 to 141 km² with an average of 58 km² (Table 1), which corresponds to an average of 61 800 elements for the domain model discretization. These catchments were selected because of the diversity of their geological and geomorphological settings. Most of them are also subject to extensive research activities for their importance in providing freshwater to the nearby cities (sites 1, 2, 3, 4, 5, 6, 15, 16, 18, 19) or flooding dynamics (sites 20, 21, 22, 23, 24). Some of these sites are also studied in collaboration with local stakeholders in issues related to water quality and river restoration (sites 8, 9, 10, 11, 12, 13, 14, 17, 18 19) or within observatories and research infrastructures (site 7: long-term socio-ecological research (LTSER) “Zone Atelier Armorique (ZAAr)”; sites 17, 18, 19: French network of Critical Zone Observatories (OZCAR) “Ploemeur-Guidel CZO”). None of these catchments present any reservoir or stream obstacle that would significantly alter the stream network. The study sites cover five major lithologies including Brioverian schist (sedimentary rock), Paleozoic sandstone and schist (sedimentary rock), plutonic rocks (mainly granite), mica schist (metamorphic sedimentary rock), and limestone (sedimentary rock). Sites have a homogeneous lithology (1:1 000 000 scale) throughout the catchment except for five sites (sites 7, 8, 17, 18, 19) that present two lithologies.

Figure 4 (a) Two-dimensional map views of distances computed along the steepest slope from simulated stream pixels to the nearest observed ones (D_so) and from the observed stream pixels to the nearest simulated ones (D_os) for the Canut catchment. Results are presented for the lowest, optimal and highest $K/R$ values explored. (b) Average distances $\stackrel{\mathrm{‾}}{D_{\mathrm{so}}}$ and $\stackrel{\mathrm{‾}}{D_{\mathrm{os}}}$ as functions of $K/R$. The shaded areas in gray around the curves correspond to the 75 m uncertainty range equal to the resolution the DEM. The optimal simulation is obtained for $K_{\mathrm{optim}}/R$ at the intersection between the two curves. At this point $\stackrel{\mathrm{‾}}{D_{\mathrm{so}}}$ and $\stackrel{\mathrm{‾}}{D_{\mathrm{os}}}$ are both equal to D_optim, and, in this case, r_optim is close to 2. The estimated optimal hydraulic conductivity K_optim is derived by using the recharge rate provided by SURFEX, and T_optim is obtained considering the mean saturated aquifer thickness (d_sat). To further illustrate the methodology, an animated figure representing 2D map views of the simulated seepage areas as a function of $K/R$ is available in the Supplement with the associated objective minimization function results in Appendix B.

3 Results

3.1 Detailed analysis of model results at a single site

Before presenting the results obtained for the ensemble of pilot sites, we first illustrate the results of the methodology at one specific site, the Canut catchment (Fig. 3, site 6). We provide details on the different steps of the numerical method and assess their performance. The main results of the calibration method are presented in Fig. 4. The dimensionless ratio $K/R$ strongly controls the spatial distribution of the hydraulic head, i.e., the saturated aquifer thickness d_sat, the shape of the groundwater table and its intersection with the surface (Fig. 4a). As $K/R$ increases, the head gradient decreases and progressively disconnects from the surface. This implies that the seepage areas become sparser, mostly organized downstream close to the catchment outlet. Inversely, lower values of $K/R$ tend to expand the seepage patterns along the valleys and depressions towards the head of the catchment. Figure 4a shows the sensitivity of $\stackrel{\mathrm{‾}}{D_{\mathrm{so}}}$ and $\stackrel{\mathrm{‾}}{D_{\mathrm{os}}}$, considering three values of $K/R$. It confirms that the distance from the observed to the simulated stream network $\stackrel{\mathrm{‾}}{D_{\mathrm{os}}}$ increases with $K/R$ and inversely that the distance from the simulated to the observed stream network $\stackrel{\mathrm{‾}}{D_{\mathrm{so}}}$ decreases with $K/R$.

High-order streams are accurately predicted in all three simulations as shown by the green pixels (Fig. 4a). Low-order streams are more sensitive and drive most of the variations in $\stackrel{\mathrm{‾}}{D_{\mathrm{so}}}$ and $\stackrel{\mathrm{‾}}{D_{\mathrm{os}}}$ as shown by the evolving red to green pixels when changing $K/R$. In other words, the calibration is controlled by the spatial extent of the streams from the valleys to the headwaters following the topographic depressions. D_optim is equal to 147 m and remains smaller than twice the resolution of the DEM indicating a close match of the observed streams (Fig. 4b). Using the DEM resolution of 75 m as an indicator of uncertainty, $K/R$ ranges between 1458 and 7082 (shaded area in Fig. 4b), corresponding to a hydraulic conductivity ranging between $\mathrm{1.1}\times {\mathrm{10}}^{-\mathrm{5}}$ and $\mathrm{4.8}\times {\mathrm{10}}^{-\mathrm{5}}$ m s⁻¹ for an estimated average recharge of 215 mm yr⁻¹. In this case, the optimal hydraulic conductivity K_optim is estimated at $\mathrm{1.8}\times {\mathrm{10}}^{-\mathrm{5}}$ m s⁻¹. We deduce the optimal transmissivity T_optim of $\mathrm{5.1}\times {\mathrm{10}}^{-\mathrm{4}}$ m² s⁻¹ considering an average saturated thickness of 28.1 m simulated by the model (Fig. 4b).

We evaluate the impact of the maximum aquifer thickness on T_optim by running the calibration procedure considering five different values of d: 5, 10, 50, 100 and 300 m. We found that the simulated stream network matches the observed one for all thicknesses (d) (Fig. 5a1, A to F). However, we found differences in the estimated T_optim (Fig. 5a2). For cases C, D, E and F, where the maximum aquifer thicknesses are greater than 30 m, the optimal transmissivity T_optim remains constant at around $\mathrm{4.0}\times {\mathrm{10}}^{-\mathrm{4}}$ m² s⁻¹. For cases A and B with smaller thicknesses (<30 m), T_optim reaches much larger values of $\mathrm{4.1}\times {\mathrm{10}}^{-\mathrm{3}}$ and $\mathrm{1.8}\times {\mathrm{10}}^{-\mathrm{3}}$ m² s⁻¹, respectively. Such divergences come from the breakdown of the Dupuit–Forchheimer assumption. Small thicknesses bring the flow lines closer to the surface and widen the seepage areas (Bresciani et al., 2014), effects that must be offset by substantially higher hydraulic conductivities and transmissivities to lower the water table.

Figure 5 Sensitivity analysis of the method for the Canut catchment (site 6) to (a) the DEM resolution, (b) the density and extent of the steam network displayed by different stream network products, and (c) the aquifer thickness. Panels (a1, b1, c1) show the downslope flow path distances of the simulated stream pixels projected onto the observed reference stream network for $K_{\mathrm{optim}}/R$. Panels (a2, b2, c2) show r_optim as a function of T_optim.

We also investigate the sensitivity of the resolutions of both the DEM (Fig. 5b) and the observed stream network (Fig. 5c) to the estimations of T_optim. The resolution of the DEM has only a minor influence on the estimation of optimal transmissivity, while the resolution of the stream network has a major impact. Figure 5b1 shows that the simulated stream network corresponds well to the observed network for the different DEM resolutions. However, the estimated optimal transmissivities vary significantly across the different cases (Fig. 5b2). For cases C, D and E, the estimated T_optim values remain close to each other, ranging from $\mathrm{5.0}\times {\mathrm{10}}^{-\mathrm{5}}$ to $\mathrm{5.6}\times {\mathrm{10}}^{-\mathrm{5}}$ m² s⁻¹. For case F, T_optim reaches a value of $\mathrm{2.9}\times {\mathrm{10}}^{-\mathrm{4}}$ m² s⁻¹, while for A and B, it takes values of $\mathrm{1.9}\times {\mathrm{10}}^{-\mathrm{3}}$ and $\mathrm{9.7}\times {\mathrm{10}}^{-\mathrm{4}}$ m² s⁻¹, respectively. For cases C to F, the r_optim criterion (Eq. 3) remains close to 1 (from 0.5 to 1.8), smaller than the threshold of 2 (Eq. 4). However, for the 5 and 25 m resolutions tested (cases A and B), the distances $\stackrel{\mathrm{‾}}{D_{\mathrm{so}}}$ and $\stackrel{\mathrm{‾}}{D_{\mathrm{os}}}$ are highly sensitive to the mismatch between an increasingly accurate DEM and a coarsely defined stream network. The main factor driving $\stackrel{\mathrm{‾}}{D_{\mathrm{so}}}$ and $\stackrel{\mathrm{‾}}{D_{\mathrm{os}}}$ is no longer the hydraulic conductivity but the mismatch between the DEM and the observed stream network, with r_optim values becoming larger (respectively, 46.5 and 7.7 for the 5 and 25 m resolutions tested). These results emphasize that DEMs with resolutions that are too fine, here 5 and 25 m, cannot be used with the observed stream network selected in this study, at least at the current stage of the methodological development. Resolutions of 75 m and coarser, however, lead to consistent estimations of the hydraulic conductivity confirming the validity of the modeling approach.

We have systematically tested the method using six different stream network products issued by global, national and local databases (Fig. 5c1). These products display important differences in the extent and densities of the stream network coming from their origin and scale of observations (Appendix A). For case A, the global-scale product HydroRIVERS (Lehner et al., 2013) locates rivers away from the topographic valleys of the DEM, with consequently a D_optim value more than 10 times larger than DEM_res. For cases B to F, the criterion r_optim (Eq. 3) remains smaller than 2, the hydrographic network is well captured and the method estimates T_optim consistently (Fig. 5c2). T_optim varies over 1 order of magnitude from 2.0 x 10⁻⁴ to $\mathrm{3.8}\times {\mathrm{10}}^{-\mathrm{3}}$ m² s⁻¹ and logically tends to decrease when the density and extent of the mapped stream network increases. Indeed, for a fixed recharge rate, lower transmissivities raise the groundwater table and broaden the headwater streams (Fig. 5c, E and F). Conversely, higher values of T contract the hydrographic networks with streams located mainly at lower elevations (Fig. 5c1, A and B). The extent of the stream system and the first-order stream locations appear to be highly sensitive to the estimated transmissivity confirming its capacity to inform the hydraulic conductivity.

3.2 Application to the ensemble of catchments

The method has been applied to the 23 other catchments (Fig. 3) with the same DEM resolution of 75 m, the product for the observed reference stream network and an aquifer thickness of 30 m. Both $\stackrel{\mathrm{‾}}{D_{\mathrm{so}}}$ and $\stackrel{\mathrm{‾}}{D_{\mathrm{os}}}$ were found to systematically intersect defining optimal D_optim (Fig. 6) and $K_{\mathrm{optim}}/R$ values (Table 1). For 19 sites, the value of D_optim is less than 2 pixels (Fig. 6), showing good consistency between the simulated and the observed stream network, and $K_{\mathrm{optim}}/R$ varies between 1397 to 19 687. Considering d_sat computed by the model, T_optim values range over 1 order of magnitude from $\mathrm{2.9}\times {\mathrm{10}}^{-\mathrm{4}}$ to $\mathrm{3.2}\times {\mathrm{10}}^{-\mathrm{3}}$ m² s⁻¹ (Fig. 6), resulting in K_optim values between $\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{5}}$ and $\mathrm{1.1}\times {\mathrm{10}}^{-\mathrm{4}}$ m s⁻¹ (Table 1). The model captures correctly the features of the observed stream network even in the presence of singular topographical features such as extended depressions or sharp changes in slope (Gauvain et al., 2021; Schumm et al., 1995). This is especially the case at site 7 where the seepage along the foot slope issued by a steep slope transition (6 % on 1000 m of length) located along a lithological contact is well represented both by the model and in the observations as a significant and perennial groundwater spring and/or wetland (Vautier et al., 2019; Kolbe et al., 2016).

Figure 6 D_optim and r_optim criteria as functions of T_optim estimated for the 24 sites. The optimal transmissivity T_optim is obtained by considering K_optim, and the mean saturated aquifer thickness d_sat is computed. The shaded area corresponds to sites with r_optim>2. The DEM resolution is 75 m and the aquifer thickness is 30 m. The error bars correspond to the estimated T_optim considering the DEM resolution as an uncertainty indicator.

4 Discussion

4.1 A new calibration method for the assessment of effective catchment-scale hydraulic properties

We have presented a process-based groundwater modeling approach to assess effective catchment-scale hydraulic properties (K and T) from the sole information on the density and spatial extent of the stream network. The proposed method (1) simulates the stream network by physically representing the 3D groundwater flows, (2) quantifies the mismatch between observed and simulated stream networks through distance-based spatial indicators and (3) calibrates the subsurface hydraulic properties by minimizing a performance criterion. We compare it with existing approaches and indicators.

Previous modeling approaches were mostly based on TOPMODEL applications (Blazkova et al., 2002; Güntner et al., 2004; Franks et al., 1998) or hillslope-scale flows (Luo et al., 2010). They had the advantage of simplicity but were limited in simple cases where flows are topography-driven. Like the 1D hillslope-scale approach of Stoll and Weiler (2010), our 3D distributed- and process-based groundwater approach remains valid in conditions where the groundwater table is not a strict replicate of the topography becoming “recharge controlled” (Haitjema and Mitchell-Bruker, 2005). It also accounts for the discontinuities of the seepage structure coming from the irregular connections of the subsurface flows with the surface (Godsey and Kirchner, 2014; Whiting and Godsey, 2016; Warix et al., 2021). It follows numerous field observations and synthetic experiments showing the strong influence of 3D groundwater flow organizations on the location of seepage zones (Goderniaux et al., 2013; Fleckenstein et al., 2006; Gauvain et al., 2021; Dohman et al., 2021).

Several indicators have been proposed to compare the spatial patterns of the modeled and observed saturated areas mostly based on cell-by-cell and cell–neighborhood approaches. This includes the likelihood measure (Franks et al., 1998; Blazkova et al., 2002), the Kappa goodness-of-fit statistic (Stoll and Weiler, 2010) and the Euclidean distance between cells (Güntner et al., 2004). These indicators are essentially local and readily accessible from information on topography and stream networks. We propose more integrative indicators based on the distance between the observed and simulated stream networks computed along the steepest slope between them, as does the IDPR (Network Development and Persistence Index), to identify zones predominantly favorable to infiltration or runoff (Mardhel et al., 2021). The advantages of this procedure are to account for the topographical structure within the definition of the distances and to constrain the comparison to the best compromise between the over- and undersaturation, mainly driven by $\stackrel{\mathrm{‾}}{D_{\mathrm{so}}}$ and $\stackrel{\mathrm{‾}}{D_{\mathrm{os}}}$, respectively.

4.2 Comparison of estimated hydraulic conductivities with previously published values

As shown in Fig. 6, the method predicts a distribution of T values that stands within 1 order of magnitude despite the broad range of lithological units investigated. Overall, our estimates of hydraulic properties are consistent with values found in previous studies conducted for similar sites and lithologic settings (Roques et al., 2016; Mougin et al., 2008; Dewandel et al., 2021; Cornette et al., 2022; Leray et al., 2012). We compared the estimated hydraulic conductivity, K_optim, with local to regional values found in the literature K_literature (Fig. 7). We focused on the comparison with local estimates from hydraulic tests or numerical groundwater models for two of the studied sites (Le Borgne et al., 2006; Kolbe et al., 2016; Jiménez-Martínez et al., 2013) or on those compiled in regional syntheses according to the lithology (Laurent et al., 2017; BRGM, 2018). Figure 7 shows a good agreement between our results and the other values. More specifically, the local values extracted for catchments 7 and 19 slightly underestimate the one from the literature by about 33 % (Fig. 7, diamonds). The hydraulic conductivities derived from regional synthesis remain within the same order of magnitude with a limited overestimation of a factor of 2 (Fig. 7, circles). This slight overestimation might result from testing methods as well as the fact that local hydraulic tests are often carried out close to transmissive geological features of major interest for water supply. Our results were also compared to those compiled in the global-scale database GLHYMPS (Huscroft et al., 2018). From the database, we derived equivalent hydraulic conductivities for each of the catchments. They range over 4 orders of magnitude, systematically lower by 1 to several orders of magnitude than our estimates (Appendix C). As shown in previous studies (de Graaf et al., 2020; Tashie et al., 2021), we find that the hydraulic conductivity dataset compiled in GLHYMPS may be locally underestimated.

Figure 7 Comparison of K_optim obtained for the 24 catchments with values found in the literature K_literature grouped into two categories according to the scale of investigation (local vs. regional). Diamonds: values from hydraulic tests and/or groundwater modeling compiled in regional syntheses according to lithologies (BRGM, 2018; Laurent et al., 2017). Circles: K obtained from local hydraulic tests for site 19 (Le Borgne et al., 2006; Jiménez-Martínez et al., 2013) and groundwater modeling for site 7 (Kolbe et al., 2016). The values provided in transmissivity by the literature are translated into hydraulic conductivities by applying the same aquifer thickness of 30 m.

Figure 8 For the five sites corresponding to r_optim>2, representation of the observed stream network on top of the simplified geological map (1:1 000 000 scale), with the downslope flow path distances of the simulated seepage areas projected to the observed streams. For site 8, differences are larger on the plutonic rocks. For the other sites, the white square identifies the area where differences are the largest.

The results also show a range of values consistent with those given by classical textbooks for the investigated lithologies (Freeze and Cherry, 1979; Domenico and Schwartz, 1990). Crystalline rocks characterized by weathering and fractures (Roques et al., 2014; Dewandel et al., 2006) have a lower hydraulic conductivity than sedimentary rocks, here represented by limestones with a karstic systems (IGN and OFB, 2019) that, as expected, display higher conductivities. Lower conductivities suggest a high water table inducing a larger spatial extent of the stream network, confirmed by local knowledge, with a much higher observed drainage density for crystalline sites compared to limestone sites (Table 1). Although our results show evidence that effective conductivities are related to variations in dominant lithologies, it is also clear that other reported factors like erosion, bedrock weathering and fracturing may tend to homogenize the hydraulic properties in similar erosion and/or weathering settings (Luo et al., 2016; Yoshida and Troch, 2016; Jefferson et al., 2010; Litwin et al., 2022).

4.3 Sensitivity to input and/or model parameters: related improvements for broader applicability

The D_optim indicator (Eq. 2) provides information on the level of uncertainty. It depends on the quality of the data (observed reference stream network and the DEM) and on the model assumptions. Five sites (8, 18, 21, 22 and 23) over the 23 sites studied display a low match between the simulated and observed stream networks (r_optim>2), with D_optim values ranging from 190 to 581 m (Fig. 6). Figure 8 maps the simulation results for these five sites. For site 18, the model predicts a seepage zone induced by a topographic depression representing potential ponds, lakes and wetlands, which are not in the observations. For sites 21 and 22, the main differences are located at a non-reported subsurface flow in the observed stream network within a karstic system (IGN and OFB, 2019); however, this is represented well by the model due to a topographic depression along this area. For site 23, the simulated flow at the bottom of the DEM valley appears to be parallel to the observed flow. In this case, it seems to be a registration error in the alignment of the stream location data within the DEM. For site 8, differences come principally from the model and, more specifically, from the assumption of a uniform hydraulic conductivity. For this site with lateral lithologic heterogeneity, we found that the model underestimates the extent and density of the stream network in the part with dominant plutonic rocks and overestimates them in the schists. At site 8, the IDPR (Mardhel et al., 2021) indicates that the granitic area is less permeable than the schist area and generally displays the limestone sites 21 and 22 primarily dominated by infiltration, consistent with our results.

The high sensitivity of applied $K/R$ to both the density and spatial extent of the observed stream network highlights the requirement of using high-quality stream products. Several issues may arise. First, river maps available in national to global databases are often incomplete compared to local databases compiled by stakeholders in direct field observations, leading to an overestimation of the effective hydraulic properties. Second, artificial channels, drainage ditches and any other departure from the geomorphic equilibrium may alter the stream network system and lead to an overestimation of the hydraulic properties in highly impacted zones. Third, the resolution of the DEM and reference stream network must be close. Nevertheless, the observed stream network layer could be adjusted to better match the DEM resolution.

Major limitations and improvements may also arise from the assumptions of the hydrogeological model. The proposed methodology can be used with other parameterizations and model conceptualizations. The model domain can be extended and the boundary conditions modified to better represent potentially longer and deeper regional groundwater circulations. If the information is available, the method could be tested with heterogeneous recharge at the catchment scale. The geometry of the aquifer could be adapted, applying flat- or irregular-bottom aquifers based on geophysical measurements (Pasquet et al., 2022) or depth to bedrock databases (Shangguan et al., 2017; Hengl et al., 2017; Pelletier et al., 2016). Exponentially decreasing conductivity with depth can be applied to the model in order to estimate the initial value of K at the surface or its characteristic decay depth. At the current stage of the method, catchment-scale lithological heterogeneities can be considered by applying the methodology independently on sub-areas characterized by a homogeneous lithology. For example, at the studied site 8, the application of the methodology on granite-dominated sub-catchments should result in lower K estimates than in the schist areas. Localized heterogeneities including weathering, fractures, faults and other discontinuities cannot be identified. They should be explicitly introduced into the model and characterized by other methods.

5 Conclusions

Global syntheses compiling accurate predictions of hydraulic properties of the subsurface are critically needed to predict water resource availabilities (Fan et al., 2019) in ungauged catchments (Sivapalan et al., 2003; Hrachowitz et al., 2013) and to assess the impact of hillslope- and catchment-scale hydrology on global change predictions (Taylor et al., 2013). Besides the climatic forcing data, requiring only a stream network map to calibrate a groundwater flow model built from a DEM, the approach presented in this article addresses this challenge, specifically for ungauged basins. Under the assumption that the transmissivity (hydraulic conductivity integrated over the saturated aquifer thickness) controls the extension and density of the hydrographic network, the approach calibrates the effective hydraulic properties against the stream network. The resolution of the stream network and DEM should be consistent. We showed that the spatialized performance criterion based on the distances between the simulated and observed stream network achieves an equilibrium between over- and undersaturation of the underlying groundwater system, as well as an equilibrium between over- and underestimation of the stream network extent. The resulting estimated K values are consistent with local values found in the literature.

A major advantage is the ease of deployment and transferability of the methodology to other catchments. Although the proposed methodology is limited to unconfined aquifers and particularly suited to contexts with significant subsurface–surface interactions, where groundwater primarily feeds streams, it aims to be deployed at multiple spatial scales by taking advantage of databases compiling topographic, hydrologic and climate information. Such deployment improving subsurface characterization from surface information would leverage the current development in crowdsourcing (Etter et al., 2020) and innovations in remote sensing (Biancamaria et al., 2016) that now provide high-resolution surface DEM products (Hawker et al., 2022; Yamazaki et al., 2017) and mapping of hydrographic networks (Grill et al., 2019; Yamazaki et al., 2019) distinguishing the perennial and intermittent streams (Fovet et al., 2021; Messager et al., 2021).

Appendix A: Hydrographic network product information

A1 Global database

HydroRIVERS (available on this website: https://www.hydrosheds.org/page/hydrorivers, last access: 1 September 2023) (Linke et al., 2019) is derived from HydroSHEDS (Lehner et al., 2013), a mapping product that provides stream information for regional- and global-scale applications, based on a grid resolution of 15 arcsec (approximately 500 m at the Equator).

A2 National database

The BD TOPAGE (available on this website: https://bdtopage.eaufrance.fr, last access: 1 September 2023) database classifies streams as perennial or intermittent based on historical photogrammetric reconstructions.

A3 Local database

The local stream and wetland inventory maps are based on observations and field surveys validated by the “Schéma d'Aménagement et de Gestion des Eaux” (SAGE) (available on this website: https://sdage-sage.eau-loire-bretagne.fr/home.html, last access: 1 September 2023).

Appendix B: Objective function of the calibration criteria

Figure B1 For the Canut catchment (site 6), (a) 2D map views of simulated seepage areas and nearest downslope flow path distances (simulated to observed and observed to simulated) for the two $K/R$ values at the bounds of the uncertainty (lower and higher) and (b) the objective function based on the developed performance criteria obtained.

Appendix C: Comparison with the GLobal HYdrogeology MaPS (GLHYMPS 2.0) permeability

Figure C1 Comparison of K_optim obtained for the 24 catchments with values from the literature, K_GLHYMPS, provided by the GLHYMPS 2.0 global permeability map (Huscroft et al., 2018) and averaged over each catchment.