Improvement of the KarstMod modelling platform for a better assessment of karst groundwater resources

Hydrological models are fundamental tools for the characterization and management of karst systems. We propose an updated version of KarstMod, software dedicated to lumped-parameter rainfall–discharge modelling of karst aquifers. KarstMod provides a modular, user-friendly modelling environment for educational, research, and operational purposes. It also includes numerical tools for time series analysis, model evaluation, and sensitivity analysis. The modularity of the platform facilitates common operations related to lumped-parameter rainfall–discharge modelling, such as (i) setup and parameter estimation of a relevant model structure and (ii) evaluation of internal consistency, parameter sensitivity, and hydrograph characteristics. The updated version now includes (i) external routines to better consider the input data and their related uncertainties, i.e. evapotranspiration and solid precipitation; (ii) enlargement of multi-objective calibration possibilities, allowing more flexibility in terms of objective functions and observation type; and (iii) additional tools for model performance evaluation, including further performance criteria and tools for model error representation.

1 Introduction

Karst aquifers constitute an essential source of drinking water for about 9.2 % of the world population (Stevanović, 2019), and it is estimated that one-quarter of the world population depends on freshwater from karst aquifers (Ford and Williams, 2013). Karst aquifers contain an important volume of freshwater, but only 1 % of their annually renewable water is used for drinking-water supply (Stevanović, 2019). Understanding the functioning of karst aquifers and developing operational tools to predict the evolution of freshwater resources is therefore a major challenge for the hydrological science community (Blöschl et al., 2019). To this day, the number of tools dedicated to karst hydrogeology is limited; these tools are also mostly developed for academic purposes and are not user-friendly. Nonetheless, such tools are required for a better assessment of groundwater vulnerability, as well as for sustainable management of the groundwater resources (Elshall et al., 2020), and should be handled by stakeholders without programming-skill requirements.

KarstMod is an adjustable modelling platform (Mazzilli et al., 2019) dedicated to lumped-parameter rainfall–discharge modelling, allowing for (i) simulation of spring discharge, piezometric head, and surface water discharge (Bailly-Comte et al., 2010; Cousquer and Jourde, 2022; Sophocleous, 2002); (ii) analysis of the internal fluxes considered in the model; and (iii) model performance evaluation and parametric sensitivity analysis. In this paper, we present the new features incorporated into KarstMod: (i) external routines to better consider the input data and their related uncertainties, i.e. evapotranspiration and solid precipitation; (ii) enlargement of multi-objective calibration possibilities, allowing more flexibility in terms of objective functions and observation type, with the possibility of including surface water discharge in the calibration procedure; and (iii) model performance evaluation, including additional performance criteria, as well as additional tools for model error representation, such as the diagnostic efficiency plot (Schwemmle et al., 2021). Also, we present two case studies to illustrate how KarstMod is useful in the framework of the assessment of karst groundwater resources and its sensitivity to groundwater abstraction. Section 2 is devoted to the presentation of the background and motivations to improve the functionalities of the platform, while Sect. 3 presents the key features of KarstMod. Section 4 illustrates the application of rainfall–discharge modelling using KarstMod within the Touvre (western France) and the Lez (southern France) karst systems, which both constitute strategic freshwater resources and ensure drinking-water supply.

2 Background and motivations

2.1 Challenges in karst groundwater resources

Karst aquifers are affected by the combination of different components of global change, such as (i) the effects of climate change which are particularly pronounced in the Mediterranean area (Dubois et al., 2020; Nerantzaki and Nikolaidis, 2020), (ii) increasing groundwater abstraction (Labat et al., 2022), and (iii) changes in land cover and land use (Bittner et al., 2018; Sarrazin et al., 2018). Therefore, the assessment of karst groundwater resource sensitivity, specifically in terms of quantity, requires operational tools to estimate the sustainable yield of karst aquifers but also to predict the impacts of climatic or anthropogenic forcing on groundwater resources in the long term (Sivelle et al., 2021). To address these issues, different modelling approaches have been developed (Jeannin et al., 2021), such as, among others, fully distributed models (Chen and Goldscheider, 2014), semi-distributed models (Doummar et al., 2012; Dubois et al., 2020; Ollivier et al., 2020), and lumped-parameter models (Mazzilli et al., 2019) including semi-distributed recharge (Bittner et al., 2018; Sivelle et al., 2022b). Among these, lumped-parameter models are recognized as major tools to explore the ability of conceptual representations to explain observations in karst systems (Duran et al., 2020; Frank et al., 2021; Poulain et al., 2018; Sivelle et al., 2019) and to manage karst groundwater resources (Cousquer and Jourde, 2022; Labat et al., 2022; Sivelle et al., 2021; Sivelle and Jourde, 2020).

2.2 Challenges in lumped-parameters modelling in karst hydrology

Lumped-parameter models consist of a functional approach that analyses a hydrogeological system at the catchment scale and describes the transformation from rainfall into discharge using empirical or conceptual relationships. Therefore, parameter values or distributions cannot be determined directly from catchment physical characteristics or in situ measurements, except for the discharge coefficient in relation to the spring, which can be estimated based on recession curve analysis. Instead, model parameter values must be estimated by history-matching. In a general way, rainfall–discharge models in karst hydrology are calibrated considering spring discharge measurements. Former studies have shown in interest in considering hydrochemical observations (Chang et al., 2021; Hartmann et al., 2013; Sivelle et al., 2022a), but such an approach requires further methodological development before being included in KarstMod. To date, KarstMod allows for a consideration of complementary observations only with piezometric heads and surface water discharge (Cousquer and Jourde, 2022).

