Assimilation of Piezometric Data to Calibrate Parsimonious Daily Hydrological Models

Abstract

Low water levels are a seasonal phenomenon, which can be long, short, and more or less intense, affecting entire watercourses. This phenomenon has become a concern for many countries who seek better understanding of the processes that affect it and learn how to optimally manage water resources (pumping, irrigation). Consequently, a lumped rainfall model at daily time step (GR) has been defined, calibrated, and regionalised over French territories. The input data come from SAFRAN, the distributed mesoscale atmospheric analysis system, which provides daily solid and liquid precipitation and temperature data throughout the French territory. This model could be improved, in particular to more accurately simulate the hydrological response of watersheds interacting with groundwater. The idea is to use piezometric data from the ADES bank, available in France, and to use it for the calibration phase of the hydrological model. The analysis was carried out across ten French catchments that are representative of various hydrometeorological behaviours and are located in a diverse hydrogeological context. Each catchment must be represented by a piezometer that closely represents the main aquifer that interacts with the basin. This piezometer is located on part of the watershed that is most covered in terms of its drainage network, and closest to its outlet. Different signal processing methods are used to characterise the relationship between the fluctuation of river flow, piezometric levels and rainfall time series. Potential processing methods will be carried out in the temporal domain. To quantify groundwater table inertia and that of the catchment area, correlograms were calculated from daily chronicles of flows and piezometric levels. A cross-correlatory analysis was set up to see, in more detail, the correlations between the flow rates (especially base flows) and piezometric level time series. This type of analysis makes it possible to study relationships between various observations, and tests were carried out to take this information into account during the phase of the calibration of hydrological model parameters. These different analyses will hopefully help us to use piezometric data to consolidate the quality and robustness of the modelling.

1. Introduction

Water resource management is one of the major issues of our time. Many factors, both natural and man-made, influence water resources and, more generally, the water cycle. Groundwater storage through aquifers is a key element of the water cycle, especially during droughts. Aquifers store water during rainy periods and contribute to stream flows during non-rainy periods. Thus, the measurement of water levels in aquifers (through piezometric data) provides information on the water content of soils and their possible availability for low water support in a watershed. Prior to addressing many issues related to water resources (irrigation, drinking water supply, navigation, hydroelectricity, respect of ecosystems), both underground and surface waters, and in the perspective of an integrated management of this resource, understanding how environments interact is essential, to better simulate them. Hydrological modelling generally takes into account different water cycle components, such as precipitation, evaporation, infiltration, and runoff. The most widely available data used by modellers are precipitation (rain gauge networks, weather radars) and flow measurements at the watershed outlet (hydrometric network). Evaporation is generally deduced from formulations combining climatic data, such as temperature, radiation, and wind [1]. Moisture and piezometry data are sometimes considered in physically based modelling [2,3,4] and in hydrogeological models [5]. In contrast, this information is rarely considered in conceptually modelling the rainfall–flow relationship.

Awareness of the importance of interactions between groundwater and rivers appeared in the 1960s [6,7], but started with the study of karsts [8]. Since the 2000s, many studies have been carried out on this subject, particularly with the adoption by the European Parliament of the Water Framework Directive (WFD). Among others, the aim of this law is to assess the sound status of water bodies requiring that exchange dynamics between surface water and groundwater be considered. However, hydrological interactions between surface water and groundwater remain very complex [9] because exchanges take place at various scales of space and time. Indeed, these interactions are very complex processes induced by geomorphology, hydrogeology, and climatic conditions [9].

In the context of hydrological modelling, piezometry is information that can be used to evaluate the different states reflecting a model’s water content at a given time step.

An increasing number of numerical weather prediction models include a representation of groundwater [10,11]. Ruelland’s study [12] also aims to simulate the relationship between climate forcing and the dynamics of groundwater levels and runoff in the upper Elqui watershed (5660 km², Chile) with the conceptual daily model HydroStrahler. Nevertheless, such representations are not detailed enough to be used to monitor or forecast groundwater resources. This is why some dedicated approaches aim to provide groundwater level forecasts at the well scale with lumped models [13].

Physically based distributed models allow for the use of piezometric data, such as calibration or validation data [14]. Specific models have also been developed to only simulate piezometry data, without considering flows [15]. Another study, which covers a wide domain, corresponding to a major part of the US (6.3 billion of km²), was carried out by Maxwell et al. [16]. A three-dimensional hydrogeological model (ParFlow) was used at a 1 km grid resolution in a steady-state run. This model has four layers over the first metre of soil and then a fifth layer from 1- to 100-metre depths. The computation time was one week on a high-performance computer for a steady-state simulation. Thus, while the study confirms the possibility of running a three-dimensional groundwater model at fine resolution over a very large territory, it is still difficult to consider its application for operational water management purposes.

