# On the value of hydrochemical data for the interpretation of flow and transport processes in the Baget karst system, France

Beatrice Richieri, Daniel Bittner, Vianney Sivelle, Andreas Hartmann, David Labat, Gabriele Chiogna

► To cite this version:

Beatrice Richieri, Daniel Bittner, Vianney Sivelle, Andreas Hartmann, David Labat, et al.. On the value of hydrochemical data for the interpretation of flow and transport processes in the Baget karst system, France. Hydrogeology Journal, 2024, 32 (5), pp.1537-1555. 10.1007/s10040-024-02801-2. hal-04818745

HAL Id: hal-04818745

https://hal.science/hal-04818745v1

Submitted on 5 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# On the value of hydrochemical data for the interpretation of flow and transport processes in the Baget karst system, France

Beatrice Richieri<fnref n="1" />,<fnref n="2" />,<fnref n="3" /> ID · Daniel Bittner<fnref n="3" />,<fnref n="4" /> · Vianney Sivelle<fnref n="5" /> · Andreas Hartmann<fnref n="3" /> · David Labat<fnref n="6" /> · Gabriele Chiogna<fnref n="2" />

Received: 13 September 2023 / Accepted: 16 May 2024 / Published online: 12 June 2024
© The Author(s) 2024

## Abstract

Continuous hourly time series of hydrochemical data can provide insights into the subsurface dynamics and main hydrological processes of karst systems. This study investigates how high-resolution hydrochemical data can be used for the verification of robust conceptual event-based karst models. To match the high temporal variability of hydrochemical data, the LuKARS 2.0 model was developed on an hourly scale. The model concept considers the interaction between the matrix and conduit components to allow a flexible conceptualization of binary karst systems characterized by a perennial spring and intermittent overflow as well as possible surface water bypassing the spring. The model was tested on the Baget karst system, France, featuring a recharge area defined by the coexistence of karst and nonkarst areas. The Morris screening method was used to investigate parameter sensitivity, and to calibrate the model according to the Kling-Gupta Efficiency (KGE). Model verification was performed by considering additional hydrochemical constraints with the aim of representing the internal dynamics of the systems, i.e., water contributions from the various compartments of the conceptual model. The hydrochemical constraints were defined based on high-temporal resolution time series of $SO_4^{2-}$ and $HCO_3^-$.

The results of this study show that the simulation with the highest KGE among 9,000 model realizations well represents the dynamics of the spring discharge but not the variability of the internal fluxes. The implementation of hydrochemical constraints facilitates the identification of realizations reproducing the observed relative increase in the flow contribution from the nonkarst area.

**Keywords** Karst · Conceptual models · Rainfall-runoff · Hydrochemistry · France

## Introduction

Karst systems cover 10–15% of the Earth’s surface and 35% of the European continental area (Goldscheider et al. 2020), representing a major worldwide source of freshwater providing drinking water to 10–25% of the world’s population (Ford and Williams 2007; Stevanović 2019). Hydrological models support the understanding of the system’s functioning and are fundamental to ensuring the sustainable water management of karst water resources (Hartmann et al. 2014). However, modelling of karst systems is still a difficult task due to the complex interaction between the matrix and conduit domains (Hartmann et al. 2013) and the difficulty in observing and measuring flow and transport processes in the subsurface (Berthelin and Hartmann 2020). Sivelle et al. (2022b) investigated the relevance of excess air and, thus, hydrostatic pressure for the assessment of matrix-conduit exchange. As a result, large uncertainty characterizes the internal dynamics of karst systems, i.e., fast discharge

<figure><img src="image_2.png" /></figure>

<footnote>Beatrice Richieri
beatrice.richieri@gmail.com</footnote>

<footnote n="1">¹ TUM School of Engineering and Design, Technical University of Munich, Arcisstr. 21, 80333 Munich, Germany</footnote>

<footnote n="2">² GeoZentrum Nordbayern, Department Geography and Geosciences, Friedrich-Alexander-Universität Erlangen-Nürnberg (FAU), Schlossgarten 5, 91054 Erlangen, Germany</footnote>

<footnote n="3">³ Institute of Groundwater Management (IGW), Technical University of Dresden, Bergstr. 66, 01069 Dresden, Germany</footnote>

<footnote n="4">⁴ Department for River Basin Management, Erftverband, Am Erftverband 6, 50126 Bergheim, Germany</footnote>

<footnote n="5">⁵ National Research Council, Research Institute for Geo-Hydrological Protection, Perugia, Italy</footnote>

<footnote n="6">⁶ Université Toulouse 3-Géosciences Environnement Toulouse-CNRS-UPS-IRD, 14 Avenue Edouard Belin, 31400 Toulouse, France</footnote>
through the conduit, infiltration into the matrix and water contributions from different geological formations inside the recharge area (Chang et al. 2021).

Hence, hydrological models need to be calibrated due to the lack of knowledge in model parameter values and process understanding (Le Moine et al. 2008). Moreover, model validation is necessary to assess model reliability and utility for further applications (Andréassian 2023; Klemeš 1986). Typically, model calibration and validation are done using measured discharge time series at a monitoring site and computing performance metrics. Among the several existing metrics (Bennett et al. 2013; Ferreira et al. 2020; Moriasi et al. 2007), the Kling-Gupta Efficiency (KGE) is considered suitable for capturing the entire flow regime (Gupta et al. 2009). However, a single metric is generally unable to properly evaluate all model characteristics and solely relying on optimal values of an objective function to assess a model has been criticized (Gupta et al. 2008; Leins et al. 2023). Indeed, acceptable values of a model performance indicator do not necessarily mean that the hydrological model is reliable due to the problem of equifinality (Cinkus et al. 2022), i.e., the existence of multiple optimal parameter sets that reproduce the observed values (Chiogna et al. 2024). Different models, characterized by different structures or parameter values, can result in satisfactory simulated spring discharge (Mudarra et al. 2019). Previous works have shown that models validated with a single spring monitoring site may not capture the internal functioning of karst system and result in large equifinality (Hartmann et al. 2017).

The use of hydrochemical data allows for a better understanding of the internal dynamics of karst systems. Besides the use of conservative tracers such as stable water isotopes (Wang et al. 2021; Winston and Criss 2004), major ions and electrical conductivity (EC) are also used to investigate the spatial and temporal variability in the hydrological response of karts systems (Barbieri et al. 2005; Chang et al. 2021; Gil-Márquez et al. 2017; Hartmann et al. 2013). However, to take full advantage of the benefits of hydrochemical data to better capture the variability of the internal response of a system, the hydrochemical data should be at a resolution at least as detailed as the resolution of discharge observations, in particular for event-based models. Nevertheless, hydrochemical data are typically collected at weekly or coarser temporal resolutions due to high analysis costs and time constraints.

The aim of this study is to further contribute to ongoing research into the coupling of hydrological (e.g., spring discharge) and hydrochemical (e.g., high-resolution ion concentration time series) information for the verification of conceptual event-based models. The method proposed by Richieri et al. (2023) is used to retrieve high-resolution hydrochemical information from continuous measurements of electrical conductivity (EC). Moreover, the semidistributed LuKARS 2.0 (land use change modelling in KARSt systems) model, based on the original model from Bittner et al. (2018), is further developed by using it for the first time at an hourly timescale.

LuKARS represents the discharge observed at the spring as the sum of water contributions draining different geological areas (hydrotopes) characterized by different response times. In this work, this approach is validated considering that different parts of the catchment can be characterized by different chemical signatures. The new version of the model represents the different compartments of a karst system as buckets, i.e., epikarst, matrix and conduits, and aims to provide an adaptable representation of binary karst systems with a perennial spring and intermittent overflow. In particular, the model provides an updated conceptualization of the interaction between the matrix and conduit in comparison to Bittner et al. (2018) to allow a more flexible and consistent conceptualization of the subsurface system following the approach proposed in KarstMod (Mazzilli et al. 2019; Sivelle et al. 2023). To constrain model parameters, the KGE performance criteria are combined with the conceptual information derived from high-resolution hydrochemical data, i.e., major ion concentrations. The hypothesis of the present work is that high-temporal-resolution hydrochemical data can support the verification of the conceptual model, whereas model structures simulating the spring discharge without consideration of hydrochemistry may feature a comparable KGE but may not match the internal dynamics of a system.

The model was tested on the Baget system, an intensively studied karst watershed in France (Labat et al. 1999; Sivelle et al. 2022a) which is characterized by an absence of inter-basin groundwater flow (Mangin 1975; AL Khoury et al. 2023). Baget shows a heterogeneous geology combining limestone bedrock and an extensive outcrop formation of black flysch (karst and non-karst areas) (Ulloa-Cedamanos et al. 2020). Time series at hourly resolution of the contributions of $HCO_3^-$ and $SO_4^{2-}$ to the total EC, called weight factors and derived by Richieri et al. (2023), contribute to the description of the internal dynamics of the Baget system, i.e., the response of water contributions from different geological areas in the catchment to different flow conditions.

The manuscript is structured as follows. First, the study area and the available dataset are presented. Then the new features of LuKARST 2.0 are described in detail, whereas the components of the model that remain unchanged are provided in the electronic supplementary material (ESM). Finally, the Morris method is used for the sensitivity analysis of the model parameters and an envelope of behavioral simulations is defined based on KGE. Finally, this study shows that only a limited subset of behavioral simulations can explain the observed hydrochemical dynamics.

## Materials and methods

This section provides information about the Baget karst watershed and the hydrochemical dataset used to select an appropriate conceptual model and parameter values.

<figure><img src="image_3.png" /></figure>
The hydrochemical data were collected at the Las Hountas spring for the period 30/03/2022–7/04/2022. This section also describes the previous investigations into the internal dynamics of the Baget karst system based on hydrochemical data at high temporal resolution (Richieri et al. 2023).

## Study area

The Baget karst system is a binary karst system characterized by a perennial spring and intermittent overflow, located 10 km southwest of the city of Saint-Girons, in the Ariege administrative department, France (Fig. 1a,b). It has a recharge area of ~13 km² and does not receive water from adjacent catchments (Mangin 1975; AL Khoury et al. 2023). The vegetation of the recharge area is dominated by a

fir-beech forest, with only a small agricultural plot covering 3% of the recharge area (Sivelle et al. 2022a). Vegetation is particularly dense along the north-facing slope, whereas grassland occupies some of the south-facing slope. The Baget system is under the influence of the Atlantic oceanic climate with a mean annual air temperature of 12.3 °C, an average annual rainfall close to 1,700 mm, and no influence of snow melt processes (Padilla et al. 1994). The annual precipitation distribution is bimodal, with peaks occurring in December and February (Ulloa-Cedamanos et al. 2020).

The catchment features heterogeneous geology, dominated by calcareous lithologies originating from the Jurassic and Lower Cretaceous. In addition to the limestone bedrock, a large outcrop of relatively impermeable black flysch partially covers the south-facing slope (Fig. 1c–e). The extent of the black flysch was obtained using the BD Charm-50

