mHM is a spatially explicit hydrological model (Samaniego et al., 2010; Samaniego et al., 2021). It considers a wide range of surface processes such as canopy interception, snow accumulation and melting, soil moisture dynamics, infiltration and surface runoff, evapotranspiration, baseflow, discharge attenuation, and flow routing (Samaniego et al., 2010; Kumar et al., 2013). mHM has been successfully applied to several catchments across the globe and has been demonstrated to outperform many hydrological models. It has been applied to basins ranging in size from 4 to 550,000 km² and at variable spatial resolutions (or grid size) (between 1 km and 100 km) (see www.ufz.de/mhm for more details).

In this study, mHM delivers the groundwater recharge. mHM is established following Zink et al. (2017). The simulation period is the same as the one which is used for the groundwater model (1/1/1991—31/12/2013) and daily time steps are chosen. Streamflow from four river gauging stations is implemented to define the inflow through four boundaries where the model perimeter does not agree with the basin watershed (inlet gauges in Figure 2). The quality of this model is supported by the good fitting between the modelled and the observed river flow (Figure 3A) at the outlet gauging station (outlet gauge in Figure 2). The water flow percolating below the unsaturated zone is taken as the recharge for the groundwater model (Figure 3B).

Groundwater flow is simulated using OpenGeoSys - OGS (Kolditz et al., 2012), an open-source scientific software platform for the numerical simulation of thermo-hydro-mechanical/chemical processes in porous media. OGS provides a flexible numerical framework (using primarily the Finite Element Method (FEM)) for solving multifield problems in porous and fractured media for applications in geoscience and hydrology. It has been recently applied in the regional groundwater modeling (Jing et al., 2018; 2019).

The model is divided vertically into two layers whose geometry is defined using information from lithology, geology, and soil characteristics (Figure 4).

The hydrogeological units of the upper layer are defined by using the Hydrogeological Map of Europe (Duscher et al., 2015) and the Global Lithological Map (Hartmann and Moosdorf, 2012). These maps contain information about aquifer productivity, lithology, and degree of consolidation. The hydrogeological units of the lower layer are derived from the Geological Map of Europe (Asch, 2003). Most materials in the upper layer are unconsolidated or partially consolidated and their thickness is defined according to Shangguan et al. (2017). If the thickness (Shangguan et al., 2017) is lower than 10 m, it is assumed that the bedrock (consolidated) outcrops or reaches the top of the saturated zone. At these locations, the materials in the upper layer are also chosen according to the Geological Map of Europe, resulting in the same material for the upper and lower layer. The thickness of the formations that constitute the lower layer is not known, thus, a constant thickness of 500 m is assumed. 500 m is sufficient to minimize the impact of the lower boundary condition (no flow) in shallow water processes. The river network is implemented in the upper layer and its geometry is computed using a digital elevation model. The position and shape of major lakes and springs is obtained from the Hydrogeological Map of Europe (Duscher et al., 2015) (Figure 5).

The hydrogeological model is made up by prism elements (Sachse et al., 2015) whose average horizontal size is approximately 1,000 m and the vertical one is equal to the layer’s thicknesses. The horizontal size is not constant throughout the entire modeling domain. It is refined around wedge-shaped river intersections. Note that the geometry of the river network and the modeled formations is simplified for avoiding too small elements by reducing the number of vertexes that define the rivers. The chosen minimum distance between consecutive vertexes is 1,000 m. Steady and transient state simulations and calibrations are carried out. The transient state model simulates/calibrates from 1991 to 2013 with monthly time steps, while the steady-state model uses averaged data from the same time period.

Three different types of BC are adopted. 1) No-flow BCs are implemented in the lateral boundary assuming that it agrees with the basin’s topographic watershed boundary. There are only four local exceptions in small valleys located in the Alpine region (red dots in Figure 5). In these valleys, in which the lateral boundary does not agree with the watershed, the piezometric head is prescribed according to the topography. No-flow BCs are also assumed at the bottom of the modeling domain. 2) Dirichlet BCs are implemented according to the topography to simulate the river network, the lakes, and the springs. Finally, 3) Neumann BCs are applied at the upper boundary of the model to implement a distributed groundwater recharge, which is previously computed with mHM.

Data from piezometers located at less than 1,000 m from rivers, streams, or lakes are neglected. Otherwise, these observations would have been mapped on the same location as the prescribed head in surface water bodies due to a minimum element size of roughly 1000 m. This would have caused ambiguity and would have forced the calibration to fail. Thus, although more than 900 observation points are available, only data from 200 piezometers are used. Note that some of these 200 piezometers do not contain measurements during the simulation period or the time series is incomplete, which reduces the number of usable observations to 137.