In order to simulate piezometry while modelling stream flows, it is important to define the state in the model structure to which it will be possible to match the piezometry. For example, in the Gardenia model [17,18], the “groundwater” reservoir is identified with the aquifer’s storage capacity and calibrated using piezometric data. The hydrometeorological modelling platform (AquiFR) also represents the main hydrological processes which occur within watersheds from precipitations to groundwater flows. The AquiFR system includes three hydrogeological modelling software programmes: the EauDyssée hydrogeological numerical platform [19], the MARTHE (Modelling Aquifers with Rectangular cells, Transport and Hydrodynamics) groundwater flow software programme [20] and the EROS (set of rivers organised in sub-basins) lumped model software programme used for karstic systems [21].

Most of the developed structures of conceptual hydrological models, therefore, include weaknesses to correctly simulate flows coming from groundwater and resulting exchanges [22,23,24]. Indeed, they have been designed to approximate the hydrological processes which underlie the generation of flows, their only control method being observed flow records. The basin’s behaviour can then be described using interconnected reservoirs, representing the basin’s various recharge, storage, and discharge systems. Several links are then formulated between groundwater flow models and their observation through soil moisture and piezometric levels.

The objective of this study is, therefore, to take into account the piezometric information available in watersheds and integrate it into the calibration phase of a conceptual hydrological model, as a way of understanding low-water issues. The objective is to add an additional control from observations of subsurface behaviour to improve hydrological model calibrations and limit possible equifinality problems [25]. To meet this objective, simulations of the subsurface compartment will be viewed via the available piezometric data in the watershed, without modifying the model structure.

2. Data Presentation

Study Sites

Ten natural watersheds, spread over the territory of metropolitan France, were selected to test the methodology on a diversified panel of watersheds (Figure 1).

These watersheds have a limited anthropic influence. They have a sufficiently long length of observed flow records, at least 25 years of available data, and a low rate of gaps (10% maximum) so that the modelling results are not overly sensitive to sampling problems. In addition, the length of the available piezometer record should match the length of the flow record as closely as possible, with a minimum number of gaps. This sample of basins comprises different basin areas, a diversified hydrological regime, and is located in various geological and hydrogeological contexts. This allowed for watersheds with varying degrees of influence from the subsurface compartment. Differences between different rainfall regimes correspond to relatively low summer water and high winter water, depending on regions and climates (oceanic, temperate, or continental). This map and its layers were produced using the GIS platform, ArcGIS, which is a complete GIS system to fully exploit in 3D geographical data.

To study relationships between water cycle data (rain, surface water, groundwater) and to use them in hydrological modelling, we collected rainfall (noted P), flow (noted Q) and piezometry (noted H) records present across our ten study sites.

❖

Meteorological data were obtained from the SAFRAN reanalysis developed by Meteo France. The SAFRAN reanalysis is based on an optimal interpolation method of meteorological data at the scale of massifs, defined from a division of France into climatologically homogeneous zones. It provided daily liquid and solid precipitation and temperature records, available on a regular 8 km × 8 km grid across France [26,27]. These data are available for the 1970–2020 period. The evapotranspiration data used in modelling were calculated using Oudin’s formula [1], depending only on air temperature.

❖

Hydrometric data were presented as daily average flow records and were taken from the national data bank for hydrometry and hydrology (Hydro French database) [28]. Chronicle durations vary, depending on basins, over a reference period chosen between 1 January 1970 to 31 December 2020. This period was chosen because it is sufficiently broad, with wet years, such as 1977, and dry years, such as 1976, 1989, 2003, and 2011.

❖

Groundwater data were extracted from the ADES (National Groundwater Data Access Portal), a french public site, which collects quantitative (water-level chronicles) and qualitative (geochemistry) data on groundwater in metropolitan France and overseas. It is an essential tool for optimal management of groundwater resources, to enhance our understanding of groundwater evolutions and contribute solutions for local, national pressing societal and European requirements. Values at variable time steps were aggregated at daily time steps so it could be compared with rainfall and flow data.

3. Data Processing

All the data processing presented in this section will be performed with R language, a free software environment for statistical computing and graphics.

3.1. Computation of a Base Flow and Base Flow Index (BFI)

To study the exchange of water with an underground compartment, one must determine the total flow share that contributes to it. Flow record hydrographs have, therefore, been broken down into a slow component, called the base flow (Qb), and a fast component. This second component is used to evaluate the storage capacity of a watershed.

Many numerical methods separate the base flow. Several numerical filters appear to be suitable for this purpose, such as the Institute of Hydrology’s [29] local sliding window minimum, or recursive numerical filtering [30]. Although these different approaches do not accurately reflect the base flow, the Lyne and Hollick digital filter does provide a simple, objective, and reproducible base-flow estimate [30].

The method used follows the same principle as the recursive filter, assuming that precipitation has a non-negligible influence on base flow [31]. The filter equation (Equation (1)) is then governed by a single, so-called filter parameter (β), which was determined by Nathan and McMahon and takes into account the flow on the day under study and the previous day.

The authors also hypothesised that running the function once in the chronic direction, a second time in the reverse direction, and a final time in the normal direction, results in a more responsive base flow (Figure 2) that relates less to total flow during recession periods.