<figure><img src="image_4.png" /><figcaption>**Fig. 1** Overview of the study site. **a** The location of the study area in Europe. **b** The location of the Antichan rain gauging station in relation to the study area. **c** Geological map of the Baget catchment with the location of the Las Hountas spring (modified after the BD Charm-50 geology map from the French Geological Survey). **d** Zoomed in area of the geological map showing the locations of the Las Hountas spring, at which the water samples were collected, and the outlet of the recharge area. **e** Geological subsurface cross-section in the location indicated with a black line (d). Elevation is provided in meters above mean sea level (m AMSL) (modified from Debroas 2009)</figcaption></figure>
geology map (French Geological Survey) by measuring the corresponding covered area. Its estimation of ~30% of the catchment is comparable to the 25.2% estimated by Ulloa-Cedamanos et al. (2020). Spring flow contributions derived from different geological areas of the catchment have undergone different dissolution processes. Calcium carbonate (CaCO₃) associated with the limestone is dissolved by carbonic acid (H₂CO₃) resulting in water typified by elevated Ca²⁺ and HCO₃⁻. The black flysch contains pyrite, where oxidation processes release sulfuric acid (H₂SO₄), which in turn dissolves CaCO₃ with the products SO₄²⁻, Ca²⁺ and HCO₃⁻. Consequently, the water contributions from the limestone bedrock and black flysch provide different components to the chemical signature at the spring (Richieri et al. 2023). A more detailed overview of the different sources affecting the chemical composition at the spring is given in section 'Prior investigation into hydrochemical signals'.

The karstic watershed is characterized by one perennial spring, called Las Hountas. Las Hountas is situated 110 m away from the outlet of the catchment (Fig. 1d) and is representative of a part of the total response of the system. During precipitation events, the downstream part of the Baget catchment recharges the Las Hountas spring. Only during high flood events the karst conduits in the upper part of the catchment start to discharge water to the surface, actively contributing to the streamflow generation in the Lachein stream, which is usually dry (Sivelle et al. 2020). Groundwater discharge to streamflow in the upper part of the surface channel occurs ~50 days/year (Mangin 1975). The surface water of the Lachein stream bypasses the Las Hountas spring and directly reaches the outlet of the catchment (Ulloa-Cedamanos et al. 2021). This is evident by comparing the measured discharge at the outlet and at Las Hountas spring: the discharge at the spring shows a plateau at ~0.6 m³/s (Fig. S1 in the ESM).

## Data collection

The hourly precipitation data used as an input for the model were recorded at the meteorological station of Antichan, ~8 km away from the spring (Fig. 1b). The observed discharge time series at the Las Hountas spring was derived by recording water level with an Aqua TROLL 200 device (In-Situ Inc., United States) and applying the rating curve from Mangin (1975), which was adjusted to represent the current state of the cross section during additional field measurements performed in 2021 and 2022. Evapotranspiration was not considered because this study focuses on the simulation of peak flows in response to intense precipitation events during which evapotranspiration effects are negligible.

Table 1 contains information about the temporal resolution, number of samples and statistics of the specific EC (µS/

<table><caption>**Table 1** Statistics of the used datasets, including temporal resolution (TR), number of samples (n), mean, standard deviation (SD), minimum (MIN) and maximum (MAX) values of the specific electrical conductivity EC (µS/cm), major solute species (mg/L), pH and water temperature T (°C) measured at the Las Hountas spring for the events 30/03/2022–7/04/2022</caption><thead><tr><th>Statistic</th><th>EC (µS/cm)</th><th>Ca (mg/L)</th><th>Mg (mg/L)</th><th>Na (mg/L)</th><th>K (mg/L)</th><th>HCO<sub>3</sub><sup>-</sup> (mg/L)</th><th>SO<sub>4</sub> (mg/L)</th><th>NO<sub>3</sub><sup>-</sup> (mg/L)</th><th>Cl (mg/L)</th><th>pH</th><th>T (°C)</th></tr></thead><tbody><tr><td>TR</td><td>Quarter-hour</td><td>Event-based</td><td>Event-based</td><td>Event-based</td><td>Event-based</td><td>Event-based</td><td>Event-based</td><td>Event-based</td><td>Event-based</td><td>Event-based</td><td>Quarter-hour</td></tr><tr><td>n</td><td>830</td><td>122</td><td>122</td><td>122</td><td>122</td><td>122</td><td>122</td><td>122</td><td>122</td><td>122</td><td>830</td></tr><tr><td>Mean</td><td>332.2</td><td>64.5</td><td>4.5</td><td>1.1</td><td>0.6</td><td>193.6</td><td>10.8</td><td>1.6</td><td>1.5</td><td>8.2</td><td>9.9</td></tr><tr><td>SD</td><td>7.6</td><td>1.9</td><td>0.3</td><td>0.1</td><td>0.1</td><td>4.4</td><td>2.5</td><td>0.2</td><td>0.08</td><td>0.2</td><td>0.2</td></tr><tr><td>MIN</td><td>316.5</td><td>61.1</td><td>3.7</td><td>0.9</td><td>0.3</td><td>185.5</td><td>7.9</td><td>1.2</td><td>1.3</td><td>7.2</td><td>5.7</td></tr><tr><td>MAX</td><td>358.2</td><td>68.9</td><td>5.2</td><td>1.4</td><td>1.0</td><td>205.6</td><td>19.4</td><td>1.9</td><td>1.9</td><td>8.3</td><td>10.2</td></tr></tbody></table>

<figure><img src="image_5.png" /></figure>
cm), water temperature, pH and major solute species total concentration (i.e., each present both as free ion and as part of complexes as described in Richieri et al. 2023), i.e., Ca, Mg, HCO₃, SO₄, NO₃, Cl, Na and K (mg/L), which were measured at the Las Hountas spring during the event-based sampling campaign from 30/03/2022 to 7/04/2022. The specific EC was measured with a time interval of 15 min by means of the In-Situ Aqua TROLL 200 device and reported at the standard temperature of 25 °C. The water samples were collected by a 6712 ISCO sampler (Teledyne ISCO, United States), which was connected to the EC probe to automatically start sampling above a water level threshold (30 cm). The sampling frequency was hourly during the rising curve of the hydrograph, every 2 h during the recession curve, and then every 3 h near baseflow conditions. The ISCO sampler was installed inside a shelter that always provided shade. The samples were collected every day in plastic vials with zero headspace, filtered through a 0.22-µm membrane and stored in a refrigerator at 4 °C for ~1 week prior to analysis. For each sample, an aliquot for cation analysis was acidified with nitric acid (HNO₃) to prevent complexation and precipitation (Weiss 2020; Ulloa-Cedamanos et al. 2020). The solute concentrations, provided at the standard temperature of 25 °C, were analyzed by the laboratory Geosciences Environment Toulouse, France. The ICP-OES (inductively coupled plasma optical emission spectrometry) was used to quantify Ca, Mg, Na and K; ion chromatography was used to quantify NO₃, SO₄ and Cl; titration analysis was used to quantify HCO₃⁻.

To verify the accuracy of chemical quantification, the charge balance of each sample was computed using the software PHREEQC. The charge mass balance error ranged from 0.8 to 4.9% with a median of 2.4%. Since all the samples show a charge error lower than ±5%, the results of the laboratory analysis were considered reliable (Parkhurst and Appelo 2013). Figure S2 in the ESM shows the time series of the computed charge balance (%) together with the sum of cations (TC) and anions (TA) (mEq/L).

## Prior investigation into hydrochemical signals

Richieri et al. (2023) previously investigated the hydrological functioning of the Baget karst system by comparing continuous water level recordings and EC measurements with high-temporal resolution major solute concentrations for multiple precipitation events, occurring in October 2021, November 2021 and November 2022. EC dynamics confirmed the complex hydrological response behavior of the Las Hountas spring, which is characterized by a simultaneous increase in water level and EC (flushing and piston effects) during precipitation events that followed dry periods, and by dilution processes during peak spring discharge periods. Richieri et al. (2023) also used the individual

contributions of individual major ions i to the total observed EC, called weight factors $f_i$, to identify the varying contributions of water derived from different areas of the watershed, i.e., limestone and black flysch, at different flow conditions. The individual contributions of each ion to the total EC were computed for each sample with PHREEQC, by considering ion molar conductivity, molar concentration and the electrochemical activity coefficient as well as pH and temperature of the sample (Table 1). A complete description of the equations used to compute the weight factors is provided in Richieri et al. (2023). Equation (1) describes the dissolution process characterizing the limestone formation. The dissolution of one mole of calcium carbonate (CaCO₃) by carbonic acid (H₂CO₃) produces one mole of Ca²⁺ and two moles of HCO₃⁻. Equation (2) describes the dissolution process characterizing the black flysch containing pyrite, whose oxidation releases strong acids, i.e., sulfuric acid (H₂SO₄). The reaction between one mole of CaCO₃ and half a mole of H₂SO₄ produces one mole of Ca²⁺, one mole of HCO₃⁻ and half a mole of SO₄²⁻. Thus, the dissolution of CaCO₃ by both H₂CO₃ and H₂SO₄ leads to a lower alkalinity (HCO₃⁻) than what would be observed in the case of only dissolution by H₂CO₃. Consequently, it was shown that the increase and decrease in the weight factors of SO₄²⁻ and HCO₃⁻ ($f_{SO₄²⁻}$ and $f_{HCO₃⁻}$), respectively, indicate a relative increase in the water draining the black flysch, which was observed to be simultaneous to flushing and piston effects. On the contrary, the increase in $f_{HCO₃⁻}$ indicates a larger water contribution from the limestone bedrocks during dilution processes and baseflow conditions (Richieri et al. 2023).

$$CaCO_3 + H_2CO_3 \rightarrow Ca^{2+} + 2HCO_3^- \quad (1)$$

$$CaCO_3 + 0.5H_2SO_4 \rightarrow Ca^{2+} + HCO_3^- + 0.5SO_4^{2-} \quad (2)$$

It is pertinent to note, that despite the black flysch containing Na silicate minerals (Na₂SiO₃) (Ulloa-Cedamanos et al. 2021), the origin of HCO₃⁻ is considered to be dominated by CaCO₃ dissolution due to the low Na⁺ content of water samples (Table 1). The variability in silicate weathering input could be further investigated by means of isotopes (Hagedorn and Whittier 2015; Spence and Telmer 2005), which were not available for the present study. The role of Na silicate weathering by H₂CO₃ (Ulloa-Cedamanos et al. 2021) is therefore considered negligible in the Baget catchment. In addition, a minor contribution of SO₄²⁻ might be related to a gypsum formation within the catchment. However, as the area of gypsum bedrock only represents 0.2% of the recharge area (Ulloa-Cedamanos et al. 2020), the present study considers the mass flux of SO₄²⁻ to be controlled by the discharge from the black flysch formation, which covers 30% of the catchment.
Despite being often considered a conservative tracer, Cl⁻ is not used in this study for the investigation of system dynamics for two reasons. Firstly, the computation of the weight factors fᵢ at high temporal resolution is affected by high uncertainty in the case of ions with low concentration such as Cl⁻ (Table 1; Richieri et al. 2023). Second, this work focuses on the investigation of water draining different geological formations in the catchment, which does not contain any halite deposit or other geogenic sources of Cl⁻ (Ulloa-Cedamanos et al. 2021).