Another challenge concerns the evaluation of the water fluxes within the soil–vegetation–atmosphere continuum. Bittner et al. (2021) computed several models to evaluate the fluxes related to interception, evapotranspiration, and snow process. The results show significant uncertainties related to input data, as well as potential compensation between the various uncertain processes. In some cases, snowmelt is a controlling factor in the water balance (Doummar et al., 2018; Liu et al., 2021); thus, a suitable snowmelt estimation is required to improve hydrological model performance (Çallıet al., 2022). Therefore, two meteorological modules have been added to KarstMod: (i) a “snow routine” and (ii) a routine to compute the potential evapotranspiration (PET, mm d⁻¹), denoted the “PET routine”. The two additional modules allow us to better account for snow and evapotranspiration processes.

3 Implementation

The updated version of KarstMod implements additional features to enhance the rainfall–discharge modelling practices. First, we describe the additional modules (snow and PET routines) for a better meteorological forcing estimation. Then, we introduce the additional tools proposed for (i) the setup and calibration of the model structure, (ii) model performance evaluation, and (iii) uncertainty consideration. Fig. 1 shows a screenshot of the KarstMod software.

Figure 1 Screenshot of the KarstMod software with (a) model structure, (b) data import, (c) model parameters, (d) run parameters, (e) calibration results, (f) command bar, and (g) results and graphs.

3.1 Meteorological modules

3.1.1 Snow routine

KarstMod allows using either observation-based precipitation time series P (mm d⁻¹) or estimated precipitation time series P_sr (mm d⁻¹) using a snow routine. The latter is similar to the one used by Chen et al. (2018) – without the radiation components – which has been successfully used for improving the simulation of karst spring discharge in snow-covered karst systems (Chen et al., 2018; Cinkus et al., 2023a). It consists of a modified HBV–snow routine (Bergström, 1992) for simulating snow accumulation and melt over different sub-catchments based on altitude ranges. Each sub-catchment is defined by two values that the user must input: (i) the proportion among the whole catchment (sum must be equal to 1) and (ii) the temperature shift, related to the altitude gradient. The different estimated precipitation values ($P_{\mathrm{sr}}^{\ast }$, mm d⁻¹) associated with the sub-catchments are calculated and summed to produce the estimated precipitation time series P_sr, which corresponds to a single variable representative of the catchment. P_sr thus gives the water leaving the snow routine and is equivalent to the recharge into the first compartment of the model (compartment E in KarstMod). The snow routine workflow requires both air temperature T (°C) and precipitation P (mm d⁻¹) time series. P is considered to be snow when T in the sub-catchment is lower than the temperature threshold T_s (°C). Snow melts when the temperature exceeds the threshold according to a degree-day expression. The snowmelt is a function of the melt coefficient MF (mm °C⁻¹ d⁻¹) and the degrees above the temperature threshold. Runoff starts when the water level exceeds the liquid-water-holding capacity of snow CWH (−). The refreezing coefficient CFR (−) stands for refreezing liquid water in the snow when snowmelt is interrupted (Bergström, 1992). The output of the snow routine consists of a redistributed precipitation time series P_sr. The four parameters of the snow routine (i.e. T_s, MF, CWH, and CFR) can be considered in the parameter estimation procedure, as well as in the sensitivity analysis. The snow routine features can be activated from the model structure area (Fig. 1a). Figure 2 shows the general workflow implemented in the snow routine. $P_{\mathrm{sr}}^{\ast }$ is estimated for each time step t based on the precipitation P and air temperature T time series for each sub-catchment i. The total snow routine output P_sr is calculated as a weighted sum of $P_{\mathrm{sr}}^{\ast }$ time series:

$$\begin{array}{}\text{(1)}& P_{\mathrm{sr}}=\sum _i^N{P}_{{\mathrm{sr}}_i}^{\ast }\times p_i,\end{array}$$

where p_i is the proportion of the sub-catchment i in relation to the complete catchment area, such as $\sum p_i=\mathrm{1}$, and N is the total number of sub-catchments. The snow routine allows for an estimation of $P_{\mathrm{sr}}^{\ast }$ according to Algorithm 1.

Figure 2 Snow routine workflow.

3.1.2 Potential evapotranspiration routine

An additional module allows us to compute the potential evapotranspiration PET (mm d⁻¹) based on Oudin's formula (Oudin et al., 2005). The PET routine can be activated from the model structure area (Fig. 1a). The PET routine affects only compartment E. The latter stands for soil and epikarst storage zone, where the water is available for actual evapotranspiration AET (mm d⁻¹) and flows toward infiltration or surface discharge. Infiltration occurs when the water level in the compartment is greater than a given threshold Emin; otherwise, the compartment is considered to be under-saturated and does not produce infiltration. In this case, the water in compartment E is still available for evapotranspiration. KarstMod allows us to consider evapotranspiration in four separate ways (Fig. 3):

• a.

  The water transfer in the soil–atmosphere continuum can be pre-processed by the user. In this case, the given precipitation time series consists of the effective precipitation P_eff (mm d⁻¹), derived from precipitation P (mm d⁻¹) and actual evapotranspiration AET (mm d⁻¹) with Eq. (2). The evapotranspiration flux is not activated in the model structure selection panel in KarstMod (Fig. 1a).

  $$\begin{array}{}\text{(2)}& P_{\mathrm{eff}}=P-\mathrm{AET}\end{array}$$