The equation for this method is as follows:

$$\begin{gathered} Q_t=\beta \times Q_{t-1}+\frac{\left(1+\beta \right)}{2}\times \left(Q_t-Q_{t-1}\right) \\ with Q_t as the fast flow at time t ; Q_{t-1} as the previous flow, \\ and \beta as the fixed parameter. \end{gathered}$$

(1)

As seen in the Nathan and McMahon study [31], a value of β = 0.925 is retained.

To illustrate this treatment, Figure 2 presents the hydrograph, separated according to the recursive digital filtering method for the Aa basin for the year 2008, as well as daily rainfall, total discharge and piezometry chronicles.

The ratio between average base flow (Qb) and average total flow (Q) is calculated and called the Base Flow Index (BFI). This index was developed in the United Kingdom [32] as a way of classifying soils according to their hydrological response in both low-flow and flood studies. The formula for calculating this index (Equation (2)) is as follows:

$$\begin{gathered} BFI=\frac{{{\displaystyle \sum }}_n^{i=1}{Q}_{bi}}{{{\displaystyle \sum }}_n^{i=1}{Q}_i} \\ with Q_b, the base flow determined by the previous separation method and Q \\ the total flow. \end{gathered}$$

(2)

A criterion of permeability was also extracted from the French hydrogeological referential of the BDLISA database. The hydrogeological repository BDLISA (Aquifer System Boundary Database) is a national database that allows for locating data related to groundwater. It aims to provide a mapping of hydrogeological entities on the whole territory. It has been developed by the BRGM (Geological and Mining Research Office), the French public establishment of reference in the applications of Earth sciences to manage the resources and risks of the soil and subsoil. BDLISA is used to characterise the watersheds. Watersheds, with a high permeability, are subject to infiltration, i.e., water is able to flow towards the deeper watershed layers, and thus potentially feed the water table.

BFI values and the permeability index in

Table 1. Different geographical and hydroclimatic characteristics determined at study watersheds.

Code	BV1 - A3792010	BV2 - E3511220	BV3 - E4035710	BV4 - E5400310	BV5 - H5172010	BV6 - H6221010	BV7 - H7162010	BV8 - J4313010	BV9 - J4742010	BV10 - J5102210
Watercourse	Sauer	Lys	Aa	Canche	Saulx	Aisne	Serre	Steir	Ellé	Scorff
Outlet	Beinheim	Lugy	Wizeres	Brimeux	Vitry-en-Pertrois	Givry	Pont à Bucy	Guengat	Arzano	Plouay
Area (km2)	502.2	86.7	391.6	917.2	2109.1	2887.6	1637.4	179.9	575.7	299.3
Average annual precipitation (mm/yr) (PA)	808	1017	975	948	983	897	807	1191	1130	1126
Average annual ETP (mm/yr) (ETPA)	676	614	611	622	650	654	647	668	655	662
Average annual flow (mm/yr) (QA)	234	384	416	416	391	340	259	672	530	519
Runoff Index (QA/PA) (%)	28	38	43	44	40	38	32	56	47	46
BFI	0.67	0.83	0.84	0.9	0.53	0.52	0.78	0.6	0.6	0.61
Geological context	Sandstone	Chalk	Chalk	Chalk	Chalk/Sand	Chalk/Sand	Chalk/Clay	Gneiss/Granite	Gneiss/Granite	Gneiss/Granite
Permeability	Low	High	High	High	Low	Low	Medium	Low	Low	Low

Table 1: Characteristics of the 10 study watersheds.

help identify watersheds which have the strongest influence on the underground compartment. These are BV2, BV3, BV4 and BV7, which have BFI values higher than 0.75 and with “medium” to “high” permeability. There is also a high variability in the runoff index (between 28% and 56%), which is not related to the BFI. A sample of variable areas (100 to 3000 km²) are subject to relatively homogeneous climatology (annual rainfall between 800 and 1191 mm and an average ETP between 610 and 668 mm) but feature varied hydrological functions in terms of the BFI and the runoff index (QA/PA).

3.2. Method for Processing Temporal Signals: Cross-Correlations

Bivariate correlation analyses are based on the study of the correlation between two sets of data, which is classically done in hydrogeology, by studying the rain-flow relationships [33]. To study exchange relationships between surface water and groundwater, a cross-correlation analysis was set up between the flow data at the watershed outlet and the piezometric data in the watershed or nearby. The analysis allowed us to study existing relationships between two temporal variables, according to different temporal shifts. This correlation shows the influence of a first series, called “input data”, on a second series called “output data”. If we consider two-time series x(t) and y(t), the correlogram will have the expression (Equation (3)):

$$r_{xy}\left(k\right)=\frac{C_{xy}\left(k\right)}{{\sigma }_x\ast {\sigma }_y} ;where C_{xy}\left(k\right)=\frac{1}{N}{{\displaystyle \sum }}_{t=1}^{N-k}\left(x\left(t\right)-x_m\right)\left(y\left(t+k\right)-y_m\right)$$