The Baget catchment is relatively unpolluted with an absence of SO₄²⁻, NO₃⁻ and Cl⁻ originating from anthropogenic activities. NO₃⁻ was observed to be associated with organic decomposition and to increase in autumn. Indeed, the seasonal increase in precipitation causes the leaching of NO₃⁻ from the soil (Ulloa-Cedamanos et al. 2020), whereas the observed SO₄²⁻ and Cl⁻ are entirely derived from geogenic sources and precipitation, respectively (Ulloa-Cedamanos et al. 2021).

## Model development

This section describes the LuKARS 2.0 model concept, together with the modifications done within this study with respect to the original model from Bittner et al. (2018). Subsequently, the Morris screening sensitivity analysis and the calibration procedure are presented. Finally, the use of high-resolution hydrochemical data for the selection of the model structure and parameters is described. The model was tested for the Las Hountas spring and calibrated

and validated for the periods 1/03/2022–29/03/2022 and 30/03/2022–30/04/2022, respectively.

## Description of LuKARS 2.0 model concept

LuKARS is a semidistributed model developed by Bittner et al. (2018). The model divides the catchment into hydrotopes, which are defined as independent units exhibiting similar hydrological behavior and soil characteristics. The original model from Bittner et al. (2018) was later modified by Sivelle et al. (2022a) to be coupled with KarstMod (Mazzilli et al. 2019). In this study, LuKARS 2.0 was developed and applied to the Baget catchment starting from the original model (Bittner et al. 2018) considering two hydrotopes in the upper compartment, which represent the main geological formations within the studied catchment, i.e., limestone and black flysch. LuKARS 2.0 adds three new features to LuKARS: (1) the model was modified from a daily to an hourly time step and the parameters adjusted accordingly; (2) a transfer between the matrix and conduit was implemented in the lower compartment of the model; and (3) a drainage from the conduit was implemented to represent water bypassing the spring at high flow conditions (Fig. S1 in the ESM). Figure 2 and Table 2 show the model concept and provide a description of the model parameters, respectively. The model concept was developed based on the understanding of system functioning, more precisely on the dynamics of the water contributions draining the different geological formations (Richieri et al. 2023). Since the hydrological response of the contribution from the black flysch was observed to be fast with no apparent influence on the

<figure><img src="image_6.png" /><figcaption>**Fig. 2** LuKARS 2.0 model concept including the implemented transfer between the matrix bucket and the conduit bucket and the drainage from the conduit bucket for the case of the Las Hountas spring. The</figcaption></figure>

two hydrotopes are defined based on the main geological formations present in the recharge area, i.e., karst bedrock and black flysch, and their response to rainfall

<figure><img src="image_7.png" /></figure>
<table><caption>**Table 2** Overview of the model parameters of LuKARS 2.0, including *symbol*, *unit*, *description*, and *range* used for the Morris analysis</caption><thead><tr><th>Area</th><th>Symbol</th><th>Unit</th><th>Description</th><th>Range</th></tr></thead><tbody><tr><td>Catchment</td><td>Ra</td><td>(km<sup>2</sup>)</td><td>Total recharge area</td><td>9–15</td></tr><tr><td rowspan="6">Karst (HYD1)</td><td>l<sub>hyd_1</sub></td><td>(m)</td><td>Mean distance of hydrotepe 1 to the spring</td><td>1,500–6,500</td></tr><tr><td>k<sub>hyd_1</sub></td><td>(mm<sup>2</sup>/h)</td><td>Discharge coefficient for the fast flow from hydrotepe 1 to the conduit (Q<sub>hyd_1</sub>)</td><td>1–10,000</td></tr><tr><td>E<sub>min_1</sub></td><td>(mm)</td><td>Lower storage threshold for hydrotepe 1</td><td>8–20</td></tr><tr><td>E<sub>max_1</sub></td><td>(mm)</td><td>Upper storage threshold for hydrotepe 1</td><td>21–200</td></tr><tr><td>α<sub>_1</sub></td><td>(-)</td><td>Exponent for the fast flow from hydrotepe 1 to the conduit (Q<sub>hyd_1</sub>)</td><td>0–2</td></tr><tr><td>k<sub>is_1</sub></td><td>(1/h)</td><td>Discharge coefficient for the infiltration from hydrotepe 1 to the matrix (Q<sub>is_1</sub>)</td><td>10<sup>-7</sup>–10<sup>-3</sup></td></tr><tr><td rowspan="6">Black flysch (HYD2)</td><td>l<sub>hyd_2</sub></td><td>(m)</td><td>Mean distance of hydrotepe 2 to the spring</td><td>1,000–4,000</td></tr><tr><td>k<sub>hyd_2</sub></td><td>(mm<sup>2</sup>/h)</td><td>Discharge coefficient for the fast flow from hydrotepe 2 to the conduit (Q<sub>hyd_2</sub>)</td><td>1–1,000</td></tr><tr><td>E<sub>min_2</sub></td><td>(mm)</td><td>Lower storage threshold for hydrotepe 2</td><td>1–5</td></tr><tr><td>E<sub>max_2</sub></td><td>(mm)</td><td>Upper storage threshold for hydrotepe 2</td><td>6–30</td></tr><tr><td>α<sub>_2</sub></td><td>(-)</td><td>Exponent for the fast flow from hydrotepe 2 to the conduit (Q<sub>hyd_2</sub>)</td><td>0–2</td></tr><tr><td>k<sub>MC</sub></td><td>(mm/h)</td><td>Discharge coefficient for the transfer between matric and conduit (Q<sub>MC</sub>)</td><td>10<sup>-5</sup>–10<sup>-1</sup></td></tr><tr><td rowspan="5">Lower compartment</td><td>a<sub>MC</sub></td><td>(-)</td><td>Exponent for the transfer between matric and conduit (Q<sub>MC</sub>)</td><td>1–3</td></tr><tr><td>C<sub>loss</sub></td><td>(mm)</td><td>Storage threshold for the conduit</td><td>10<sup>-2</sup>–10<sup>0</sup></td></tr><tr><td>k<sub>CS</sub></td><td>(mm/h)</td><td>Discharge coefficient for the flow from the conduit to the spring (Q<sub>CS</sub>)</td><td>10<sup>-3</sup>–10<sup>1</sup></td></tr><tr><td>a<sub>CS</sub></td><td>(-)</td><td>Exponent for the flow from the conduit to the spring (Q<sub>CS</sub>)</td><td>1–4</td></tr></tbody></table>

baseflow, hydrotepe 2 was defined as such to provide only fast flow to the conduit but no infiltration to the matrix. On the other hand, hydrotepe 1 represents the limestone karst formation, which contributes to both the discharge to the conduit and infiltration in the matrix. The fast flow $Q_{hyd}$ (m³/s) from a certain hydrotepe gets active when the water level in that hydrotepe reaches the upper storage threshold $E_{max}$ (mm) and continues until the water level goes down to the lower storage threshold $E_{min}$ (mm) (hysteresis). Nevertheless, the infiltration $Q_{is}$ (m³/s) is always active and linearly correlated with the water level in the hydrotepe. Finally, the water in the conduit will be transferred to the spring with the linear function from Mazzilli et al. (2019) (ESM, Eq. S6). The implemented modifications to LuKARS are described here, whereas the model equations from the original model (Bittner et al. 2018) are reported in the ESM.

*Transfer between the matrix and conduit.* A transfer between the matrix and conduit $Q_{MC}$ (m³/s) was implemented in the lower compartment of the model to represent the dual behavior typical of karst systems. The transfer was defined based on the approach of Mazzilli et al. (2023) as a function of the recharge area $Ra$ (km²) discharge coefficient $k_{MC}$ (mm/h), exponent $a_{MC}$ (-) and dimensionless water levels in the matrix $M$ and conduit $C$ (-). $Q_{MC}$ (m³/s) was implemented as shown in Eq. (3), with $t$ indicating the current time step, abs the absolute value and sgn the sign of the subtraction of the dimensionless water level in the conduit $C$ from that in the matrix $M$ (-). Positive values of $Q_{MC}$ (m³/s) means that the current direction of flow is from the matrix to the conduit. To avoid numerical instabilities, Eq. (3) is solved using the analytical solution for the inter-compartment coupling (Mazzilli et al. 2023).

$$Q_{MC,t} = Ra \times k_{MC} \times \text{sgn}(M_t - C_t) \times \text{abs}(M_t - C_t)^{a_{MC}} \quad (3)$$

*Discharge from the conduit.* A drainage $Q_{\text{Loss}}$ (m³/s) from the conduit was used to represent the plateau in discharge and the water bypassing the spring (Fig. S1 of the ESM), which results from the natural drainage through the conduits at high flow conditions. The drainage was implemented to start when the water level in the conduit $C_t$ (where the subscript $t$ indicates the time step at which a quantity is computed) gets higher than the storage threshold $C_{\text{loss}}$ (mm) and to remain active as long as the water level does not drop below it. Equation (4) shows the implementation of the drainage $Q_{\text{Loss},t}$ (m³/s), with $Ra$ (km²) the recharge area and $dt$ (h) the hourly time step. The water level is checked at each time step $t$ and the excess water exits immediately from the system without reaching the spring.

$$ \begin{array}{l} \text{if } C_t > C_{\text{loss}} : Q_{\text{Loss},t} = (C_t - C_{\text{loss}}) * \frac{Ra}{dt} \text{ and } C_t = C_{\text{loss}}, \\ \text{else : } Q_{\text{Loss},t} = 0 \end{array} \qquad (4) $$

## Sensitivity analysis (Morris screening)

The identification of the most sensitive model parameters and potential nonlinear interactions was done with the elementary effects method, also known as Morris screening (Campolongo et al. 2007; Morris 1991). Morris screening is a one-at-a-time (OAT) method to perform global sensitivity analysis which is widely applied in literature and it is often used in the case of models characterized by a
large number of input parameters for which variance-based sensitivity analysis would be computationally demanding (Campolongo et al. 2007; Jaxa-Rozen and Kwakkel 2018; Merchán-Rivera et al. 2022; Smith 2013). The Morris method considers that each model parameter set $q = [q_1, ..., q_n]$ varies across a discrete number $l$ of values, called levels, forming an $n$-dimensional $l$-level grid $\Gamma_l$. The method relies on the average of the elementary effects over parameter space to provide a measure of global sensitivity. The elementary effect $d_i(q)$ of the $i$th input parameter quantifies the approximate local sensitivity at the point $q$ and is defined as in Eq. (5) (Smith 2013)