• b.

  User-defined PET can be given as input in KarstMod for the evapotranspiration time series. Using Emin, the user can simulate water-holding capacity and non-linear behaviour of karst recharge.

• c.

  User-defined AET can be given as input data in KarstMod for evapotranspiration time series instead of PET. Then, KarstMod computes an estimation of effective precipitation by limiting the evapotranspiration to water content available in compartment E. The simulated AET can then be lower than the user-defined AET. Such a configuration may help in identifying the potential inaccuracy of user-defined AET for the modelling purpose but is not recommended for model setup and parameter estimation.

• d.

  The new feature in KarstMod consists of the PET routine which estimates the PET with Oudin's formula (Oudin et al., 2005) (Eq. 3). It needs a T time series and two parameters to be estimated, which can be considered in the parameter estimation procedure, as well as in the sensitivity analysis.

  $$\begin{array}{}\text{(3)}& \begin{array}{rl}\mathrm{PET}& =\left({\displaystyle \frac{R_{\mathrm{e}}}{\mathit{\lambda }\times \mathit{\rho }}}\right)\times \left({\displaystyle \frac{T+K\mathrm{2}}{K\mathrm{1}}}\right)\text{if}\\ T& +K\mathrm{2}>\mathrm{0}\text{else}\mathrm{PET}=\mathrm{0},\end{array}\end{array}$$

  where R_e is the extraterrestrial radiation (MJ m⁻² d⁻¹) depending only on the latitude and the Julian day, λ is the latent heat flux (taken to be equal to 2.45 MJ kg⁻¹), ρ is the density of water (taken to be equal to 1000 kg m⁻³), and T is the mean daily air temperature (°C). K1 (°C) and K2 (°C) are constants to adjust over the catchment for the rainfall–discharge model (Oudin et al., 2005). In KarstMod, both K1 and K2 can be considered in the parameter estimation procedure, as well as in the sensitivity analysis.

Figure 3 The four ways to account for evapotranspiration in KarstMod. The user can provide either (a) a self-computed effective precipitation (P−AET) as a single input time series, (b) both P and PET time series, (c) both P and AET, or (d) both P and T time series.

3.2 Setup and calibration of the model structure

The modular structure proposed in KarstMod is based on a widely used conceptual model which separates karst aquifers into an infiltration zone and a saturated zone or low and quick flows through the unsaturated and saturated zones (Fleury et al., 2007, 2009; Guinot et al., 2015; Mazzilli et al., 2019; Sivelle et al., 2019). Based on this conceptual representation, the platform offers four compartments organized as a two-level structure: (i) compartment E (higher level) and (ii) compartments L, M, and C (lower level). A priori, the higher level represents the infiltration zone or the soil and epikarst. At the lower level, compartments L, M, and C stand for the different sub-systems of the saturated zone or low and quick flows of the whole hydrologic system. The various model structures and their governing equations are presented in Mazzilli et al. (2022, 2019). Also, KarstMod allows us to perform hydrological modelling at both daily and hourly temporal resolutions (Sivelle et al., 2019).

The user can activate (or deactivate) the various compartments (E, L, M, and C) within the “model structure” panel (Fig. 1a). The solid and faded colours represent the activated and the inactivated features, respectively. The fluxes and their activation thresholds, as well as the exponent of the discharge law α (in the case of a non-linear discharge law such as α ≠1), are managed based on the “model parameters” panel (Fig. 1c). The user can account for pumping Q_pump (water coming out of the compartment), as well as sinking stream Q_sink (water coming into the compartment). Such an option is available only if the user provides the required time series (Fig. 1b).

The user must provide warm-up, calibration, and validation periods (Fig. 1d). The warm-up period must be set to be independent of initial conditions to avoid bias in the parameter estimation procedure (Mazzilli et al., 2012). Then, a calibration period (i.e. the period in which the parameters are estimated to reduce the predictive errors) and a validation period (i.e. a period separated from the calibration period) can be defined to run the split-sample test procedure (Klemeš, 1986). For calibration purposes, KarstMod proposes several widely used performance criteria ϕ: the Pearson's correlation coefficient R_p (Freedman et al., 2007), the Spearman rank correlation coefficient R_s (Freedman et al., 2007), the Nash–Sutcliffe efficiency NSE (Nash and Sutcliffe, 1970), the volumetric error VE (Criss and Winston, 2008), the modified balance error BE (Perrin et al., 2001), the Kling–Gupta efficiency KGE (Gupta et al., 2009), and a non-parametric variant of the Kling–Gupta efficiency KGENP (Pool et al., 2018). To compute a multi-objective calibration procedure, the user can create their objective function Φ as a weighted sum of several objective functions:

$$\begin{array}{}\text{(4)}& \mathrm{\Phi }=\sum _{i=\mathrm{1}}^N{\mathit{\omega }}_i\times {\mathit{\varphi }}_i\left(U\right),\end{array}$$