(3)

where:

• $r_{xy}\left(k\right),cross-correlation function$;
• $C_{xy}\left(k\right),cross-correlogram$;
• $k, time lag$;
• $x\left(t\right),input chronicle$;
• $y\left(t\right),output chronicle$;
• ${\sigma }_x et {\sigma }_y,standard deviations of time series x \left(t\right) and y\left(t\right).$

Piezometric data chronicles can be used as an indicator of system dynamics, since they are potentially influenced by the input variable [34]. Correlation analysis results are presented in

Table 2. Summary of Correlation Values r (P, H), r (Q, H) and r (Qb, H) for the Ten Study Watersheds and Associated Time Lag.

Watershed	BV1	BV2	BV3	BV4	BV5	BV6	BV7	BV8	BV9	BV10
r (P, H)	0.17	0.14	0.13	0.14	0.13	0.08	0.07	0.31	0.25	0.19
Time lag r (P, H) (j)	8	92	118	24	9	107	28	8	12	27
r (Q, H)	0.27	0.8	0.74	0.81	0.76	0.54	0.71	0.86	0.83	0.72
Time lag r (Q, H) (j)	145	6	11	2	5	30	8	1	4	9
r (Qb, H)	0.38	0.88	0.78	0.81	0.84	0.75	0.74	0.91	0.87	0.84
Time lag r (Qb, H) (j)	140	−6	0	−19	3	16	−4	−6	−2	4

and Figure 3.

3.3. Selecting the Most Representative Water Table Piezometer—Watercourse Exchange

To fulfil the main objective of this study, in other words, to improve hydrological modelling in a low-water context by taking into account data from the underground compartment (piezometry), the question arises as to whether a single or several piezometers should be studied. As our hydrological modelling is global, piezometric information must also be global. In this case, it will be provided either by a single representative piezometer, or by an average of information from several sites. For the sake of simplicity, we have chosen to select data from only one piezometric station per watershed. This was also the case in the Ruelland study [12], which confronted (via simulation) basic flow dynamics with piezometric dynamics on the Elqui Basin (Chile). However, we will examine the influence this single piezometer may have had on study results, by using information of other piezometers.

A preselection of piezometers was made, according to a cartographic criterion: the Euclidean distance separating the watershed outlet from different piezometers. Only the piezometers closest to the watershed outlet were preselected.

A second criterion was based on recent BRGM studies. We can quote that which relates to the exploitability of groundwater resources in France in the context of resistance to drought of major aquifers with a free water table [35], but also one which relates to determining piezometric indicators for quantitative resource management [36], as well as the study on establishing a piezometric reference network to monitor the impact of climate change on groundwater on a national scale [37]. Another useful study appears in another BRGM -FR report. Here, an atlas of piezometers with a proven groundwater–river relationship was established for each region. The atlas coupled information on the direction of groundwater flow with piezometer distances from their river, classified according to the Strahler order (1 to 6), as defined by the BD-Carthage [38].

Finally, an additional criterion based on cross-correlations was added to choose the reference piezometer of a basin. It involved looking at the correlations between hydro-climatic variables (rainfall (P), discharge (Q) and base flow (Qb)) with piezometric levels (H). Maximum cross-correlation coefficient values with piezometry and their associated time lag are shown in Table 2.

Table 2 shows that base flow is the variable with the highest correlations with piezometry. In fact, with the total flow, we potentially add a bias in the determination of the cross-correlations with piezometry through the peaks of floods, which do not have the same dynamics as the flows outside the rainy period. As for rainfall correlations, their relationship is too weak to be taken into consideration. Temporal shifts between the piezometer level and the station base flow are also variable and can be positive or negative, depending on the slope of the groundwater flow and on the direction of variation (emptying or filling) of other cumulative phenomena.

The correlation criterion used is, therefore, the correlation value between piezometry and the base flow.

Thus, the piezometer chosen as the most representative of the watershed will be the one that best meets the following cartographic and numerical criteria:

⮚

The smallest distance to the main watercourse (highest Strahler order);

⮚

Belonging to the main water body in which the studied watershed is included;

⮚

The highest correlation value between piezometry and base flow (r(Qb, H));

⮚

The smallest distance to the watershed outlet.

Figure 3 shows an example of the selection method conclusions for the most representative piezometer (called reference) on the BV3 (the Aa river). Different piezometric stations are represented on the watershed localisation map. The selected piezometric station is shown in red.

The graph on the left (1) shows the cross-correlation curves between base flows and piezometry for the piezometers closest to the Aa watershed. In red we find the curve corresponding to the piezometer selected as the most representative of groundwater–river exchanges. The green gradient shows r (Qb, H) values of different piezometric stations according to their distance from the watershed outlet. The darker the colour, the closer piezometers are to the watershed outlet; the lighter the colour, the further away they are. Some piezometers have stronger correlations with base flows; some are further away, others are closer to the outlet, with higher or lower time lags. This selection is a compromise between this correlation and the piezometer position in relation to the basin and its water body.