$$d_i(q) = \frac{f(q_1, \dots, q_{i-1}, q_i, q_{i+1}, \dots, q_n) - f(q)}{\Delta} = \frac{f(q + \Delta e_i) - f(q)}{\Delta} \quad (5)$$

where $e_i$ is a vector of zeros with one on the $i$th components and $\Delta$ is the stepsize, which is chosen from the set shown in Eq. (6) (Smith 2013)

$$\Delta \in \left\{ \frac{1}{l-1}, \dots, 1 - \frac{1}{l-1} \right\} \quad (6)$$

To obtain a global sensitivity measure, Campolongo et al. (2007) and Morris (1991) proposed the mean $\mu^*$ and variance $\sigma$, respectively, of the finite-dimensional distribution $G_i$ associated with the absolute value of the elementary effect $d_i(q)$, which is derived from randomly sampling $q$ within $\Gamma_l$. Considering $r$ sampling points, these metrics associated with the $i$th parameter and the $j$th sample can be expressed as shown in Eqs. (7) and (8).

$$\mu_i^* = \frac{1}{r} \sum_{j=1}^{r} |d_i^j(q)| \quad (7)$$

$$\sigma_i^2 = \frac{1}{r-1} \sum_{j=1}^{r} (d_i^j(q) - \mu_i)^2 \quad \text{with } \mu_i = \frac{1}{r} \sum_{j=1}^{r} d_i^j(q) \quad (8)$$

where $d_i^j$ is the elementary effect associated with the $i$th parameter and $j$th sample.

High values of $\sigma$ indicate possible interactions between the model parameters and/or of parameter nonlinearity. The parameters with $\sigma/\mu^*$ smaller than 0.1 or between 0.1 and 1 are almost linear or monotonic, respectively. The parameters with $\sigma/\mu^*$ larger than 1 are characterized by marked nonmonotonic nonlinearities or interactions with other parameters (Sanchez et al. 2014). $\mu^*$ and $\sigma$ are constructed by considering $m$ trajectories of $n+1$ points in the parameter space. The total number of realizations accounted for in the metrics calculations is defined by $(n+1)*m$. For a detailed explanation of the method, one can refer to the works of Morris (1991) and Campolongo et al. (2007) as well as to that of Smith (2013).

In this study, Morris was run for the calibration period (1/03/2022–30/03/2022) by using the sensitivity analysis library in Python (SALib) and considering the 17 model parameters reported in Table 2, 500 trajectories and 100 grid levels, for a total of 9,000 realizations. The length of the calibration period was selected considering the aim of building an event-based model. The number of realizations was considered sufficient since the method with 1,000 trajectories and 200 grid levels showed convergence of the results. Table 2 shows the parameter ranges used for the Morris analysis. Since it was the first application of LuKARS at an hourly scale, it was required to initially investigate the parameter ranges by means of an independent set of 10,000 Monte Carlo simulations. In addition, the parameter ranges for the two hydrotopes were selected to be consistent with the observations of the water contribution from the two different geological formations. A lower upper storage threshold was set for the black flysch, i.e., hydrotepe 2, to represent its faster activation and, thus, the relative increase in the water contribution at the beginning of heavy precipitation events (Richieri et al. 2023). The extents of hydrotepe 1 and hydrotepe 2 were taken fixed and equal to 70 and 30% of the entire catchment, respectively. Each realization performed for the Morris method was evaluated by computing the KGE. Equation (9) shows the equation for the computation of KGE, with $r_c$ the linear correlation between the observations and simulations, $\sigma_{sim}$ the standard deviation in simulations, $\sigma_{obs}$ the standard deviation in observations, $\mu_{sim}$ the simulation mean and $\mu_{obs}$ the observation mean (Knoben et al. 2019).

$$KGE = 1 - \sqrt{(r_c - 1)^2 + \left(\frac{\sigma_{sim}}{\sigma_{obs}} - 1\right)^2 + \left(\frac{\mu_{sim}}{\mu_{obs}} - 1\right)^2} \quad (9)$$

Among the 9,000 model realizations used for the sensitivity analysis, the realization with the highest KGE was considered as the selected one, while the subset of realizations with a KGE larger than 0.5 was defined as behavioral (e.g., Beven and Freer 2001). To assess the uncertainty of the model, each flow component—$Q_{hyd_1}$, $Q_{hyd_2}$, $Q_{is_1}$, $Q_{MC}$, $Q_{Closs}$ and $Q_{CS}$ (m$^3/s$)—derived by using the selected parameters values was compared to the distribution of the behavioral simulations, i.e., interquartile and 10–90 percentile envelopes.

## Molel selection considering hydrochemical constraints

Once the selected and behavioral parameter sets were found, they were used to validate the model for the period 30/03/2022–30/04/2022. Thereafter, it was necessary to check whether the selected simulation represented not
only the discharge at the spring but also the internal
dynamics of the system, i.e., water contributions from the
limestone and black flysch. For this purpose, the flows
from hydrotope 1 (limestone) and 2 (black flysch) to the
conduit together with the transfer between the matrix and
conduit were compared to the observed hydrochemical
data recorded during the increase in the water level—
occurred in the 2 days 30/03/2022 and 1/04/2022—for the
event 30/03/2022–7/04/2022. Therefore, three additional
criteria were defined to constrain, among the multiple
Morris realizations with KGE larger than 0.5, those also
matching the hydrochemical information. The criteria
were identified based on the weight factors (Richieri
et al. 2023) computed for the time series of HCO₃⁻ and
SO₄²⁻ observed during the event in April 2022 (Table 1)
and are described in detail within ‘Model selection using
hydrochemical constraints’. It is important to note that
the weight factors were computed using PHREEQC by
considering for each water sample the measured pH
(Table 1). In addition, despite HCO₃⁻ not being a con-
servative species, the alkalinity in the case of pH between
7.2 and 8.3 (Table 1) is typically considered to derive
from HCO₃⁻ alone (Boyd 2020). The hydrochemical
constraints were defined based on HCO₃⁻ and SO₄²⁻ as
these ions were previously used for the investigation of
the internal dynamics of the Baget system (Richieri et al.
2023). In addition, both HCO₃⁻ and SO₄²⁻ were proved
to be present at Las Hountas spring at sufficiently high
concentration to be retrieved at high temporal resolution
by applying the EC decomposition method from Richieri

et al. (2023). Indeed, the EC decomposition method
(Richieri et al. 2023) allows for the retrieval of accurate
concentration time series at high temporal resolution for
those solutes present with a sufficiently high concentra-
tion, thus significantly contributing (~>6%) to the total
EC, whereas the reconstructed time series of solutes with
low concentration, e.g., Cl⁻ (Table 1), are subject to large
uncertainty.

Results

This section presents the results of the investigation. Firstly,
the results of the Morris analysis are presented in terms
of sensitivity and nonlinear interaction among the model
parameters. Secondly, the results of the model calibration
and validation considering the KGE of the spring discharge
are shown together with the respective uncertainty bands.
Finally, it is shown how the hydrochemical criteria were
used to further constrain the model by considering not only
the discharge at the spring but also the expected internal
flow dynamics.

Sensitivity analysis

The results of the sensitivity analysis based on the Mor-
ris screening method are evaluated by comparing the mean
μ* and standard deviation σ of the distribution function of
each parameter. Figure 3a shows the value of μ* for each
parameter in descending order, while Table 3 contains the

<caption>**Fig. 3** Results of the Morris screening sensitivity analysis.</caption>

a Mean of the absolute values of the elementary effects (μ*) for each model parameter as the index of their sensitivity.
b Mean of the distribution of the absolute values (μ*) against the standard deviation of the distribution (σ) for each model parameter; the areas of the graph corresponding to monotonic and/or linear interactions and nonlinear interactions (Sanchez et al. 2014) are shown with a dashed gray line

<figure><img src="image_8.png" /></figure>
<table><caption>**Table 3** Overview of the results of the Morris analysis, including for each *parameter* the *sensitivity ranking*, the *selected parameter value* and the *median and standard deviation for the behavioral simulations* based on KGE</caption><thead><tr><th>Area</th><th>Parameter</th><th>Sensitivity ranking</th><th>Selected parameter value</th><th>Behavioral median</th><th>Behavioral standard deviation</th></tr></thead><tbody><tr><td>Catchment</td><td>Ra (km²)</td><td>8</td><td>12.9</td><td>12.0</td><td>1.78</td></tr><tr><td rowspan="5">Karst (HYD1)</td><td>l<sub>hyd_1</sub> (m)</td><td>6</td><td>6,050</td><td>4,430</td><td>1,350</td></tr><tr><td>k<sub>hyd_1</sub> (mm²/h)</td><td>5</td><td>4,140</td><td>3,640</td><td>2,990</td></tr><tr><td>E<sub>min_1</sub> (mm)</td><td>10</td><td>8.24</td><td>12.5</td><td>3.54</td></tr><tr><td>E<sub>max_1</sub> (mm)</td><td>2</td><td>124</td><td>142</td><td>45</td></tr><tr><td>α<sub>_1</sub> (-)</td><td>4</td><td>1.19</td><td>1.39</td><td>0.46</td></tr><tr><td rowspan="5">Black flysch (HYD2)</td><td>Log(k<sub>is_1</sub>) (1/h)</td><td>9</td><td>-5.34</td><td>-4.74</td><td>-3.65</td></tr><tr><td>l<sub>hyd_2</sub> (m)</td><td>13</td><td>3,300</td><td>2,610</td><td>900</td></tr><tr><td>k<sub>hyd_2</sub> (mm²/h)</td><td>11</td><td>274</td><td>405</td><td>287</td></tr><tr><td>E<sub>min_2</sub> (mm)</td><td>15</td><td>2.58</td><td>2.94</td><td>1.19</td></tr><tr><td>E<sub>max_2</sub> (mm)</td><td>12</td><td>27.8</td><td>18.61</td><td>6.46</td></tr><tr><td rowspan="5">Lower compartment</td><td>α<sub>_2</sub> (-)</td><td>14</td><td>1.8</td><td>1.03</td><td>0.58</td></tr><tr><td>Log(k<sub>MC</sub>) (mm/h)</td><td>16</td><td>-4.03</td><td>-3.14</td><td>-1.69</td></tr><tr><td>a<sub>MC</sub> (-)</td><td>17</td><td>2.33</td><td>2.01</td><td>0.59</td></tr><tr><td>C<sub>loss</sub> (mm)</td><td>3</td><td>0.37</td><td>0.69</td><td>0.6</td></tr><tr><td>k<sub>CS</sub> (mm/h)</td><td>1</td><td>0.56</td><td>1.42</td><td>2.46</td></tr><tr><td></td><td>a<sub>CS</sub> (-)</td><td>7</td><td>1.48</td><td>2.39</td><td>0.94</td></tr></tbody></table>