where ω_i is the weight affected in relation the objective function ϕ_i(U), with $\sum _{i=\mathrm{1}}^N{\mathit{\omega }}_i=\mathrm{1}$, and U is a general notation for the observations used for parameter estimation purposes. In the KarstMod modelling platform, U corresponds to either spring discharge Q_s, piezometric-head measurements Z_obs (available for compartments E, L, M, and C), or surface water discharge Q_loss from compartment E. Also, the objective function ϕ can be computed based on transformed U to avoid high water level bias in terms of quadratic error. The following transformations are available in KarstMod: $\mathrm{1}/U$, $\sqrt{U}$, and $\mathrm{1}/\sqrt{U}$. Therefore, the user can use any combination of the objective function ϕ, the observations U, and the variable transformations. Depending on the modelling purpose, the user must refer to the literature to define the suitable objective function (Bennett et al., 2013; Ferreira et al., 2020; Hauduc et al., 2015; Jackson et al., 2019).

The model is calibrated using a quasi-Monte-Carlo sampling procedure with a Sobol sequence sampling of the parameter space (Sobol, 1998). The procedure involves finding an ensemble of parameter sets providing an objective function Φ greater than the user-defined value. The calibration procedure is stopped when either the user-defined maximum duration of the sampling procedure t_max is reached or the user-defined number of parameter sets n_obj are collected. KarstMod offers a “run” option, allowing the model to run for a user-defined parameter set without any calibration procedure, and so allows it to investigate “by hand” the parameter space and the sensitivity of the model to specific parameters.

3.3 Model evaluation

The model performance can be evaluated for both the calibration and validation periods. It allows us (i) to ensure the robustness of model predictions, even under changing conditions (which is a key point for the assessment of climate change impact) and (ii) to avoid model over-fitting within a specific range of hydro-climatic conditions observed during the calibration period. KarstMod allows the computation of the above-mentioned performance criteria for both the calibration and validation periods. Even though the notation “validation” is disputable, such a procedure is required to evaluate both explanatory and predictive dimensions of the model structure (Andréassian, 2023). Then, KarstMod offers an ensemble of numerical tools devoted to (i) checking the model consistency, i.e. the explanatory dimension of the model (Beven, 2001; Shmueli, 2010), and (ii) evaluating the model performance, i.e. the predictive dimension of the model structure.

To check the model consistency, the simulation based on the parameter set that provides the highest objective function value can be analysed through an ensemble of graphs, such as (i) internal and external fluxes as a function of time; (ii) cumulative volumes for both observed and simulated time series for spring discharge Q_s and surface water discharge Q_loss; (iii) simulated mass balance as a function of time; and (iv) a comparison of observations and simulations for either Q_s or Q_loss with probability function plots, auto-correlograms of the spring discharge time series, and cross-correlograms of precipitation–discharge time series.

To evaluate the model performance, KarstMod offers a “model evaluation” panel, available from the graphs panel (Fig. 1g), that includes several sub-panels (from the left to the right):

• The diagnostic efficiency DE (Schwemmle et al., 2021) consists of a diagnostic polar plot that facilitates the model evaluation process, as well as the comparison of multiple simulations. The DE accounts for constant, dynamic, and timing errors and their relative contributions to the model errors. Also, the decomposition of the errors between the periods of high flows and low flows allows us to better investigate the model bias, as well as to provide a critical evaluation for impact studies, particularly for the assessment of climate change impacts. Indeed, the accurate evaluation of low-flow periods (in terms of frequency, intensity, and duration) becomes increasingly crucial for groundwater resource variability assessment.

• The available objective functions Φ are presented as a radar chart which consists of a polygon where the position of each point in relation to the centre gives the value of the performance criteria. The closer the point is to the outside of the radar chart, the better the model performs. The radar chart is made for both calibration and validation periods and for each of the calibration variables considered in the modelling (Q_s; Z_obsA with A for either E, M, C, or L compartments; and Q_loss).

• The KGE (Gupta et al., 2009) consists of a diagonal decomposition of the NSE (Nash and Sutcliffe, 1970) to separate Pearson's correlation coefficient R_p, the representation of bias β_KGE, and the variability α_KGE. Thus, the KGE is comparable to multi-objective criteria for calibration purposes (Pechlivanidis et al., 2013). The sub-panel offers (i) a bi-plot of the three KGE components and (ii) a radar plot visualization of the KGE components, allowing for the identification of potential counterbalancing errors according to these different components (Cinkus et al., 2023b). The two above-mentioned plots also include the decomposition of the KGENP (Pool et al., 2018) in terms of Spearman's rank correlation coefficient R_s, representation of bias β_KGENP, and non-parametric variability α_KGENP.

3.4 Dealing with uncertainties

Moges et al. (2021) summarize the various sources of uncertainty in hydrological models, including structural and parametric uncertainties, as well as uncertainties related to input data and observations. The latter concerns both the input (i.e. precipitation and evapotranspiration) and the output (i.e. discharge) of the modelled systems. Many references are devoted to the uncertainties related to input data and observations. As an example, Westerberg et al. (2020) include information about the discharge uncertainty distribution in the objective function and perform better discharge simulations. Also, the precipitation error can be dependent on the data time step (McMillan et al., 2011) and could impact the hydrological model performance (Ficchì et al., 2016). Lumped-parameter hydrological models consider meteorological time series representative of a whole catchment, which may require some pre-processing, particularly for snow processes since these can have a strong influence on flow dynamics. Thus, KarstMod includes variables related to both the snow routine (i.e. the redistributed precipitation time series P_sr) and the PET routine (i.e. estimated potential evapotranspiration PET) in the parameter estimation procedure. This allows us to investigate the sensitivity of the flow simulation to these input data when using snow and PET routines. Nonetheless, KarstMod does not include features to investigate the impact of observation uncertainties on parameter estimation.