The following part aims to use piezometric information from the reference piezometer to calibrate two surface hydrological models used to determine low water levels.

4. Methodology

Currently, in hydrology, we use conceptual rainfall-runoff models, which consider only one control: the flow data. Weaknesses remain in the structures of these models to correctly simulate the exchanges with the aquifers. A complementary vision consists of integrating an additional control from the underground observation by using piezometric data.

Additionally, as previously mentioned, the objective of this study is to test the contribution of piezometric information when calibrating a hydrological model in a low water context, by using an additional control variable. Another important point that we want to study indirectly is reducing possible equifinality problems.

Two conceptual models are used. The first is a very parsimonious two-parameter model developed to simulate preferentially low water levels. The second is a more parameterised model, with four parameters, which better simulates the whole hydrograph.

For these two models, we tested the influence of adding piezometric information to the calibration performance of the models. For this purpose, the objective function is reviewed in order to integrate a deviation criterion between modelling results and piezometric observations on the basin.

We then evaluated the robustness of the optimised parameter sets and the relevance of the parameterisation to reproduce low water criteria. We also tried to analyse the modelling errors in order to propose improvements in hydrological modelling.

4.1. Hydrological Models

The models used in this study are lumped and parsimonious conceptual models GRLoiEau2J [39] and GR4J [40]. These models have been used in many studies [41,42].

In the GR4J model (Figure 4), effective precipitation (Ps) and actual evapotranspiration (Es) are calculated as a function of the soil moisture storage level (S), net precipitation (P-PE), and the X1 parameter (mm), characterising the maximum capacity of the soil moisture reservoir. Reservoir percolation (Perc) is also a function of the filling rate of the soil moisture reservoir. Interbasin groundwater flows are controlled by a second parameter (X2, mm/d). If X2 is positive, the basin gains water and if X2 is negative, there is a loss of water. The actual rainfall is then divided into two flow components in the routing function: 90% is routed by a unit hydrograph, whose time base parameter is X4 (d), to a nonlinear routing reservoir of parameter X3 (mm) whose discharge leads to the flow Qr, and the remaining 10% is routed by a unit hydrograph to produce the flow Qd. The simulated flow is the sum of these two components (Qr + Qd).

The GRLoiEau2J (Figure 5) model follows the same principle as the GR4J model but using only two parameters. At the level of atmospheric exchanges, parameter A, which describes the moisture reservoir capacity, is regionalised to the entire French metropolitan territory. Balance on the basin is ensured by parameter B, a multiplier coefficient of the outflow. The model does not have a unit hydrograph; the transfer reservoir of maximum capacity C (mm) makes it possible to reproduce the variability of the daily flow.

For the two models used, snow modules are deactivated because the studied watersheds do not have significant snowfall.

4.2. Multi-Criteria and Multi-Forcing Objective Function

Piezometric information will be used for the calibration phase of the models, by introducing it into the objective function. The correlation between piezometric height (H) variations and base flow (Qb) variations will be used to modify the objective function and to adapt it to our objective of low water level simulations.

Indeed, model calibrations were carried out thanks to an objective function, to determine free parameters of the models, so as to find simulated flows which were closest to the observation. A criterion commonly used in hydrology is the Kling–Gupta Efficiency (KGE) criterion [43]:

$$KGE=1-\sqrt{{\left(r-1\right)}^2+{\left(\alpha -1\right)}^2+{\left(\beta -1\right)}^2}$$

(4)

where:

• α is the ratio of simulated and observed flow variances;
• β is the model bias;
• r is the correlation coefficient between observed and simulated flows.

Garcia et al. [44], proposed a new objective function to better simulate low-flow indices. This objective function results from a combination taken as the average between a KGE criterion calculated on flows, which gives greater weight to high flows, and a KGE criterion calculated on the inverse of flows, which gives greater weight to low flows. The authors showed that this objective function provides the best compromise for all the low-flow evaluation criteria studied.

The expression of the objective function FO (Q) is provided by the following formula:

$$FO\left(Q\right)=0.5\times KGE\left(Q\right)+0.5\times KGE\left(\frac{1}{Q}\right)$$

(5)

We have, therefore, chosen to start from this objective function FO (Q) and to add the correlogram information between piezometric level variations and base flow variations. Adding this information, in the new objective function, is realised by a KGE criterion calculated on the correlogram r(Qb, H), for time lags ranging from 0 to 200 days.

The expression of the new objective function FO (Q, H) is as follows:

$$\begin{gathered} FO\left(Q,H\right)=a\times KGE\left(Q\right)+b\times KGE\left(\frac{1}{Q}\right)+c\times KGE\left(r_{Qb,H}\right) \\ \begin{array}{rr}\hfill with& \\ \hfill : a, b, c & \mathrm{respectively}, \mathrm{the} \mathrm{weights} \mathrm{assigned} \mathrm{to} \mathrm{flows}, \mathrm{opposite} \mathrm{flows} \mathrm{and} \hfill \end{array} \\ \mathrm{the} \mathrm{correlogram} r_{Qb,H} \end{gathered}$$