corresponding sensitivity ranking. Among the model parameters, the discharge coefficient for the flow from the conduit to the spring $k_{CS}$ is the most sensitive, followed in order by the upper storage threshold for hydrotepe 1 $E_{max_1}$ and by the storage threshold for the conduit $C_{loss}$. On the contrary, the discharge coefficient $k_{MC}$ and the exponent $a_{MC}$, which are both related to the transfer between the matrix and conduit ($Q_{MC}$), are the least sensitive. In between the most and least sensitive parameters, the parameters related to the flow components from the hydrotepe 1, i.e., limestone, are more sensitive than those from the hydrotepe 2, i.e., black flysch.

Figure 3b shows the scatter plot of the computed $\sigma$ and $\mu^*$, delineating the areas of the graph corresponding to

monotonic and/or linear interactions and nonlinear non-monotonic interactions (Sanchez et al. 2014). Most of the parameters fall above the bisector, indicating strong non-linear nonmonotonic interaction among most of the model parameters.

The sensitivity of the parameters was investigated by considering the distribution of the performance of the Morris model realizations over the parameter ranges (Table 2; Fig. 4). Figure 4 shows the distribution of occurrence of the metric KGE of the model realizations for the most and least sensitive parameters, which are the discharge coefficient of the flow from the conduit to the spring $k_{CS}$ and the exponent of the transfer between the

<figure><img src="image_9.png" /><figcaption>**Fig. 4** Distribution of occurrence of the Kling-Gupta Efficiency (KGE) of the model realizations obtained by means of the Morris screening sensitivity analysis for a the most sensitive parameter, i.e., the recession coefficient of the flow from the conduit to the spring ($k_{CS}$), and b for the least sensitive parameter, i.e., the exponent of the flow transfer between the matrix and the conduit ($a_{MC}$)</figcaption></figure>
matrix and conduit $a_{MC}$, respectively (Fig. 3a; Table 3). The distribution of the realizations with KGE lower than zero, larger than zero and larger than 0.5 confirms the sensitivity rank obtained with the Morris analysis (Table 3). Indeed, for the exponent $a_{MC}$, all the three distributions of KGE are uniform over the parameter range (Fig. 4b). On the contrary, for the discharge coefficient $k_{CS}$, the realizations with KGE larger than 0 and larger than 0.5 show a probability of occurrence skewed towards the right side of the parameter range (Fig. 4a).

## Model calibration and validation

Among the 9,000 Morris realizations, the parameter sample corresponding to the highest KGE was chosen as the selected parameter set and the corresponding realization was chosen as the selected discharge. The selected parameters,

calibrated for the period 1/03/2022–30/03/2022, are shown in Table 3 and lead to a KGE of 0.89. The model was then validated with the same selected parameter values for the period 30/03/2022–30/04/2022 leading to a KGE of 0.8. Figure 5a shows the calibrated and validated discharge time series together with the observed discharge at the spring. Overall, the simulated discharge matches the rising and falling limbs of the observed hydrograph, with the exception of the recession phase between April 4 and April 18, which is underestimated by 0.08 m³/s on average. For a better visualization, Fig. S3a in the ESM shows the time series of the observed discharge, the selected discharge of the Morris realizations based on KGE and the difference between the observed and selected discharges at the Las Hountas spring.

To account for the uncertainty in the model results, the distribution of the Morris behavioral simulations with a KGE greater than 0.5 was considered for both the calibration

<figure><img src="image_10.png" /><figcaption>**Fig. 5** Simulated (red lines) flow components ($Q_{spring}$, $Q_{hyd_1}$, $Q_{hyd_2}$, $Q_{is_1}$, $Q_{MC}$ and $Q_{Closs}$) of the selected simulation together with their uncertainty bands computed from the distribution of the behavioral simulations with a KGE larger than 0.5 (gray bands). The dark grey and grey areas of the bands represent the interquartile and the 10–90% percentile ranges of the behavioral simulations, respectively, while the black line is the median of the distribution. The blue line is the discharge observed at the spring</figcaption></figure>
<caption>**Table 4** Relative width of the uncertainty bands computed from the distribution of the behavioral simulations (KGE > 0.5) and of the constrained simulations (hydrochemical constraints) for the flow components $Q_{spring}$, $Q_{hyd_1}$, $Q_{hyd_2}$, $Q_{is_1}$, $Q_{MC}$ and $Q_{Closs}$. Both the</caption>

relative width of the interquartile and 10–90% percentile ranges were normalized considering the mean of the behavioral and constrained simulations

<table><thead><tr><th rowspan="2">Flow component</th><th colspan="2">Behavioral simulations</th><th colspan="2">Constrained simulations</th></tr><tr><th>Relative width interquartile range</th><th>Relative width 10–90% percentile range</th><th>Relative width interquartile range</th><th>Relative width 10–90% percentile range</th></tr></thead><tbody><tr><td>$Q_{spring}$</td><td>0.345</td><td>0.613</td><td>0.185</td><td>0.555</td></tr><tr><td>$Q_{hyd_1}$</td><td>0.386</td><td>0.774</td><td>0.336</td><td>0.941</td></tr><tr><td>$Q_{hyd_2}$</td><td>0.822</td><td>1.44</td><td>0.396</td><td>0.844</td></tr><tr><td>$Q_{is_1}$</td><td>2.66</td><td>1.99</td><td>0.942</td><td>1.57</td></tr><tr><td>$Q_{MC}$</td><td>2.66</td><td>2.00</td><td>0.967</td><td>1.64</td></tr><tr><td>$Q_{Closs}$</td><td>0.39</td><td>2.41</td><td>0.962</td><td>3.1</td></tr></tbody></table>

and the validation periods. Figure 5 shows, for each flow component, i.e., $Q_{spring}$, $Q_{hyd_1}$, $Q_{hyd_2}$, $Q_{is_1}$, $Q_{MC}$, and $Q_{Closs}$, the interquartile and 10–90% percentile envelopes as well as the median of the 448 selected simulations. For each flow component, the selected simulated discharge falls inside the interquartile range. The simulated discharge at the spring $Q_{spring}$ and the flow from the hydrotepe 1 to the conduit $Q_{hyd_1}$ show the lowest relative width of the interquartile uncertainty band (Fig. 5a,b; Table 4). On the contrary, the flow components characterized by larger uncertainty are the infiltration from hydrotepe 1 to the matrix $Q_{is_1}$ and the transfer between the matrix and conduit $Q_{MC}$ (Fig. 5d,e; Table 4). The relative width of the 10–90% percentile ranges follows the same behavior as the interquartile ranges: $Q_{spring}$ and $Q_{hyd_1}$ show the narrowest bands (Fig. 5a,b; Table 4); whereas $Q_{is_1}$, $Q_{MC}$ and $Q_{Closs}$ are characterized by the

largest uncertainty (Fig. 5d–f; Table 4). The drainage from the conduit $Q_{Closs}$ is active only during high peak discharges, whereas it is null at low flow conditions (Fig. 5f).

## Model selection using hydrochemical constraints

After validating the model by only considering the KGE metric, it was investigated as to whether the internal dynamics of the system, i.e., water contributions from the limestone and black flysch, were represented or not. The observed temporal dynamics of the system show high temporal variability, justifying the high-temporal resolution of hydrochemical data collection. Here, daily sampling frequencies or averaged daily values would be unsuitable to represent the variability of the water contributions. For this purpose, the high-temporal resolution hydrochemistry

<figure><img src="image_11.png" /><figcaption>**Fig. 6** Water level (wl) (m), electrical conductivity EC (µcm/S) and weight factors $f_{HCO_3^-}$ (-) and $f_{SO_4^{2-}}$ (-) time series observed at the Las Hountas spring during the increase in water level—occurred in the 2 days 30/03/2022 and 1/04/2022—for the event 30/03/2022–</figcaption></figure>

7/04/2022. The weight factors are computed as the contribution of the individual free ion to the total EC (Richieri et al. 2023). The three green panels indicate the times at which there is a simultaneous increase in $f_{SO_4^{2-}}$ (-) and decrease in $f_{HCO_3^-}$ (-)
data collected in April 2022 (Table 1) were used together with the understanding of the system functioning. Richieri et al. (2023) found out that it is possible to distinguish the water contributions from the two main geological formations present in the Baget catchment by looking at the contribution to the total EC, called weight factor $f_i$, of $HCO_3^-$ and $SO_4^{2-}$. More precisely, a simultaneous increase in $f_{SO_4^{2-}}$ and decrease in $f_{HCO_3^-}$ indicated an increase in the relative contribution from the black flysch (Richieri et al. 2023). Figure 6 shows water level (wl) (m), EC (µcm/S) and weight factors $f_{HCO_3^-}$(-) and $f_{SO_4^{2-}}$ (-) time series observed at the spring during the increase in water level, which occurred between 30/03/2022 and 1/04/2022, for the event 30/03/2022–7/04/2022. Three peaks in $f_{SO_4^{2-}}$, simultaneous to a low point in $f_{HCO_3^-}$ were observed on May 30 at 8 pm, May 31 at 9.30 am and May 31 at 10 pm. However, the simulation selected among the 9,000 Morris realizations (highest KGE) does not represent the observed varying relative water contributions. Indeed, the flow from the limestone to the conduit $Q_{hyd_1}$ and the flow from the black flysch to the conduit $Q_{hyd_2}$ are almost parallel lines with a null transfer between the matrix and conduit $Q_{MC}$ (Fig. 7a).

Based on the observed hydrochemical data (Fig. 6), three hydrochemical constraints were used to identify among the Morris realizations with KGE larger than 0.5 those also respecting the internal dynamics of the system. The first criterion defines that the standard deviation of the transfer

between the matrix and conduit $Q_{MC}$ needs to be significant. Here, the presented conceptual model requires that an exchange between the matrix and conduit exists. To define a threshold value to consider the significance of the exchange, it is considered a standard deviation computed over the calibration period, which has to be larger than the mean standard deviation of the behavioral Morris simulations (0.004). In fact, many behavioral simulations lead to zero or minimal exchange between the two compartments. Other thresholds for the minimum accepted standard deviation were also tested: smaller thresholds (e.g., 0.002) had limited impact on the selected simulations respecting all three criteria, whereas larger thresholds (e.g., 0.01) led to no simulation respecting all three criteria. The second criterion states that the increase in the flow from the black flysch to the conduit ($Q_{hyd_2}$) during the rising limb of the hydrograph must be larger than the increase in the flow from the limestone to the conduit ($Q_{hyd_1}$) in the same period. Finally, the third criterion defines that, on April 1, the flow from the black flysch to the conduit ($Q_{hyd_2}$) must be lower than the flow from the limestone to the conduit ($Q_{hyd_1}$). Figure 7b shows the median and interquartile ranges of $Q_{hyd_1}$, $Q_{hyd_2}$ and $Q_{MC}$ for the seven Morris realizations respecting the three defined hydrochemical constraints. The median of the flow from the limestone ($Q_{hyd_1}$) and black flysch ($Q_{hyd_2}$) to the conduit captures the relative increases in water contribution from the black flysch observed on May 30 at 8 pm, May 31 at 9.30 am and May 31 at 10 pm, with the first peak anticipated at ~12 h (Figs. 7b and 6). The median of the transfer between