As with many environmental problems, parameter estimation in rainfall–discharge modelling consists generally of ill-posed problems, i.e. the modelling encounters issues with regard to the unicity, identifiability, and stability of the problem solution (Ebel and Loague, 2006). As a consequence, several representations of the modelled catchment may be considered to be equally acceptable (Beven, 2006). Knoben et al. (2020) evaluate the performance of 36 daily lumped-parameter models over 559 catchments and show that between 1 and up to 28 models can show performances close to the model structure with the highest performance criteria. Such results are widely covered in catchment hydrology (Dakhlaoui and Djebbi, 2021; Darbandsari and Coulibaly, 2020; Gupta and Govindaraju, 2019; Pandi et al., 2021; Zhou et al., 2021) but are still poorly investigated in karst hydrology. Indeed, the structural-uncertainty impacts on rainfall–discharge modelling in karst hydrology are not properly evaluated, whereas many studies consider several hydrological model structures to include structural uncertainty in flow simulation (Hartmann et al., 2012; Jiang et al., 2007; Jones et al., 2006; Sivelle et al., 2021). KarstMod includes more than 50 combinations of the various compartments, as well as various compartment models (i.e. compartments with linear or non-linear discharge laws and compartments with infinite characteristic time), and allows a quick implementation of the various model structures. The user can easily manage to start the modelling with one single compartment and gradually move to a more complex model structure with up to four compartments, five fluxes connected to the spring, four internal fluxes, and one flux running out of the system.

Considering each model structure, parametric equifinality can be investigated using (i) dotty plots of the values of the objective function against the parameter values; (ii) dotty plots of the values of the performance criteria used to define the aggregated objective function; and (iii) the variance-based, first-order S_i and total ST_i sensitivity indexes for the model parameters. Details concerning the computation of sensitivity indexes within KarstMod are given in Mazzilli et al. (2022, 2019).

4 Examples of application

To illustrate the KarstMod application and the use of the above-presented functionalities for the assessment of karst groundwater resources, we propose two case studies: (i) the Touvre karst system and (ii) the Lez karst system. Both karst systems consist of strategic freshwater resources for drinking-water supply (DWS) for the cities of Angoulême (western France) and Montpellier (southern France).

4.1 The Touvre karst system (La Rochefoucauld)

The Touvre karst system is one where the infiltration consists of (i) a delayed infiltration of effective precipitation into the karstic recharge area and (ii) a direct infiltration of surface water from the Tardoire, Bandiat, and Bonnieure rivers. The latter are surface streams flowing on metamorphic rocks that partly infiltrate into subterranean levels upon contact with carbonate formations, mainly composed of Middle to Upper Jurassic limestones. The springs of the Touvre, located 7 km east of Angoulême (western France), are made up of four outlets, namely the Bouillant, the Dormant, the Font de Lussac, and the Lèche (Labat et al., 2022). In the following, the Touvre Spring discharge designates the accumulated discharge of the four mentioned outlets.

The Touvre karst system constitutes a strategic freshwater resource for the DWS of Angoulême, with around 110 000 inhabitants, but also contributes to the water supply for industry and agriculture. In 2015, there were 84 pumping wells over the karstic impluvium of the Touvre karst system and around 100 more in the Tardoire, Bandiat, and Bonnieure river catchments. Based on the data provided by the Adour-Garonne Water Agency, the annual groundwater abstraction for agriculture represents 4.6 Mm³, whereas annual groundwater abstraction for DWS represents 1.1 Mm³ over the karstic impluvium of the Touvre karst system. On the three river catchments (out of the karstic impluvium), the annual groundwater abstraction represents 2.5 Mm³ for agriculture and 3.3 Mm³ for DWS through river intakes or alluvial groundwater abstraction. The total annual volume of abstracted groundwater in the area represents around 5 % of the annual volume of transit at the Touvre Spring. This is quite low compared with karst aquifers in France exploited for their groundwater resources, such as the Lez Spring (Jourde et al., 2014) and the Oeillal Spring karst catchments (Sivelle et al., 2021), where the annual groundwater abstraction volume represents, respectively, 50 % and 15 % of the annual volume of transit at the springs. Therefore, the Touvre karst system seems not to be over-exploited at the moment, but the impact of groundwater abstraction should be addressed in the context of global change to ensure sustainable management of this strategic freshwater resource.

The area is characterized by an ocean-influenced climate, with a mean annual precipitation of around 800 mm yr⁻¹, distributed over an average of 255 rainy days. The estimation is performed with Thiessen polygon methods based on 11 meteorological stations over the area (Labat et al., 2022). The mean annual potential evapotranspiration is around 770 mm yr⁻¹ according to the Penman–Monteith estimation provided by the French meteorological survey (Météo-France). The Touvre daily spring discharge shows a significant variability ranging from 3 to 49 m³ s⁻¹, with a coefficient of variation around 0.46 (Fig. 5b).

The surface streamflow rates for the Bonnieure, Bandiat, and Tardoire rivers are concentrated within the autumn and winter periods. During the summer period, the discharge in the three rivers is very low (Fig. 5c). The more significant groundwater abstraction is performed during the summer period, while the Touvre Spring discharge reaches its lowest values within the late-summer and early-autumn periods (Fig. 5c and d).