(6)

We will, therefore, test the contribution of this new objective function for the calibration of the models, while optimising the weight of each term.

4.3. Evaluation Criteria

In addition to the restitution of total flows and base flows, evaluation criteria for hydrological modelling will be calculated on low-water indices related to the severity of the low-water period. Average flows can be calculated on a daily time interval (VCNd) or monthly (QMNA). On a daily time interval, they correspond to the annual minimum of a sliding average over d days. The number of days is often chosen as 10 to smooth out the effect of measurement errors on flows or effects from anthropogenic activities [45]. On a monthly time interval, we speak of lowest average monthly flows (QMNA). This index is mainly used in France for drought management.

To evaluate the quality of the calibration of these low-water indices, as well as the total flow and base flows (Qb), a bounded version of the Nash and Stucliffe [46] criterion, called NSE′ [47], has been used:

$$\begin{gathered} NS{E}^{\prime }=\frac{NSE}{2-NSE} \\ with;NSE=1-\frac{{{\displaystyle \sum }}_{i=1}^n{\left(Q_i^{obs}-Q_i^{sim}\right)}^2}{{{\displaystyle \sum }}_{i=1}^n{\left(Q_i^{obs}-{\mu }_{obs}\right)}^2} \\ where {\mu }_{obs} is the average of the chronicles of observed and simulated values \end{gathered}$$

(7)

This criterion was calculated on the QMNA, the VCN10, the base flow (Qb), the correlation between piezometry and the base flow (r (Qb,H)) and the total flow (Q).

This provides assessment criteria on four indices to evaluate the simulation with the new objective function (FO (Q,H)).

5. Results

5.1. Analysis of Choosing Parameter Weights in the Multi-Criteria and Multi-Forcing Objective Function

To evaluate weights to be attributed to each objective function term, a first calibration was carried out on the entire observation period, by varying the value of “c” (weight attributed to the piezometry) from 0.2 to 0.8. In this case, the complement to 1 is distributed equally to the values “a” and “b.” For example, if “c” = 0.2 then “a” and “b” will be equal to 0.4. Following different calibrations, we examined the evolution of different validation criteria in

Table 3. Description of the NSE′ for different evaluation criteria.

Name	Description
NSE′ (QMNA)	Bounded version of the NSE calculated on the QMNA
NSE′ (VCN10)	Bounded version of the NSE calculated on the VCN10
NSE′ r(Qb,H)	Bounded version of the NSE calculated on the correlation between piezometry and the base flow (r(Qb,H))
NSE′ (Qb)	Bounded version of the NSE calculated on base flows (Qb)
NSE′ (Q)	Bounded version of the NSE calculated on the flows (Q)

compared to a calibration with the initial objective function (equivalent to “c” = 0).

Results are presented on Figure 6 graphs for the two models studied.

Figure 6 does not present the value of NSE’ criterion, but the improvement or deterioration by the FO (Q, H) compared to FO (Q).

Figure 6 shows that assigning excessive weight on the value of weight “c” of the FO (Q, H) (0.6 and 0.8), significantly degrades the performances of the two models being validated, for most of the validation criteria.

For the GRLoieau2J model, a lower weight, such as 0.2 or 0.3, will not suggest any real improvements in the validation criteria. A weight of 0.4 seems to be the best compromise between gain and loss on the validation criteria. This weight will, therefore, be kept for this model.

For the GR4J model, the optimal weight retained is 0.2.

The optimal weights for each model having been determined, we examined the temporal robustness of calibrations that were carried out with the new objective function (FO (Q, H)).

5.2. Temporal Robustness

Calibration-validation tests were conducted to evaluate the robustness of the calibration. The calibration algorithm used in this study consists of a two-step search procedure.

First, the parameter space was scanned before running a local search algorithm. This approach will not be discussed here but has proven to be effective for parsimonious models, such as GR4J [48,49]. Evaluating different objective functions was based on a classical split sample test scheme [50]. Flow records were divided into two independent sub-periods of equal length (P1: 1974–1994 and P2: 1995–2015). We first calibrated the model parameters on the P1 period and validated them on P2, and then swapped the two sub-periods, i.e., we calibrated on P2 and approved on P1. Two years before each period (1974–1975 for P1 and 1995–1996 for P2) were used to initialise the model.

The performance of the model is evaluated over the two validation periods. The results for all four evaluation criteria are shown in Figure 7, for the two objective functions tested and for the two hydrological models.