<figure><img src="image_12.png" /><figcaption>**Fig. 7** Simulated discharge time series for the transfer between the matrix and conduit $Q_{MC}$, the flow from hydrotepe 1 (limestone) to conduit $Q_{hyd_1}$ and the flow from hydrotepe 2 (black flysch) to the conduit $Q_{hyd_2}$ during the increase in water level—occurred in the 2 days 30/03/2022 and 1/04/2022—for the event 30/03/2022–7/04/2022. **a** Discharge time series from the Morris simulation</figcaption></figure>

selected based on the KGE of the spring discharge. **b** Median and interquartile range of the behavioral Morris simulations selected with the three hydrochemical constraints. The three green arrows indicate the times at which the increase in the contribution from the black flysch $Q_{hyd_2}$ is larger than that from the limestone $Q_{hyd_1}$
<figure><img src="image_13.png" /><figcaption>**Fig. 8** Uncertainty bands computed from the distribution of the behavioral simulations with a KGE larger than 0.5 and respecting the three hydrochemical constraints for the flow components $Q_{spring}$, $Q_{hyd_1}$, $Q_{hyd_2}$, $Q_{is_1}$, $Q_{MC}$ and $Q_{Closs}$ (gray bands). The dark grey and grey areas of the bands represent the interquartile and the 10–90% percentile ranges of the behavioral simulations, respectively, while the black line is the median of the distribution. The blue line is the discharge observed at the spring</figcaption></figure>

the matrix and conduit (Fig. 7b) is positive during the entire rising limb of the hydrograph, indicating flow direction from the matrix to the conduit.

Figure 8 shows the distribution of the constrained Morris simulations, i.e., the simulations respecting to both the constraint on the spring discharge (acceptable simulations with KGE larger than 0.5) and the three hydrochemical constraints. The relative interquartile ranges of the constrained simulations (Table 4) are reduced in comparison to the distribution of the simulations considered behavioral based on the discharge KGE (Fig. 5) for all the flow components except $Q_{Closs}$. As observed in Fig. 5, $Q_{spring}$ (Fig. 8a) and $Q_{hyd_1}$ (Fig. 8b) are characterized by the lowest relative width of the interquartile band, whereas $Q_{is_1}$ (Fig. 8d), $Q_{MC}$ (Fig. 8e) and $Q_{Closs}$ (Fig. 8f) by the largest (Table 4). The relative width of the 10–90% percentile ranges is less significantly reduced than the interquartile ranges (Figs. 5 and 8). $Q_{spring}$ (Fig. 8a) shows the lowest relative width of the percentile band, while $Q_{MC}$ (Fig. 8e) and $Q_{Closs}$ (Fig. 8f) show the largest (Table 4).

Despite the reduction in the uncertainty bands, the median of the distribution of $Q_{spring}$ does not respect the plateau of 0.6 m³/s (Fig. S1 of the ESM). $Q_{hyd_1}$ and $Q_{Closs}$ are the only flow components displaying a larger relative width of the percentile band (Fig. 8b,f) in comparison to that observed without considering the hydrochemical constraints (Fig. 5; Table 4). Finally, the 10–90% percentile bands for $Q_{hyd_2}$ and $Q_{Closs}$ are characterized by a relative width that increases proportionally with the discharge.

Even if the 10–90% percentile range for $Q_{spring}$ is overall decreased (Fig. 8a), the percentile band is reduced for the calibration and increased for the validation period (Figs. 5a and 8a). On the contrary of what was noted before applying the hydrochemical constraints (Fig. 5a), the observed discharge at the spring falls outside the 10–90% percentile band for the calibration and inside the percentile band for the validation period (Fig. 8a). This results from the fact that the 10–90% percentile range of $Q_{spring}$, respecting the hydrochemical constraints (Fig. 8a), includes realizations
which were excluded from the 10–90% percentile range of the realizations selected based only on KGE (Fig. 5a).

To compare the performance of the model in simulating the discharge at the spring before and after considering the hydrochemical constraints, Fig. S3b (ESM) shows the time series of the difference between observed discharge and constrained discharge. The constrained discharge corresponds to the realization with the highest KGE among the subset of realizations respecting the hydrochemical constraints.

## Discussion

This section reports the novelties of the new LuKARS 2.0 conceptual model and the use of hydrochemical criteria to verify the model concept and parametrization.

### New LuKARS conceptual model

LuKARS 2.0 was developed at an hourly scale starting from the original LuKARS model at a daily scale (Bittner et al. 2018). When comparing the values of the parameters which are in common to both models, the discharge coefficient of the fast flow $k_{hy}$ and of the infiltration to the matrix $k_{is}$ are the parameters showing the largest variation between daily (reported in Bittner et al. 2018) and hourly simulations (this study). This difference in magnitude can be explained by the fact that the karst system in which Bittner et al. (2018) tested LuKARS, i.e., Kerschbaum, does not show relevant concentrated recharge processes and thus is characterized by a more pronounced infiltration than the Baget system, which shows a natural subsurface drainage rapidly recharged by sinkholes along the Lachein stream (Mangin 1975; Sivelle and Labat 2019). Thus, the response of Kerschbaum can be described by means of lower $k_{hy}$ and larger $k_{is}$ than the Baget system.

To represent the plateau in the discharge time series at the Las Hountas spring (Fig. S1 of the ESM), a drain is implemented from the conduit out of the catchment. The sensitivity analysis indicates that the parameter controlling the activation of the drainage ($C_{loss}$) is the third most sensitive parameter (Table 3). The selected Morris simulation as well as the interquartile range of the behavioral simulations (Fig. 5) confirm the relevance of the drain to discharge out of the system the excess water (Fig. 5f) and thus respect the maximum threshold of 0.6 m³/s at the spring (Fig. S1 of the ESM). However, the median of the behavioral simulations coincides with the 25% quartile and is constantly null throughout the calibration and validation periods, meaning that $Q_{Closs}$ does not play a role for 50% of the realizations (Fig. 5f). The behavioral simulation respecting the

hydrochemical data shows the same behavior (Fig. 8f): the median of the $Q_{Closs}$ is constantly null, while the interquartile and 10–90% percentile bands show the large impact that the drainage from the conduit has for a share of the realizations. This indicates, as discussed subsequently, that further data at high temporal resolution, such as solutes and isotopes time series, should be considered to constrain the model parameters. In addition, time series collected at different points of the catchment would give more insight into the role of the conduits as natural drains as well as into the dynamics of the water draining from the different recharge areas.

To represent the characteristic dual response behavior of karst systems, the original LuKARS model (Bittner et al. 2018) was further modified by implementing two buckets in the lower compartment, i.e., matrix and conduit, and a transfer $Q_{MC}$ between these two. The sensitivity analysis indicates that the parameters controlling the transfer between the matrix and conduit $Q_{MC}$, i.e., the discharge coefficient $k_{MC}$ and exponent $a_{MC}$, are the least sensitive parameters (Table 3; Fig. 3a). The selected Morris simulation and the median of the realizations with KGE larger than 0.5 are characterized by no and little influence of $Q_{MC}$ on the simulated discharge at the spring, respectively (Fig. 5e). The role of the transfer between the matrix and conduit becomes more relevant after applying the hydrochemical constraints: the realizations respecting the hydrochemical data (Fig. 6) present an increase in the width of the interquartile range and a median larger than the mean standard deviation of the behavioral Morris simulations, i.e., 0.004 (Fig. 8e).

Overall, the selected Morris realization well represents the dynamics of the spring discharge $Q_{spring}$ for both the calibration and validation, with an exception for the period between April 4 and 18 in which the simulated discharge is lower than the observed discharge (Fig. 5a). The model cannot represent the increase in discharge at the spring on April 4 due to the input data, which do not show significant precipitation on April 4 that could cause the increase in simulated discharge (Fig. S1, ESM). The hourly precipitation data were recorded at the Antichan gauging station, ~8 km north-east from the Las Hountas spring, and used as spatially homogeneous input for the presented conceptual model. Due to the heterogeneity of the precipitation field, it might be that the data are not always accurate for the Baget catchment during this specific rainfall event. Consequently, the entire recession curve until April 18 systematically underestimates the observed discharge of 0.08 m³/s on average. The raster daily E-OBS precipitation data (Cornes et al. 2018) were used as a comparison with the Antichan gauging station and showed consistent observations. However, due to the daily temporal resolution and the coarse spatial resolution of the grid (0.25°) of E-OBS, that dataset could not be used for the present study.
## Sensitive parameters

The results of the sensitivity analysis can be qualitatively compared to those obtained by Bittner et al. (2020) by applying the active subspaces on the daily LuKARS. In the present paper, the parameters related to hydrotope 1, i.e., limestone, are more sensitive than those from hydروtopе 2, i.e., black flysch, probably since the limestone represents 70% of the entire catchment area and thus, assuming spatially homogeneous precipitation rate, contributes to a larger extent to the total recharge to the spring. This observation agrees with the results from Bittner et al. (2020), where the parameters related to the hydrometer with a larger area are generally more sensitive. On the contrary, the sensitivity of the other parameters differs between the present study and the previous investigation due to the different representations of the lower compartment in the model concept and the different study areas. Indeed, according to Bittner et al. (2020), the discharge coefficient regulating the infiltration to the matrix $k_{is}$ is the most sensitive parameter, whereas in the present study $k_{is}$ shows a sensitivity ranking of 9 out of 17 (Table 4).

The parameters regulating the transfer between the matrix and conduit, i.e., discharge coefficient $k_{MC}$ and exponent $a_{MC}$, are the least sensitive parameters of the presented model concept (Table 3; Fig. 3a). This is also due to the model structure. In fact, the matrix-conduit interaction can be merged as a single compartment, as, for example, in the original version of LuKARS. While this choice can be satisfactory for reproducing discharge time series, it can lead to erroneous hydrochemical interpretation of the results. Therefore, since the sensitivity is a function of the used objective function, the results of the sensitivity analysis performed with respect to KGE, and only against the discharge at the spring, may not capture the importance of all the parameters, especially when nonlinearities are present (Fig. 3b). Indeed, the selection of the Morris realizations respecting the hydrochemical constraints show the relevance of the transfer between the matrix and conduit. Those realizations characterized by a transfer between the matrix and conduit $Q_{MC}$ with a standard deviation larger than the mean standard deviation of the behavioral Morris simulations (Fig. 7b) well represent the observed hydrochemical dynamics (Fig. 6).