The groundwater model is calibrated by means of the code PEST (Doherty, 1994). PEST is a non-linear parameter estimation software that can be adapted to any simulation code (i.e., it is model-independent) and uses the Gauss-Marquardt-Levenberg method that reduces the discrepancies between computed and measured data to a minimum in the weighted least squares sense (Doherty, 1994) (see the PEST Homepage ² for more details).

Although the measured and simulated piezometric heads fitted well under steady state conditions, similar fittings and acceptable calibration results were obtained when the values of K were modified. This fact does not reflect per se a low sensitivity of the piezometric head with respect to K, it only reveals that the variations of the computed piezometric head at the observation points when K is modified are too small to affect substantially the calibration results of a regional hydrogeological model, in which differences between computed and observed piezometric heads could be considered as acceptable. From theory we know that, the maximum elevation (with respect to the hydraulic head in the rivers) of the piezometric head (h _max [L]) occurs in the groundwater divide between two rivers (Chesnaux, 2013). The equation for calculating h _max that is derived from the general Forchheimer’s solution, is given by (Bear, 1988) as follows: h max = W X 2 8 T (1) where W [L³T⁻¹] is the groundwater recharge, X [L] is the separation between the two rivers and T [L²T⁻¹] is the transmissivity ( T = K b , where b [L] is the aquifer thickness). For Eq. (1), it is evident that h _max increases with small values of K and only slightly vary for high values of K. Given that most piezometers were located in high transmissive materials, their steady-state piezometric head computed with different values of K did not vary enough to affect the calibration statistics (in the context of a regional model). As a result, the computed piezometric heads with different values of K were within the range of acceptable error, and the calibration code could not discern the best values of K. This apparently low sensitivity of the steady-state piezometric head with respect to K in regional models is consistent with previous results from de Graaf et al., 2015, where the sensitivity analyses revealed low variations of the groundwater depth with different K values.

To demonstrate this aspect further, the piezometric heads at the observation points computed with two scenarios with different values of K are compared. The values of K are increased and decreased by one order of magnitude with respect to those obtained from the final calibration (i.e., K values between the two compared scenarios differ by two orders of magnitude). The good fitting between the piezometric head computed at observation points for both scenarios (Figure 6), the high Pearson correlation coefficient R) (>0.99) and the low root mean square error (RMSE) (1.3 m) indicates that the steady-state piezometric head in the observation points is not enough sensitive to K. Thus, calibrated values of K with the steady state model were not reliable.

Although observation points located less than 1,000 m apart from rivers were neglected for the calibration, the used observation points were relatively close to rivers considering the coarse horizontal discretization of the model. For example, only one mesh node may be between a river and observation points located at 2,000 m. Thus, the computed “absolute” elevation of the piezometric head in the observation points is strongly dominated by the river hydraulic head.

If observation points are not located at model nodes, the simulation code assigns them to the nearest node. Thus, in certain cases, observation points are assigned to far nodes when a coarse discretization is used. Near the rivers, where the smallest elements are located, the distance between the actual position of the observation points and the assigned node can reach up to 500 m. Two different situations arise in these cases. On the one hand, the computed piezometric head may differ from the real one because of the hydraulic gradient. For example, assuming a realistic hydraulic gradient in flat areas of 0.001 and a maximum displacement of 500 m, the computed piezometric head would be 0.5 m higher or lower than the measured piezometric head. Thus, there is a minimum error of ±0.5 m encountered, arising between the hydraulic gradient and the displacement of the observation points. This error increases with higher hydraulic gradients than 0.001, which commonly appear in steeper terrains, and far from rivers where the size of the elements gets larger. On the other hand, when observation points are displaced, their elevation is also modified. If an observation point is assigned to a node that is located at an elevation 100 m higher than its actual elevation, it is not possible to fit the modeled piezometric head with the observed one. This effect is most relevant in steep areas because the change in elevation within a single element is higher than in flat areas. For this reason, the error is smaller in flat areas, especially when the shallow piezometers are considered, while it can reach several meters for piezometers located at steeper terrain areas (Figure 7).