Figure 4 shows the model structure for the Touvre karst system that consists of three compartments organized into two levels (Labat et al., 2022). The upper level corresponds to reservoir E and represents both the unsaturated part of the system and a temporary aquifer. This reservoir relates to the two reservoirs of the lower level: C (conduit) and M (matrix), representative of quick and slow flow dynamics, respectively. The upper level of the model structure is affected by P and ET, while the lower level of the model structure is affected by (i) groundwater abstraction and (ii) sinking river streamflow from the surface to underground. Figure 4 shows the various time series required for the hydrological modelling of the Touvre karst system. The methodology for daily time series preparation given in Labat et al. (2022) allows us to account for the influence of groundwater abstraction on the transmissive or capacitive part of the karst aquifer, as well as the influence of the concentrated and diffuse infiltration of the surface river streamflow.

Figure 4 Screenshot of KarstMod, with a focus on the panel “model structure” for the Touvre karst system. The solid lines correspond to the activated fluxes, whereas the faded coloured lines are not activated. $Q_{\mathrm{p}}^{\mathrm{M}}$ and $Q_{\mathrm{p}}^{\mathrm{C}}$ stand for groundwater abstraction that affects compartments M and C, respectively, while $Q_{\mathrm{s}}^{\mathrm{M}}$ and $Q_{\mathrm{s}}^{\mathrm{C}}$ stand for sinking flow that affects compartments M and C, respectively.

Figure 5 Daily time series for the Touvre system: (a) precipitation (P) and potential evapotranspiration (PET), (b) observed and simulated karst spring discharge (Q_Touvre obs and Q_Touvre sim), (c) observed river streamflow discharge (Q_Bonnieur, Q_Bandiat, Q_Tardoire), and (d) and (e) groundwater abstraction discharge ($Q_{\mathrm{p}}^{\mathrm{aggriculture}}$ and $Q_{\mathrm{p}}^{\mathrm{domestic}}$).

The objective of the hydrological modelling is to assess the impact of groundwater abstraction on spring discharge, more particularly during low-flow periods (Labat et al., 2022). Thus, the calibration is performed according to the KGENP, which improves the simulations during mean- and low-flow conditions using the Spearman rank correlation due to its insensitivity to extreme values (Pool et al., 2018). The sampling procedure is set up to find n_obj=5000 simulations with KGENP values greater than 0.9. Afterwards, the model is evaluated using the various features proposed in KarstMod (Fig. 6). The diagnostic efficiency plot (Fig. 6a) testifies to several elements: (i) the model seems to slightly overestimate high flow and underestimate low flow; (ii) the timing error is about 0.9, testifying to suitable flow dynamics in the model; (iii) low-flow periods contribute more to the model errors; and (iv) there is no offset in the simulated spring hydrograph. The radar chart (Fig. 6b) shows a good equilibrium between the various objective functions whose values are greater than 0.8, except for the NSE criteria (NSE = 0.75). This is the consequence of the design of these criteria, which tends to outweigh the errors during floods. Here, the NSE value is still greater than 0.7 and testifies to a “very good” fit according to Moriasi et al. (2007). Finally, the decomposition of the KGE (Fig. 6c and d) shows R_p=0.91, α=1.15, and β=1.02, testifying too accurate dynamics and low bias but variability that is slightly too high.

Figure 6 Screenshot of KarstMod with a focus on the sub-panel “model evaluation”. Application for the model evaluation on the Touvre system: (a) diagnostic efficiency plot (Schwemmle et al., 2021), (b) radar chart of the objective functions, (c) bi-plot of the KGE (square) and KGENP (triangle) components, and (d) radar chart of the KGE components.

4.2 The Lez Spring

The Lez Spring (southern France) consists of the main outlet of a karst system situated in the North Montpellier Garrigue hydrogeological unit, delimited to the west by the Hérault River and to the north and east by the Vidourle River. The geology in the area corresponds to the Upper Jurassic layers separated by the Corconne–Matelle fault (oriented N30°), leading to the two main compartments in the aquifer (Bérard, 1983; Clauzon et al., 2020). The karst aquifer is unconfined in the western compartment and is locally confined in the eastern compartment. The Lez Spring is located about 15 km north of Montpellier. It is of Vauclusian-type, with an overflow level at 65 m a.s.l. and a maximum daily discharge of approximately 15 m³ s⁻¹. The area is characterized by a typical Mediterranean climate, with dry summers and rainy autumns. Over the 2009–2019 period, the mean annual precipitation is around 900 mm yr⁻¹, distributed over an average of 133 rainy days (estimation with Thiessen polygon methods based on four meteorological stations over the area: Prades-le-Lez, Saint-Martin-de-Londres, Sauteyrargues, and Valflaunès). The mean annual potential evapotranspiration is around 900 mm yr⁻¹ according to the estimation based on Oudin's formula, with the temperature measured at Prades-le-Lez station, while the real annual evapotranspiration is around 450 mm yr⁻¹ (eddy covariance flux station of Puéchabon).