## Impact of hydrochemical information on model selection

The model was calibrated and validated considering the discharge and the hydrochemical signal at the spring, respectively, over a time period of 9 days. The short period of time was chosen in line with the performed analyses so as to have an event-based model. In addition, despite the limitations

related to calibration over a short period of time (Leins et al. 2023), the length of the calibration and validation periods were chosen based on the collected high-resolution hydrochemical data (Table 1). Indeed, the interpretation of the internal fluxes of a longer simulated period may require the consideration of reactive transport, i.e., chemical reaction with soil and aquifer materials, seasonal variation in background concentrations, and coexistence of water with different ages.

Based on the previous investigation into the system functioning of Richieri et al. (2023), the parameter ranges (Table 2) for hydrotopes 1 and 2, i.e., limestone and black flysch, respectively, were defined to enhance a faster reaction of the black flysch in case of a strong precipitation event. In fact, it was observed that when the water level at the spring rises simultaneously to the EC, the relative contribution to the EC of $SO_4^{2-}$, which originated from the pyrite contained in the black flysch (Eq. 2), increases. The faster response of the black flysch was considered by defining the lower upper and lower thresholds for the activation of the fast flow to the conduit (Table 2). However, when investigating the internal flow component of the selected Morris realization, the increase in fast flow from hydrometer 2, i.e., black flysch, observed during the precipitation event from 30/03/2022 to 7/04/2022 (Fig. 6), is not captured (Fig. 7a). Hydrochemical constraints were thus implemented to constrain those realizations respecting the dynamics of $f_{SO_4^{2-}}$ and $f_{HCO_3^-}$ (Fig. 6). The constrained simulations (Fig. 7b) capture the larger relative increases in the contribution from the black flysch and the lower relative increases in the contribution from the limestone that occurred on May 30 at 8 pm, May 31 at 9.30 am and May 31 at 10 pm. The first relatively faster increase in black flysch on May 30 at 8 pm (Fig. 7b) was anticipated with respect to the observed hydrochemical data (Fig. 6) of ~12 h, probably due to the initial water content in the soil before the rise in water level and the need of larger storage.

The verification of the model by considering hydrochemical constraints resulted in the overall reduction of the interquartile uncertainty bands with respect to the case in which the Morris realizations were evaluated based only on the KGE (Figs. 5 and 8). This is due to the application of the hydrochemical constraints, which reduce the selected simulations from 448 to 8. Therefore, the hydrochemical criteria may be too stringent for the validation of the parameter sets, but they can be used to verify whether the conceptual model applies to simulations and can, therefore, be used to test the plausibility of the conceptual model. Despite the reduction of the uncertainty bands, the median of the realizations respecting the hydrochemical constraints of $Q_{spring}$ underestimates the baseflow and exceeds the plateau of 0.6 m³/s (Fig. S1 of the ESM), which is respectively both the selected

<figure><img src="image_14.png" /></figure>
Morris realizations and the median of the Morris realizations with a KGE larger than 0.5. Finally, the relevance of the hydrochemical constraints is also seen by the fact that the percentile range for the $Q_{\text{spring}}$ of the realizations selected, based only on KGE (Fig. 5a), does not contain some realizations that are included in the percentiles of the realizations respecting the hydrochemical dynamics (Fig. 8a). Therefore, without considering the hydrochemical constraints, some parameter sets representing the hydrochemical dynamics would have been excluded.

An alternative to the approach demonstrated in this study, hydrochemical data could be used during the model calibration process. Here, the sensitive parameters of the simulations for the hydrochemical dynamics could be further investigated and constrained. Among the most sensitive parameters (Fig. 3; Table 3), patterns can be identified for $\alpha_{1-}$, $k_{\text{hyd}_1}$ (mm$^2$/h), $E_{\text{max}_1}$ (mm), $k_{\text{is}_1}$ (1/h) and $k_{\text{hyd}_2}$ (mm$^2$/h), whose values fall in ranges narrower than those initially defined for the Morris analysis (Table 2).

## Conclusion

This study presents the development of LuKARS 2.0 on an hourly scale for the Baget karst catchment. The model concept was modified in comparison to the original daily LuKARS from Bittner et al. (2018) to represent the plateau of 0.6 m$^3$/s, which characterizes the discharge at the Las Hountas spring. Thus, a drain from the conduit was implemented to discharge the excess water out of the system as soon as the water level in the conduit rose above a defined threshold. The interaction between the matrix and conduit is here updated in comparison to Bittner et al. (2018) to allow a more flexible conceptualization of karst systems. The sensitivity analysis performed using the Morris method with a total of 9,000 realizations was applied to investigate the parameter sensitivity and calibrate the model. The model was calibrated and validated for the periods 1/03/2022–29/03/2022 and 30/03/2022–30/04/2022, respectively, considering KGE as a performance metric. Among the Morris realizations, the behavioral simulations were first defined based on the KGE value of the simulated discharge at the spring. Then, three additional hydrochemical constraints were considered to represent the internal dynamics of the systems, i.e., water contributions from the limestone and black flysch present in the catchment. The hydrochemical constraints were defined based on $SO_4^{2-}$ and $HCO_3^-$ time series at high resolution and the previous work of Richieri et al. (2023).

The results of the present investigation show that the selected Morris realizations can represent the dynamics of the spring discharge but do not sufficiently account for the temporal variation of the contribution from the limestone

and black flysch, which was observed during the rising limb of the event 30/03/2022–7/04/2022. The implementation of hydrochemical constraints leads to the representation of the relative increase in black flysch, while still providing a good rendering of the spring discharge. The uncertainty quantification performed on those realizations respecting the hydrochemical data show a reduction in the relative width of the interquartile bands in comparison to before the implementation of the hydrochemical constraints. The selected simulation and the interquartile range of the behavioral simulations show the relevance of the implemented drainage from the conduit to represent a plateau in the spring discharge. However, the 10–90% percentile range of the spring discharge exceeds the defined threshold, indicating that further constraints in the parameter ranges would be required. In addition, further validation using other datasets would allow for investigating to what extent the information contained in hydrochemical data time series leads to a better representation of the system internal dynamics.

In conclusion, this study highlights the importance of considering the internal dynamics of a system when selecting hydrological models. Indeed, even if different model realizations simulated the spring discharge with a comparable KGE, not all of them were consistent with the dynamics of the measured $SO_4^{2-}$ and $HCO_3^-$ time series. Hence, hydrochemical data at high temporal resolution need to be collected to conceptualize and select hydrological models as well as to properly constrain the model parameter ranges.

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1007/s10040-024-02801-2.

**Acknowledgements** The authors thank the associate editor and the reviewers for their valuable comments and constructive feedback which significantly contributed to the improvement of this paper.

**Funding** Open Access funding enabled and organized by Projekt DEAL. This research is part of the ROCKAT project (RObust Conceptualisation of KArst Transport), funded by Deutsche Forschungs-gemeinschaft (DFG), CH 981/6-1 and HA 8113/6–1.

**Data availability** The hydrochemical data at the Las Hountas spring were collected during the fieldwork in March and April 2022. The dataset can be downloaded from the Hydroshare repository at the following link http://www.hydroshare.org/resource/d0d733f9505f4b289e84d99729fb16dd

## Declarations

**Conflicts of interests** On behalf of all authors, the corresponding author states that there is no conflict of interest.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated
otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## References

AL Khoury I, Boithias L, Bailey RT, Ollivier C, Sivelle V, Labat D (2023) Impact of land-use change on karst spring response by integration of surface processes in karst hydrology: the ISPEEKH model. J Hydrol 626(12):130300. https://doi.org/10.1016/j.jhydrol.2023.130300

Andréassian V (2023) On the (im)possible validation of hydrogeological models. C R Géosci 355:1–9. https://doi.org/10.5802/crgeos.142

Barbieri M, Boschetti T, Petitta M, Tallini M (2005) Stable isotope (2H, 18O and 87Sr/86Sr) and hydrochemistry monitoring for groundwater hydrodynamics analysis in a karst aquifer (Gran Sasso, Central Italy). Appl Geochem 20(11):2063–2081. https://doi.org/10.1016/j.apgeochem.2005.07.008

Bennett ND, Croke BFW, Guariso G, Guillaume JHA, Hamilton SH, Jakeman AJ, Marsili-Libelli S, Newham LTH, Norton JP, Perrin C, Pierce SA, Robson B, Seppelt R, Voinov AA, Fath BD, Andreassian V (2013) Characterising performance of environmental models. Environ Model Softw 40:1–20. https://doi.org/10.1016/j.envsoft.2012.09.011

Berthelin R, Hartmann A (2020) The shallow subsurface of karst systems: review and directions. In: Bertrand C, Denimal S, Steinmann M, Renard P (eds) Eurokarst 2018, Besançon: advances in karst science, April Issue. Springer, pp 61–68. https://doi.org/10.1007/978-3-030-14015-1_7

Beven K, Freer J (2001) Equifinality, data assimilation, and uncertainty estimation in mechanistic modelling of complex environmental systems using the GLUE methodology. J Hydrol 249(1–4):11–29. https://doi.org/10.1016/S0022-1694(01)00421-8

Bittner D, Narany TS, Kohl B, Disse M, Chiogna G (2018) Modeling the hydrological impact of land use change in a dolomite-dominated karst system. J Hydrol 567:267–279. https://doi.org/10.1016/j.jhydrol.2018.10.017

Bittner D, Teixeira Parente M, Mattias S, Wohlmuth B, Chiogna G (2020) Identifying relevant hydrological and catchment properties in active subspaces: an inference study of a lumped karst aquifer model. Adv Water Resour 135:103472. https://doi.org/10.1016/j.adwatres.2019.103472

Boyd CE (2020) Carbon dioxide, pH, and alkalinity. In: Water quality. Springer, Cham, Switzerland. https://doi.org/10.1007/978-3-030-23335-8_9

Campolongo F, Cariboni J, Saltelli A (2007) An effective screening design for sensitivity analysis of large models. Environ Model Softw 22:1509–1518. https://doi.org/10.1016/j.envsoft.2006.10.004

Chang Y, Hartmann A, Liu L, Jiang G, Wu J (2021) Identifying more realistic model structures by electrical conductivity observations of the karst spring. Water Resour Res 57:e2020WR028587. https://doi.org/10.1029/2020WR028587

Chiogna G, Marcolini G, Engel M, Wohlmuth B (2024) Sensitivity analysis in the wavelet domain: a comparison study. Stoch Env Res Risk Assess. https://doi.org/10.1007/s00477-023-02654-3