The model is vertically divided into three mesh layers, one representing the unconsolidated materials and two for the bedrock. For practical purposes, it is possible to consider that the vertical resolution does not really exist. Therefore, no vertical gradients can be resolved and observations are difficult to be addressed in a model node. This issue is reflected in the fact that errors are higher in deep piezometers than those in shallow ones (Figure 7). Given the relatively small thickness of the upper layer, shallow observation points are less vertically displaced (to be assigned to nodes) minimizing the error. As a result, the model cannot be calibrated by fitting the absolute piezometric head in a steady state. Thus, the calibration should be done under transient modeling conditions. However, the lack of accuracy in computing the exact elevation of the piezometric head also complicates the parametrization under transient conditions.

Two different situations may arise during the calibration process due to the effects associated with the displacement of an observation point: the computed piezometric head can be located above (case 1) or below (case 2) the measured one (Figure 8). If the computed piezometric head is located above, the value of K is increased during the calibration so to decrease the objective function for fitting the average piezometric head, but, the increase of K smooths the groundwater oscillations/fluctuations that cannot be fitted. Contrary, the value of K is decreased to improve the objective function in case 2 (the computed piezometric head is located below the measured one). In this case, the average piezometric head increases but also the magnitude of the oscillations which further complicates the fitting process. In either situation, calibrated parameters by fitting the piezometric head are not representative because they do not allow reproducing correctly the groundwater oscillations. To deal with these issues, we propose an alternative calibration strategy, detailed below.

Considering that the piezometric head increases with low values of K and that groundwater cannot be located above the surface, a steady state calibration assuming the groundwater table is located near the surface (2 m below) leads to the lower bound of K. Although its simplicity, this procedure allows establishing objectively the lower bound of K. The calibration is carried out by fitting the computed piezometric head with the surface elevation at points distributed regularly through the modeled domain. The number of points and their location is chosen arbitrarily. We consider one point every 8,000 m to undertake this calibration stage (Figure 9A). If selected points fall near rivers or streams (less than 1,000 m), they are not considered during the calibration to avoid the influence of the BCs.

The high R-value (0.987), when comparing the computed heads with those derived from topography at the selected points regularly distributed through the modeled domain, indicates that the computed piezometric head is located near the surface (Figure 9B), and thus, that the calibrated value for K can be considered as its lower bound (Table 1).

Information concerning K and the storage coefficient (S [-]) can be also obtained by analyzing base-flow data and application of spectral analyses (Houben et al., 2022). For that purpose, the time series of river baseflow is transferred into the frequency domain for the river flow data (Gelhar, 1974) or piezometric head data (Jiménez-Martínez et al., 2013). This technique is applied to daily data ranging from 1/1/1991 to 31/12/2017 from gauging stations that are strategically chosen according to lithology/geology to obtain information about the different materials (i.e., only observation points located in sub-basins mostly containing one type of material are used) (Figure 10A). Baseflow is computed from river flow data using the digital filter method proposed by Lyne and Hollick (1979). The spectrum of the piezometric head (S _hh ), the recharge (S _RR ) and the baseflow (S _qq ) are related as follows (Zhang and Schilling, 2004): S R R = a 2 1 + t C 2 ω 2 S h h (2)

and, S q q = a 2 S h h (3) where ω [T⁻¹] is the frequency, a [-] is the discharge constant or recession coefficient for a linear reservoir ( a = 3 T / L 2 , where L [L] is a characteristic distance) and t C = S / a [T] is a characteristic time. Assuming that the mean baseflow is equal to the mean groundwater recharge, it is possible to compute t C by fitting the normalized spectrum of the baseflow to one with the following function: f ω = 1 1 + t C 2 ω 2 (4) where f( ω ) is the baseflow normalized function. Computed spectra are automatically fitted to Eq. (4) by using PEST (Figure 10B). Then, the hydraulic diffusivity ( α = T / S [L²T⁻¹]) is obtained from t C taking L as the distance from the river to the watershed (Gelhar, 1974). Here we took the approximation by Dingman (1978) to estimate L from the drainage density distribution as: L = 1 / 2 D (5) where D is the drainage density [L⁻¹] (Horton, 1932).

Information concerning S is needed to calculate the T from α . S is obtained by using data from piezometers located in the same sub-basin of the considered gauging stations, and applying Eq. (3). The ratio between S _qq and S _hh allows calculating a, which is used together with t C to obtain a value for S. Some of the analyzed sub-basins do not have any piezometer that allows computing S. In these sub-basins, T is calculated considering the average S from the other sub-basins. Note that the values of the specific storage coefficient (S _S ) and K, which are required for the setup of the model, are obtained by dividing S and T by the averaged aquifer thickness of every sub-basin. In the case of unconsolidated materials, whose thickness is unknown, a low value of S _S is assumed.