Since 1854, the Lez Spring has supplied the drinking water to the city of Montpellier and its surroundings. It currently constitutes the main freshwater resource for around 350 000 people in the area. The present water management scheme allows for pumping at higher rates than the natural-spring discharge during low-flow periods while supplying a minimum discharge rate (around 0.23 m³ s⁻¹) into the Lez River to ensure ecological flow downstream and reducing flood hazards via rainfall storage in autumn (Avias, 1995; Jourde et al., 2014). The pumping plant was built in 1982, with four deep wells drilled to intercept the karst conduit feeding the spring, 48 m below the overflow level of the spring. Pumping in these wells allows for up to 0.18 m³ s⁻¹ to be withdrawn during low-flow periods (with an authorized maximum drawdown of 30 m), while the average annual pumping flow rate is about 0.10 m³ s⁻¹ (over the 2008–2019 period). Due to the pumping management of the aquifer, which supplies about 30 to 35 Mm³ of water per year to the metropolitan area of Montpellier, the discharge at the Lez Spring is often low or nil. Discharge is also measured downstream (Lavalette gauging station), where the measured discharge corresponds to the Lez Spring discharge and the main tributaries (Lirou and Terrieu streams), which flow essentially after intense Mediterranean rainfall events. As suggested in Cousquer and Jourde (2022), the surface water discharge, denoted Q_loss, can be estimated as the difference between the total discharge in Lavalette and the Lez Spring discharge.

In the present context of global change, Mediterranean karst systems already show significant decreases in spring discharge (Doummar et al., 2018; Dubois et al., 2020; Fiorillo et al., 2021; Hartmann et al., 2012; Nerantzaki and Nikolaidis, 2020; Smiatek et al., 2013), which could be aggravated with groundwater abstraction (Sivelle et al., 2021). The Lez Spring is strongly exposed to global climate change impacts: (i) the Mediterranean area has been identified as a climate change hotspot (Diffenbaugh and Giorgi, 2012) where the projected warming spans 1.8–8.4 °C according to CMIP6 and 1.2–6.6 °C according to CMIP5 during the summer period (Cos et al., 2022), and (ii) the water management scheme will have to adapt to the future needs for drinking water of the growing population in the area, as well as to changes in the freshwater consumption practices (e.g. water use restriction order). Therefore, a sustainable water management plan for the Lez Spring requires a good appreciation of the hydrological functioning and of the operational hydrological model to properly address impact studies. In this framework, KarstMod allows for choosing and calibrating a suitable model structure. This constitutes the first step for a global climate change impact study that requires prediction tools to simulate the aquifer response to various external forces.

Figure 7 shows the model structure for the Lez karst catchment (Mazzilli et al., 2011) that consists of three compartments organized into two levels. The upper level corresponds to compartment E and represents the unsaturated part of the system, including soil-water-holding capacity Emin and discharge lost from the compartment Q_loss. Compartment E is exposed to P and ET, and discharge towards the lower level of the model structure starts when the water level exceeds Emin. The lower level consists of two inter-connected compartments, M and C, allowing us to reproduce the lateral exchanges, denoted Q_MC, between the transmissive function (compartment C) and the capacitive function (compartment M) of the karst aquifer. Both the M and C compartments are considered to be bottomless, allowing us to reproduce periods of non-overflow at the Lez Spring when the mean water level in the aquifer stands below 65 m a.s.l., mainly during summer periods due to pumping in the karst conduit. Figure 8a and b show the various daily time series required for the hydrological modelling of the Lez karst system (i.e. P, ET, and Q_pump).

Figure 7 Screenshot of KarstMod, with a focus on the panel “model structure” for the Lez karst system. The solid lines correspond to the activated fluxes, whereas the faded coloured lines are not activated. Q_loss stands for the surface water discharge from the epikarst compartment, $Q_{\mathrm{p}}^{\mathrm{C}}$ stands for groundwater abstraction that affects compartment C, while Z_C stands for piezometric-head measurements considered to be representative of compartment C.

Figure 8 Daily time series for the Lez system: (a) precipitations (P) and evapotranspiration (ET), (b) groundwater abstraction, (Q_pump), (c) observed and simulated karst spring discharge (Q_Lez obs and Q_Lez sim), (d) observed and simulated piezometric head (Z_Lez obs and Z_Lez sim), (e) surface water discharge (Q_loss), and (f) simulated exchanges fluxes between compartment M and C (Q_MC).

The available hydrological observations for model calibration consist of spring discharge Q_s, piezometric-head measurements Z_c at the Lez Spring, and surface water discharge from secondary outlets and intermittent springs Q_loss (Fig. 8c, d, and e).

The surface water discharge is estimated to be the difference between discharge measured at the Lavalette station (15 km downstream of the Lez Spring) and the discharge measured at the Lez Spring, as proposed by Cousquer and Jourde (2022). Therefore, Q_loss includes all the water lost from the epikarst within several seasonally overflowing springs (i.e. Lirou Spring, Restinclière Spring, and Fleurette Spring). KarstMod allows for easy handling of the various parameter estimations depending on the considered hydrological observations (i.e. spring discharge, piezometric-head measurements, and surface discharge from the epikarst). The sampling procedure is set up to find n_obj=5000 simulations with an aggregated objective function Φ greater than 0.6. As suggested by Cousquer and Jourde (2022), using complementary hydrological observations in addition to the spring discharge allows for us to reduce the parametric uncertainties in the modelling of the Lez Spring discharge. Therefore, using a multi-objective calibration procedure implemented in KarstMod, the objective function is built as follows:

$$\begin{array}{}\text{(5)}& \mathrm{\Phi }={\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}\times \mathrm{NSE}\left(Q_{\mathrm{s}}\right)+{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}\times \mathrm{NSE}\left(Z_{\mathrm{c}}\right)+{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}\times \mathrm{NSE}\left(Q_{\mathrm{loss}}\right).\end{array}$$