Cinkus G, Mazzilli N, Jourde H, Wunsch A, Liesch T, Ravbar N, Chen Z, Goldscheider N (2022) When best is the enemy of good: critical evaluation of performance criteria in hydrological models. Hydrol Earth Syst Sci Discuss. https://doi.org/10.5194/hess-2022-380

Cornes R, van der Schrier G, van den Besselaar EJM, Jones PD (2018) An ensemble version of the E-OBS temperature and precipitation datasets. J Geophys Res: Atmos 123(17). https://doi.org/10.1029/2017JD028200

Debroas É-J (2009) Géologie du bassin versant du Baget (zone nord-pyrénéenne, Ariège, France): nouvelles observations et consequences [Geology of the Baget watershed (north-Pyrenean zone, Ariège, France): new observations and consequences]. Strata 3(46):93

de Ferreira PML, da Paz AR, Bravo JM (2020) Objective functions used as performance metrics for hydrological models: state-of-the-art and critical analysis. RBRH 25:e42. https://doi.org/10.1590/2318-0331.252020190155

Ford DC, Williams PW (2007) Karst hydrogeology and geomorphology. Wiley, Chichester, UK, 567 pp

Gil-Márquez JM, Barberá JA, Andreo B, Mudarra M (2017) Geochemical evolution of groundwater in an evaporite karst system: Brujuelo Area (Jaén, S Spain). Proc Earth Planet Sci 17:336–339. https://doi.org/10.1016/j.proeps.2016.12.085

Goldscheider N, Chen Z, Auler AS, Bakalowicz M, Broda S, Drew D, Hartmann J, Jiang G, Moosdorf N, Stevanovic Z, Veni G (2020) Global distribution of carbonate rocks and karst water resources. Hydrogeol J 28:1661–1677. https://doi.org/10.1007/s10040-020-02139-5

Gupta HV, Wagener T, Liu Y (2008) Reconciling theory with observations: elements of a diagnostic approach to model evaluation. Hydrol Process 22(18):3802–3813. https://doi.org/10.1002/hyp.6989

Gupta HV, Kling H, Yilmaz KK, Martinez GF (2009) Decomposition of the mean squared error and NSE performance criteria: implications for improving hydrological modelling. J Hydrol 377:80–91. https://doi.org/10.1016/j.jhydrol.2009.08.003

Hagedorn B, Whittier RB (2015) Solute sources and water mixing in a flashy mountainous stream (Pahsimeroi River, U.S. Rocky Mountains): implications on chemical weathering rate and groundwater–surface water interaction. Chem Geol 391:123–137. https://doi.org/10.1016/j.chemgeo.2014.10.031

Hartmann A, Antonio Barberá J, Andreo B (2017) On the value of water quality data and informative flow states in karst modelling. Hydrol Earth Syst Sci 21(12). https://doi.org/10.5194/hess-21-5971-2017

Hartmann A, Wagener T, Rimmer A, Lange J, Brielmann H, Weiler M (2013) Testing the realism of model structures to identify karst system processes using water quality and quantity signatures. Water Resour Res 49(6). https://doi.org/10.1002/wrcr.20229

Hartmann A, Goldscheider N, Wagener T, Lange J, Weiler M (2014) Karst water resources in a changing world: review of hydrological modeling approaches. Rev Geophys 52:218–242. https://doi.org/10.1002/2013RG000443

Jaxa-Rozen M, Kwakkel J (2018) Tree-based ensemble methods for sensitivity analysis of environmental models: a performance comparison with Sobol and Morris techniques. Environ Model Softw 107:245–266. https://doi.org/10.1016/j.envsoft.2018.06.011

Klemeš V (1986) Dilettantism in hydrology: transition or destiny? Water Resour Res 22(9S):177S-188S. https://doi.org/10.1029/WR022i09Sp0177S

Knoben WJM, Freer JE, Woods RA (2019) Technical note: inherent benchmark or not? comparing Nash-Sutcliffe and Kling-Gupta efficiency scores. Hydrol Earth Syst Sci 23(10). https://doi.org/10.5194/hess-23-4323-2019

Labat D, Ababou R, Mangin A (1999) Linear and nonlinear input/output models for karstic springflow and flood prediction at different time scales. Stoch Env Res Risk Assess 13(5):337–364. https://doi.org/10.1007/s004770050055
Le Moine N, Andréassian V, Mathevet T (2008) Confronting surface- and groundwater balances on the La Rochefoucault-Touvre karstic system (Charente, France). Water Resour Res 44(3):1–10. https://doi.org/10.1029/2007WR005984

Leins T, Liso IS, Parise M, Hartmann A (2023) Evaluation of the predictions skills and uncertainty of a karst model using short calibration data sets at an Apulian cave (Italy). Environ Earth Sci 82(14):351. https://doi.org/10.1007/s12665-023-10984-2

Mangin A (1975) Contribution à l'étude hydrodynamique des aquifères karstiques [Contribution to the hydrodynamic study of karst aquifers]. Annal Spéléol 29(3):283-332; 29(4):495-601; 30(1):21-124

Mazzilli N, Guinot V, Jourde H, Lecoq N, Labat D, Arfib B, Baudement C, Danquigny C, Dal Soglio L, Bertin D (2019) KarstMod: a modelling platform for rainfall-discharge analysis and modelling dedicated to karst systems. Environ Model Softw 122:103927. https://doi.org/10.1016/j.envsoft.2017.03.015

Mazzilli N, Sivelle V, Cinkus G, Jourde H, Bertin D (2023) Karst-Mod user guide, version 3.0. hal-01832693v2, French SNO Karst. http://www.sokarst.org. Accessed May 2024

Merchán-Rivera P, Geist A, Disse M, Huang J, Chiogna G (2022) A Bayesian framework to assess and create risk maps of ground-water flooding. J Hydrol 610:127797. https://doi.org/10.1016/j.jhydrol.2022.127797

Moriasi DN, Arnold JG, Liew M, Bingner RL, Harmel RD, Veith TL (2007) Model evaluation guidelines for systematic quantification of accuracy in watershed simulations. Trans ASABE 50:885–900. https://doi.org/10.13031/2013.23153

Morris MD (1991) Factorial sampling plans for preliminary computational experiments. Technometrics 33:161–174. https://doi.org/10.1080/00401706.1991.10484804

Mudarra M, Hartmann A, Andreo B (2019) Combining experimental methods and modeling to quantify the complex recharge behavior of karst aquifers. Water Resour Res 1–21. https://doi.org/10.1029/2017WR021819

Padilla A, Pulido-Bosch A, Mangin A (1994) Relative importance of baseflow and quickflow from hydrographs of karst spring. Groundwater 32:267–277. https://doi.org/10.1111/j.1745-6584.1994.tb00641.x

Parkhurst DL, Appelo CAJ (2013) Description of input and examples for PHREEQC version 3: a computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations. US Geological Survey Techniques and Methods, Book 6, Chapter A43. US Geological Survey, 497 pp. http://pubs.usgs.gov/tm/06/a43. Accessed May 2024

Richieri B, Bittner D, Hartmann A, Benettin P, van Breukelen BM, Labat D, Chiogna G (2023) Using continuous electrical conductivity measurements to derive major solute concentrations in karst systems. Hydrol Process 37:e14929. https://doi.org/10.1002/hyp.14929

Sanchez DG, Lacarrière B, Musy M, Bourges B (2014) Application of sensitivity analysis in building energy simulations: combining first- and second-order elementary effects methods. Energy Build 68:741–750. https://doi.org/10.1016/j.enbuild.2012.08.048

Sivelle V, Labat D (2019) Short-term variations in tracer-test responses in a highly karstified watershed. Hydrogeol J 27:2061–2075. https://doi.org/10.1007/s10040-019-01968-3

Sivelle V, Renard P, Labat D (2020) Coupling SKS and SWMM to solve the inverse problem based on artificial tracer tests in karstic aquifers. Water 12 (4). https://doi.org/10.3390/w12041139

Sivelle V, Jourde H, Bittner D, Richieri B, Labat D, Hartmann A, Chiogna G (2022a) Considering land cover and land use (LCLU) in lumped parameter modeling in forest dominated karst catchments. J Hydrol 612:128264. https://doi.org/10.1016/j.jhydrol.2022.128264

Sivelle V, Pérotin L, Ladouche B, de Montety V, Bailly-Comte V, Champollion C, Jourde H (2022b) A lumped parameter model to evaluate the relevance of excess air as a tracer of exchanged flows between transmissive and capacitive compartments of karst systems. Front Water 4:930115. https://doi.org/10.3389/frwa.2022.930115

Sivelle V, Cinkus G, Mazzilli N, Labat D, Arfib B, Massei N, Cousquer Y, Bertin D, Jourde H (2023) Improvement of the KarstMod modeling platform for a better assessment of karst groundwater resources. Hydrol Earth Syst Sci Discuss. https://doi.org/10.5194/hess-2023-17

Smith RC (2013) Uncertainty quantification: theory, implementation, and applications. SIAM, Computational Science and Engineering, Philadelphia, PA

Spence J, Telmer K (2005) The role of sulfur in chemical weathering and atmospheric CO2 fluxes: evidence from major ions, δ13CDIC, and δ34SSO4 in rivers of the Canadian Cordillera. Geochim Cosmochim Acta 69(23):5441–5458. https://doi.org/10.1016/j.gca.2005.07.011

Stevanović Z (2019) Karst waters in potable water supply: a global scale overview. Environ Earth Sci 78(23):1–12

Ulloa-Cedamanos F, Probst A, Dos-Santos V, Camboulive T, Granouil-lac F, Probst J-L (2021) Stream hydrochemical response to flood events in a multi-lithological karstic catchment from the Pyrenees Mountains (SW France). Water 13(13):1–25. https://doi.org/10.3390/w13131818

Ulloa-Cedamanos F, Probst J-L, Binet S, Camboulive T, Payre-Suc V, Pautot C, Bakalowicz M, Beranger S, Probst A (2020) A forty-year karstic critical zone survey (Baget catchment, Pyrenees-France): lithologic and hydroclimatic controls on seasonal and inter-annual variations of stream water chemical composition, pCO2, and carbonate equilibrium. Water 12(5):1227. https://doi.org/10.3390/w12051227

Wang Z, Guo X, Kuang Y, Chen Q, Luo M, Zhou H (2021) Recharge sources and hydrogeochemical evolution of groundwater in a heterogeneous karst water system in Hubei Province, Central China. Appl Geochem 136:105165. https://doi.org/10.1016/j.apgeochem.2021.105165

Weiss J (2020) Basics of ion chromatography: separation mechanisms, stationary phases, detection methods, and applications. Thermo Fisher, Dreieich, Germany

Winston WE, Criss RE (2004) Dynamic hydrologic and geochemical response in a perennial karst spring. Water Resour Res 40(5). https://doi.org/10.1029/2004WR003054

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.