For the GRLoiEau2J model (points in red), low water indices (NSE’QMNA and NSE’VCN 10) are globally better reproduced using the new objective function (FO (Q, H)). The same tendency was observed for NSE’r (Qb,H), where the simulated correlation between the piezometry and the base flow are better reproduced with the new objective function (FO (Q,H).

For the other validation criteria (NSE’Q and NSE’Qb), there are no significant differences between the two objective functions.

Moreover, Figure 8 shows a decreasing trend of GRloiEau2J parameter B and an increasing trend of parameter C. We now know that, when parameter C increases, it implies that we will better reproduce the slow flows. This is consistent with gains observed on the restitution of low water levels.

For the GR4J model, results are more moderate than those observed previously. Indeed, as shown in Figure 9, we see that the implementation of the new objective function FO (Q, H) will provide a better restitution of some evaluation criteria (quite weak), and that a strong degradation of these same criteria was observed. The GR4J model is, therefore, more sensitive to using the new objective function, despite the lower weight attributed to parameter C. Figure 8 shows a stronger evolution of GR4J parameters (X1, X2, X3 and X4), especially the X1 control parameter of the production function.

To conclude on comparing two objective functions for these two models, a summary of the different evaluation criteria is presented in Figure 9.

Adding piezometry data to the new objective function improves performance for the two-parameter GRLoiEau2J model, as shown in Figure 9.

For the GR4J model, the new objective function FO (Q, H) does not offer a significant improvement over the flow objective function (FO (Q)).

The new objective function (FO (Q, H)) thus shows a better restitution of various validation criteria, in particular low-water indices, for the GrLoieau2J model, which is the least parameterised model. This improvement was not found for the GR4J model, as this model is more complex and less parsimonious. Moreover, we observed that weight choices attributed to piezometry in the new objective function must be adapted according to the model used. The more complex the model, the lower the weight assigned to the piezometry control criterion (c) should be.

Significant results observed with the model with few parameters show that ultimately, a different set of parameters can improve calibration criteria without degrading them too much. The contribution of the piezometric information seems to be of interest for this simple model structure, which has more difficulty in representing specific watershed functions, when a more complex structure has additional parameters to take into account. This is the case of the subsurface exchange parameter included in the GR4J model, which proposes nuances to the production function, between what is contributed (lost) by rainfall (evaporation) and what can be contributed or lost by subsurface exchanges.

The improvement observed during the calibration of the GRLoiEau2J model raises the question of the stability of the parameters with respect to the new objective function. An equifinality test will, therefore, be carried out and discussed in the following section. This equifinality test will be associated with the influence of the piezometer choice to determine the impact of this hypothesis of selecting the most representative piezometer on the parameters: how will the parameters perform when a nearby piezometer which does not meet all our selection criteria is chosen?

6. Discussion

6.1. About Equifinality and FO

One of the problems associated with hydrologic model calibration is that of equifinality [25]. Equifinality is present when a number of different parameter sets produce equivalent model performance. This implies that any parameter value determined by calibration will depend on the other parameter values in the model.

Here, we analysed the equifinality of both models with respect to the two studied objective functions. To do so, we performed several model calibrations by successively imposing a parameter at a fixed value. We chose to set the parameter that represents the transfer. For the GRLoiEau2J model, the transfer parameter is parameter C, successively imposed at ten values between 50 and 20,000 mm. The B parameter is thus optimised for these ten imposed C values. For the GR4J model, the X3 parameter was successively imposed at ten values between 20 and 3000 mm. Parameters X1, X2 and X4 were thus optimised for these ten imposed X3 values.

We obtained results for ten optimisations by setting the parameter that represented the watershed flow transfer. Optimisation results correspond to objective function values (FO) and to validation criteria values for the calibration (NSEQ).

Figure 10 and Figure 11 show these different results, for the GRLoiEau2J and GR4J models, respectively. When curves present a plateau of values, it corresponds to equivalent performances for different fixed-parameter values. This helps visualise the equifinality and the discriminatory nature of a given objective function.

The GRLoiEau2J model presents a lower equifinality than the GR4J model, in terms of curves making it easier to determine optimal parameter values. The curves present a more marked maximum than for the GR4J model, which presents wider ranges of equivalent parameters (value plateaus). We also observe that the optimum can be obtained for different parameters depending on the FO used. For example, maximum plateaus for the two objective functions are not identical for some watersheds (BV1, BV2 and BV7). This optimum difference in the FO (Q) and FO (Q, H) justifies the benefit of using the FO (Q, H) for these three watersheds, which may lead to an improvement in the restitution of their validation criterion (Figure 7). For the other seven watersheds, the optimal parameters between the two FOs are identical.

Similar hydrological behaviours are observed for the nested basins in our sample. This is, for example, the case for BV3 and BV4, which will have the same values of NSEQ but also of the two objective functions, while keeping identical optimal values for parameter C.

With the GrLoieau2J model, we also observe that the three catchments with a high permeability (BFI > 0.8 for catchments 2, 3 and 4) present very marked curves to determine the optimum, which then move away from each other according to the FO tested. However, the bias introduced by considering piezometric heights in the FO leads to similar flow restitutions (identical NSEQ). The GR4J model does not highlight this observation.

For the GR4J model, as shown in Figure 11, the two objective functions follow the same evolutionary trend as the NSEQ in the form of plateaus of values while having relatively non-identical optimums of parameter values in calibration.

It is also interesting to note that the set of parameters which present a calibration criterion (FO) optimum is not necessarily the one that would lead to a validation criterion (NSEQ) optimum. Moreover, NSEQ evolution curves are identical, except for BV5 with the GR4J model. The models, therefore, show identical performances for the two objective functions during the calibration of the latter.

6.2. Analysis of Modelling Errors in Calibration

Modelling errors resulting from calibrations with the two tested objective functions are analysed according to observed flow classes and groundwater level fluctuation classes (Figure 12).

From the values of each watershed, four intervals are defined, separated by the quartiles Q1, Q2 and Q3, respectively, associated with the frequencies 25, 50 and 75%, making it possible to divide flows into four classes: [min, Q1]; [Q1, Q2]; [Q2, Q3] and [Q3, max] with, respectively, min and max the minimum and maximum values of the data series.

Figure 12 summarises the results obtained on our watershed sample during the calibration of both models for each of the two objective functions.

Overall, for both models, we do not see a trend in the evolution of modelling errors for the water table fluctuation classes, regardless of the objective function. Errors distributed by flow classes decrease when the flow values increase. This evolution trend is found for GRLoieau2J only with the objective function FO (Q, H). Indeed, for the objective function FO (Q), modelling errors decrease for low and high flows, but increase for medium flows. For GR4J, there is no such distinction between the two objective functions; errors remain similar in terms of median values and follow the same decreasing trend for increasing flow values.

We thus observe a bias of the two models in overestimating low flows and underestimating high flows, in particular when the groundwater–river relationship is taken into account.

6.3. Influence of the Piezometer Choice on the Calibration Parameters of the GRLoiEau2J Model

We were interested in results that would have been obtained by choosing another piezometer. Indeed, other piezometers, similar to that chosen, might have been chosen because they met the criteria for our method. The method presented was then applied by taking into account the data from these different piezometers, in order to calculate the parameters of the GRLoieau2J model. Results are presented in Figure 13 and classify selected piezometers into two categories: “selected,” i.e., piezometers located in the basin with a greater distance to the outlet, and “not selected,” which groups piezometers which are no longer located in the basin but which are simply in the water body with greater distances to the outlet.

Overall, the performance obtained following calibration with other piezometers is slightly worse. By looking at validation criteria results, and those of calibration, Figure 13 shows a slight decrease, depending on the distance of the piezometers from the reference selection criteria (Ref).

Parameter B is sensitive to the distance of a piezometer from the outlet (“Selected”) as well as to the choice of a piezometer outside the criteria (“not selected”). Parameter C is less sensitive to distance, but very sensitive to the selection of a “not selected” piezometer.

In the logic of improving the modelling performance in a low water context, taking into account parameter sensitivity will be even more important. Indeed, in the case of a “not selected” piezometer, whose C parameter significantly decreases, this implies that the model will tend to propose faster flows and be less spread out in time. We also observed a greater variability of parameter C, which shows that the model struggles when taking into account piezometric information that may be more strongly decoupled from base flows. Indeed, for piezometers in the “not selected” criteria and with the greatest variability of parameters, correlations r(Qb, H) are in some cases relatively low, up to 0.1.

7. Conclusions

The aim of this study was to contribute to an additional control derived from observations of subsurface behaviour to current modelling, in order to limit the equifinality problem.

The first tests aimed at improving the calibration of two simple hydrological models (GRLoiEau2J, GR4J), which took into account subsurface data using a new objective function, were conclusive only for the simplest model with two parameters, the GRLoiEau2J model. For the four-parameter model, the GR4J model, flow restitution, as well as low-water index restitution, did not contribute any significant gain, and was even degraded in some cases.

These results can be linked to the simplicity of the GRLoiEau2J model, a very parsimonious global model with two free parameters, B and C. The results of the equifinality analysis also showed a very low equifinality for the two parameters B and C. In this case, the model cannot converge to different solutions without imposing an additional constraint. This constraint is brought here by taking into account piezometric information, making it possible to slightly orientate the model calibration to put more weight on modelling the groundwater–river exchange. For the GR4J model, which has more equifinality, introducing an additional constraint might have led to a different optimisation. Surprisingly, this was not the case. The calibration with the two tested objective functions led to the same parameterisation. In fact, the introduction of a constraint in the underground exchanges did not seem necessary for the GR4J model, which already has a parameter that allows it to take into account underground exchanges.

The difficulty of the GR4J model for simulating low water indices can be explained by the presence of a control parameter for exchanges (X2) within the modelling chain. Indeed, in other studies [51,52], similar conclusions were found. For less parsimonious models (GR4J, GR5J, GR6J) than GRLoiEau2J, which integrate one or more parameters that regulate underground and inter-basin exchanges to their structure, adding piezometric information does not seem to provide relevant information. These models are already able to model different watershed losses or gains with the underground compartment due to their exchange parameters (X2, X5). On the other hand, for a more parsimonious model, introducing a constraint in the optimisation function can modify the parameterisation to better take into account exchanges with the underground compartment. This can be considered as an important result to improve the performance of very parsimonious models.