The calibration procedure leads to an optimal Φ=0.65 decomposed as ϕ Q_s =0.70, ϕ Z_c =0.57, and ϕ Q_loss =0.70 within the calibration period. Model performance evaluation based on the validation period shows suitable model performance for both spring discharge and the piezometric head, with ϕ Q_s =0.54 and ϕ Z_c =0.79, but poor model performance according to the surface water discharge, with ϕ Q_loss =0.36. After this, the results can be evaluated using the various features proposed in KarstMod (Fig. 9). The results show higher model performances for Q_s and Z_c than for Q_loss. The model performance appears to be quite satisfactory with regard to the variable of interest when it comes to assessing the impact of the water management scheme on the groundwater resources within the Lez aquifer.

The simulated exchange fluxes between compartments M and C (Fig. 8f) show consistent dynamics compared to the observations. Indeed, during periods of high flow, the exchange fluxes are oriented from compartment C to compartment M (i.e. Q_MC<0). Significant precipitation events lead to rapid rises in the piezometric head; saturation of the transmissive part of the aquifer; and, finally, the establishment of overflow at the Lez Spring (i.e. Q_s>0), as well as other overflowing springs (i.e. Q_loss>0). Conversely, during the periods of low piezometric heads (i.e. both Q_s and Q_loss are nil), the simulated exchange fluxes are oriented from compartment M to compartment C (i.e. Q_MC>0). Such flow exchanges between the capacitive and transmissive parts of karst aquifers have been evidenced using KarstMod based on other karst environments (Duran et al., 2020; Frank et al., 2021; Labat et al., 2022; Sivelle et al., 2019).

Figure 9 Screenshot of KarstMod, with a focus on the sub-panel “model evaluation”. Application for the model evaluation based on the Lez system. The panel is composed such that (i) each row corresponds to the variable for calibration (Q_S, Q_loss, and Piezo_C), and (ii) each column corresponds to (a) diagnostic efficiency plot, (b) radar plots (one should note that VE and BE are not computed according to the piezometric time series), (c) decomposition of KGE (square) and KGENP (triangle), and (d) radar plot of the KGE decomposition.

5 Conclusions

KarstMod consists of a useful tool for the assessment of karst groundwater variability and sensitivity to anthropogenic pressures (e.g. groundwater abstraction). This tool is devoted to promoting good practices in hydrological modelling for learning and occasional users. KarstMod requires no programming skills and offers a user-friendly interface, allowing any user to easily manage hydrological modelling. As a first step, KarstMod can be used to explore the ability of conceptual representations to explain observations such as discharge or piezometric heads in karst systems. More advanced use of KarstMod is also possible as it provides a complete framework for (i) primary analysis of data, (ii) comparison of various model structures, (iii) evaluation of hydrological model performance, and (iv) first assessment of parametric uncertainties. The research community is increasingly using KarstMod to address various challenges in karst hydrology, from understanding hydrological processes to practical applications such as the evaluation of groundwater management plans or even the assessment of the impact of groundwater abstraction and climate change on karst groundwater resources.

Future developments of KarstMod might include (i) a consideration of the spatial heterogeneity in recharge processes, which is essential when considering snowmelt and land cover (Sivelle et al., 2022a); (ii) the simulation of electrical conductivity (Chang et al., 2021), major-ion concentration (Hartmann et al., 2013), or natural tracers such as air excess (Sivelle et al., 2022a); and (iii) the assessment of structural uncertainty (Cousquer et al., 2022). KarstMod should tend toward an open-source research software to avoid duplication of efforts in karst hydrological modelling. Also, a Python version is required for a better connection to an additional framework for sensitivity analysis, such as the SAFE toolbox (Pianosi et al., 2015), and for model calibration procedures such as particle swarm optimization (Eberhart and Kennedy, 1995; Lee, 2014). Finally, the development of the KarstMod modelling platform will benefit from better transparency and repeatability with an open-source approach, as observed with other numerical tools (Pianosi et al., 2020).

Appendix A: Abbreviations

AET	Actual evapotranspiration (mm d−1)
CFR	Refreezing coefficient (−)
CWH	Liquid-water-holding capacity of snow (−)
DE	Diagnostic efficiency (Schwemmle et al., 2021)
ET	Evapotranspiration (mm d−1)
KGE	Kling–Gupta efficiency (Gupta et al., 2009)
KGENP	Non-parametric Kling–Gupta efficiency
	(Pool et al., 2018)
MF	Melt coefficient (mm °C−1 d−1)
P	Precipitation (mm d−1)
Peff	Effective precipitation (mm d−1)
Psr	Precipitation computed with the snow routine
	(mm d−1)
$P_{\mathrm{sr}}^{\ast }$	Precipitation for a single sub-catchment computed
	with the snow routine (mm d−1)
PET	Potential evapotranspiration (mm d−1)
Rp	Pearson's correlation coefficient
Rs	Spearman rank correlation coefficient
NSE	Nash–Sutcliffe efficiency (Nash and Sutcliffe, 1970)
nobj	Targeted number of parameter sets
QA	Water discharge considered for the flow
	component A (m3 s−1)
T	Air temperature (°C)
Ta	Active temperature for snowmelt (°C)
Tn	Active temperature for refreezing (°C)
tmax	Maximum duration for sampling the
	parameter space (s)
Ts	Temperature threshold (°C)
U	Observations considered for parameter estimation
VE	Volumetric error (Criss and Winston, 2008)
ZA	Water level considered for element A (m a.s.l.)
ϕ	Performance criteria
Φ	Objective function