The power spectrum of the selected gauging stations is smoothed by applying a filtering tool to facilitate the fitting (Figure 10B). The power spectra of all considered gauging stations are shown in the Supplementary Material.

Baseflow is an integrative water flux provided by the different geologic formations in the basin, and the flow toward rivers is primarily along high conductivity flow paths. Consequently, values of T (and K) derived from spectral analyses of the baseflow are biased towards higher transmissivity values. Consequently, the results from spectral analyses are assumed to be upper bounds of K (Table 1). In the same manner, the value of S derived from spectral analysis (Table 1) is used to constrain the upper bound of S _S (Jiménez-Martínez et al., 2013). In our case, (6·10⁻³ m⁻¹) is used as the upper bound. Information about the lower bound is not available. Thus, values of 10⁻⁴ m⁻¹ and 10⁻⁵ m⁻¹ are assumed for the unconsolidated and consolidated materials. A lower value is used for the consolidated materials since they are mainly from confined aquifers in the bedrock.

Four alternative calibration techniques are conducted to circumvent the complications emerging from the standard parametrization approach. The first method (M1) consists in fitting the piezometric head anomalies without considering the absolute position of the piezometric head. The difference between the computed and observed piezometric heads is corrected after each iteration taking as reference the initial head difference between them. In this way, difficulties associated with the computation of the absolute position of the piezometric head are avoided. The other three methods are based on fitting statistical parameters of the piezometric head time series. In method 2 (M2) the model is calibrated by fitting the variance of the measured against computed piezometric heads. The third method (M3) matches the Pearson correlation coefficient R of the computed and measured piezometric head at each observation point. Finally, in the fourth method (M4), the model is jointly calibrated using the variance and the Pearson correlation coefficient R.

The fitting between measured and computed piezometric heads with the four calibration methods gave acceptable results at all locations of the modeled basin (Figure 11). Given the number of piezometers used for the calibration, R and RMSE (root mean square error) values are used to evaluate its quality. R and RMSE values are computed for all the observation points (Table 2). They are also calculated with and without applying a bias, accounting for the change in the elevation of the piezometers induced by the numerical model. R and RMSE are firstly calculated considering every observation point without applying a bias correction. In this case, R and RMSE are around 0.95 and 19 m for all calibration methods. Both values are acceptable considering the coarse resolution, the range of simulated piezometric heads (from 309 to 733 m. a.s.l), and especially, the spatial discretization of the model. R and RMSE values are improved substantially when the shift in the elevation of the observation points is eliminated. In this case, R values are higher than 0.99 while RMSE values range between 1.13 and 1.14 m. The indicators (R and RMSE) are so similar among different methods such that it is not possible to discern the best calibration method.

R and RMSE are also calculated individually for each piezometer (Table 2). Slight differences can be observed when comparing the distribution of individual R obtained with each one of the calibration methods (Figure 12). Values of individual R range between 0.45 for the calibration method M1 to 0.5 for the calibration method M3. The highest R values are obtained for the calibration methods M2 and M3. RMSE values oscillate between 0.34 (M1) to 0.32 (M2), which are acceptable considering the average individual variance. The individual analysis of the observation points provides relevant information concerning the best calibration method. Although differences are not too large, the best values of R and RMSE are obtained with the calibration method M2 (i.e., using the variance as reference).

The computed parameters for each one of the calibrations (Table 3) are similar for most of the hydrological units and agree with the materials. There are only some exceptions in which the calibrated value varies among calibrations. Probably, it occurs because, exceptionally, the fitted variable (or variables), which varies for each calibration method, is not sensitive to the calibrated parameter. Even though several challenges arose from the large domain size and coarse discretization, different calibration methods led to very similar and reasonable estimates of hydraulic parameters with small errors. Although results are similar for all the calibration methods, it is important to highlight that, in our case, due to the “master-slave” paradigm used by the parallel PEST (Doherty, 2010; Wang et al., 2019), the computation time decreases when statistical parameters are used because input files are much lighter than that used when all the piezometric head observations are considered.

The steady state results are computed and compared with the measured ones by using the calibrated parameters with the method M2 (Figure 13A). R and RMSE values are >0.98 and 9.9 m, respectively, when only the shallow piezometers are regarded, while R and RMSE are 0.93 and 23 m, respectively, when all the piezometers are considered. The piezometric head distribution considering the parameters derived with method M2 is realistic and, as expected, is controlled by the topography (Figure 13B).