Impact of seasonal snow on the recharge of a mountain karst aquifer under climate change: the Dévoluy case study (Southern Alps, France)

Seasonal snow strongly influences groundwater recharge in mountain aquifers, yet its role in mid-altitude karst systems under climate warming remains poorly quantified. We investigated the Dévoluy karst aquifer (Southern French Alps) to assess how snow controls recharge and how spring discharge may respond to rising temperatures. Using the KarstMod platform, we developed a rainfall-snow-discharge model incorporating a degree-day snow routine to partition precipitation between rainfall and snow, and simulate the snowmelt. The model was calibrated and validated over four contrasting years (two low-snow, one high-snow, and one very high-snow year). Results show that accounting for snow processes is essential to reproduce the observed discharge dynamics, highlighting the dominant role of snow accumulation and melting in controlling both flow timing and magnitude in karst environments. We then tested the sensitivity of the karst spring discharge to temperature perturbations. Under +2 and +4 °C warming scenarios, simulated winter flows increase while snowmelt peaks occur earlier, resulting in earlier and more severe summer low-flow periods. August discharge decreases by 26 % to 42 %, respectively, compared to present conditions. These findings demonstrate the critical role of seasonal snow in regulating recharge in mid-altitude karst aquifers and highlight that ongoing warming will substantially reduce summer water availability in mountain regions.

1 Introduction

Mid-altitude mountains play a major role in the water supply providing a significant proportion of the water to downstream plains (Dettinger et al., 2004; Fayad et al., 2017; Jódar et al., 2020; Taucare et al., 2020; Viviroli et al., 2007). In Europe, intermediate altitudes (1500 to 2500 m) have usually a strong snow component in their precipitation and seasonal snowpack constitutes a temporary water reserve. The melting of snow in the spring increases the flow of rivers, which supports the water availability throughout the catchment during dry periods (Barnett et al., 2005). Snowpack and snowmelt evolution also act as a regulator of recharge, delaying it over time. But this latter effect is difficult to quantify, as snowpack evolution depends on numerous factors such as deposition thickness, weather conditions, slope exposure and soil infiltration capacity (Bayard et al., 2005; Eeckman et al., 2025; Meeks and Hunkeler, 2015; Reisch and Toran, 2014).

The evolution of snow cover is essentially dependent on precipitation and temperature, making these systems particularly sensitive to global warming (Eckert et al., 2024). A reduction in snow cover or earlier melting can significantly alter hydrological regimes, particularly during summer low-flow periods or droughts when water demand is the highest (Zierl and Bugmann, 2005). As in other mid-mountain basins, this impact is particularly significant in the Alps, where water resources are extensively used for drinking water supply, irrigation and hydroelectricity (Andrew and Sauquet, 2017; Branche, 2016; Magand et al., 2015).

The role of snowpack and snowmelt evolution in the availability of surficial water resources was successfully addressed by hydrological modelling (Garçon, 1996; Valéry et al., 2014). However it remains poorly studied for hydrogeological systems and especially for karsts while carbonate rocks are widely distributed in the Alps (Chen et al., 2017; Goldscheider et al., 2020). These aquifers present a high degree of structural heterogeneity (De Waele and Gutierrez, 2022; Ford and Williams, 2007). They are formed by a network of interconnected conduits enclosed in a more or less porous and fissured matrix, allowing both fast and slow flows. This duality of flow velocities complexifies hydrological dynamics, in particular the restitution of discharge at the outlets.

The impact of snow on the hydrodynamics of karst systems is of growing interest. For example, using natural and artificial tracers, Jódar et al. (2020) have shown that snow cover facilitates diffuse recharge through slow infiltration, and that karst organization and flow are strongly linked to this aspect. Using modelling, Çallı et al. (2022) have shown that the proper representation of snow processes at basin scale significantly improves karst flow simulation. Some studies have coupled rainfall-discharge models with climate projections to investigate the effects of climate change on karst hydrological regimes (Doummar et al., 2018), enabling projections of discharge evolution under different warming scenario. Fan et al. (2023) examined a mountainous karst system and highlighted how the timing of snowmelt affects maximum spring discharge under climate change scenarios. These studies demonstrate that snowpack and snowmelt are key drivers of karst recharge and discharge, and that their proper representation is essential for accurately simulating hydrological responses in snow-influenced karst system.

Hydrological modelling is an efficient tool for studying the recharge of karst aquifers by snowmelt. Among the multiple available approaches, reservoir-type rainfall-snow-discharge models are commonly used in surface hydrology and hydrogeology. Models such as GR4J (Perrin et al., 2003) and its snow routine CemaNeige (Valéry, 2010; Valéry et al., 2014), HBV (Bergström, 1992) or MORDOR (Garavaglia et al., 2017; Garçon, 1996) have proven their efficiency in mountain environments. However, few models specifically dedicated to karst environments include a snow routine, with the exception of luKARS (Bittner et al., 2018) and KarstMod (Mazzilli et al., 2019; Sivelle et al., 2025). To get past this limitation, some studies use precipitation preprocessing to estimate the snow component (Çallı et al., 2022; Fan et al., 2023; Hartmann et al., 2014). These approaches are generally based on degree-day models (Martinec and Rango, 1986), which estimate melting as a function of air temperature. Degree-days models are therefore parsimonious, using a conceptual modelling to simplify processes in the critical zone at the atmosphere-surface interface. They can be applied at the catchment scale, but give improved results with sub-catchments based on elevation bands (Ruelland, 2024; Valéry et al., 2014).

Despite these advances, the understanding of hydrological dynamic of mid-altitude karst systems influenced by snow remains limited. The strong spatio-temporal variability of snow accumulation and melting, coupled with the heterogeneity of karst aquifers, makes it difficult to accurately quantify discharge dynamics and their evolution under climate change. This knowledge gap is particularly critical in mountain regions where water resources are largely dependent on snow storage and release.

To address this issue, we take advantage of the latest release of the KarstMod platform (Sivelle et al., 2025) and apply a conceptual rainfall-snow-discharge model with a snow routine to investigate how snow accumulation and melt influence karst discharge, and how rising temperatures may affect flow regimes. The structure of the model is constrained by preliminary study of the hydrographs at the karst outlet and it incorporates a degree-day snow routine applied to sub-catchments.

The Dévoluy massif, located in the French Southern Alps, provides an ideal case study as it represents a typical mid-altitude karst system strongly influenced by seasonal snow, yet it has never been investigated through hydrological modelling before. Four years of daily observations of precipitation, air temperature and discharge are available. Altitudes range from 870 m a.s.l. at the Gillardes Spring – the main outlet – to 2789 m a.s.l. at the Obiou summit. Snowpack is generally present above 1700 m a.s.l. from December to April.

In this paper we aim to (1) simulate the daily discharge of a karstic spring in a mid-altitude Alpine environment using the KarstMod V3 platform (Sivelle et al., 2025), (2) quantify the contribution of seasonal snow to spring discharge, and (3) test the sensitivity of a karst spring discharge to temperature perturbations consistent with TRACC-type (climate change adaptation trajectory) warming trajectories (French TRACC guideline, Soubeyroux et al., 2024). The originality of this study is based on (i) a mid-altitude Alpine karst system with a strongly seasonal snowpack, (ii) the use of a parsimonious rainfall-snow-discharge conceptual model, (iii) a scenario design based on simple temperature perturbations, and (iv) a quantitative focus on seasonal metrics, focusing on the summer period.

2 Dévoluy case study and data

2.1 Case study: geology and hydrogeology

The Dévoluy massif is an isolated mountain range bordered to the west by the Trièves plateau and the Grand Buëch valley, to the north and east by the Drac valley, and to the south by the Petit Buëch valley (Fig. 1). It culminates at the Obiou (2789 m a.s.l.), and displays the characteristic morphology of a perched syncline, with a central depression located around 900 m a.s.l. surrounded by steep ridges over 2000 m a.s.l. The basin opens out to the north through the Souloise gorges, the lowest altitude of the massif at 870 m a.s.l., where the karstic Gillardes Spring is outflowing.

The karstified rocks forming the aquifer lie unconformably on a strongly folded and faulted substrate consisting of Jurassic and Lower Cretaceous limestones and calcareous marls (Ferry and Grosheny, 2013; Porthault, 1974, 1976). The Dévoluy karstic aquifer is developed within Upper Cretaceous (Campanian-Maastrichtian) calcareous turbidites, locally known as Senonian (Bonhomme, 1972; Luparini, 1975). This formation is composed of multi-decimetric finely bedded limestone interbedded with ribbed or nodular cherts. Their thickness varies here between 400 and 800 m (Parizot et al., 2025).

These Senonian deposits are themselves unconformably overlain by Eocene conglomerates, Eocene-Oligocene nummulitic limestones, marls and sandstones, and Quaternary units (Grosjean et al., 2017; Huet et al., 2023; Meckel et al., 1996; Monjuvent, 1979). While the carbonate deposits at the base of the Nummulitic are connected and contribute to the karst aquifer, the overlying marly limestone layers form semi-impermeable horizons that impact infiltration to the karstic aquifer.

A major structural feature, the Dévoluy median thrust, divides the massif into two main N–S synclines (Gidon and Pairis, 1976; Meckel et al., 1996; Parizot et al., 2025). To the east, the structure corresponds to the SE–NW-oriented Saint-Étienne syncline, while to the west it aligns with the SW–NE Saint-Disdier syncline. This double N–S organization overlies a broader E–W syncline whose limbs rise both southward toward the Pic de Bure (2709 m a.s.l.) and northward above the Gillardes Spring. This large-scale fold (∼10 km) produces converging north- and south-dipping strata toward the Saint-Disdier area, giving the Dévoluy its characteristic morphology of a northward-tilted basin. This configuration directs surface and groundwater flows toward the Souloise gorges, the main outlet of the massif (Fig. 1). The Dévoluy median thrust does not prevent hydraulic communication between the two synclines (Fig. 1), as showed by dye tracing experiments (Zappelli et al., 2018) and by recent 3D structural modelling of the Senonian karst reservoir (Parizot et al., 2025). The current groundwater drainage pattern results from successive phases of karst evolution during the Neogene and Quaternary (Rispal et al., 2025). Karstification is well developed in the massif, affecting both the Upper Cretaceous limestones and parts of the Eocene units. The oldest speleogenetic features, found around 2400 m a.s.l., date to the late Middle Miocene (Jagercikova et al., 2021; Lemot et al., 2023) whereas the main networks feeding the Gillardes Spring are likely younger. The present outlet is partly blocked by deposits from the last glaciation, and the downstream conduits remain inaccessible. Submerged passages identified below the spring (between 870 and 770 m a.s.l.) indicate a recent rise in the base level and point to a late-stage reorganization of the Dévoluy karst system (Audra and Palmer, 2013).

At lower elevation in the center of the massif, the presence of an Oligocene low-permeability marly cover above the karstic aquifer limits infiltration and allows the development of two rivers in each syncline branch: the Souloise river to the east and the Ribière river to the west. Both rivers converge at Saint-Disdier, before joining the Souloise gorges. At higher altitudes, there is relatively little soil and vegetation cover and the karst aquifer outcrops, thus limiting runoff and favoring fast and deep infiltration (Bonhomme, 1972; Luparini, 1975).

Numerous dye tracings have shown that almost all groundwater flows converge to the Gillardes Spring (Bonhomme, 1972; Luparini, 1975; Zappelli et al., 2018). Only to the south, a limited part of the Dévoluy massif, in the Bure plateau, flows southwards feeding several springs located in the Sigouste valley (Maire, 1980; Zappelli et al., 2018).

Located on both sides of the Souloise river at an altitude of 870 m a.s.l., the karstic Gillardes Spring is composed of two main outlets (the Petites Gillardes Spring and the Grandes Gillardes Spring) with a mean interannual discharge of 4.45 m³ s⁻¹ over the studied period. The catchment area drawn on Fig. 1 extends over 176 km²; it includes the watershed of the Souloise and Ribière rivers and the recharge area of the Gillardes Spring. The Bure Plateau is excluded because it is drained southward. The Dévoluy massif is a perched karst aquifer forming an isolated syncline (Fig. 1), topographically and geologically disconnected from adjacent massifs and aquifers. This specific setting ensures that the catchment boundary is well constrained with minimal uncertainty, and makes significant underground water exchanges with surrounding systems unlikely.

Figure 1 Geological map of the Dévoluy massif (SE France) and W–E structural cross section (Modified from Parizot et al., 2025). Hydrological information is also represented: rivers, Gillardes Spring, catchment limits, SAFRAN grid.

2.2 Datasets

In this mountainous area, to catch the spatial variability at the regional scale of the Dévoluy massif, we used meteorological data provided by the SAFRAN model of Météo-France (Vidal et al., 2010). The SAFRAN grids are 8×8 km long, with 4 meshes covering the study area (# 7481, 7482, 7565, and 7566). The average value of the 4 selected meshes was computed and used as daily input to the KarstMod model from 2015 to 2019. The meteorological inputs consist of precipitation (mm d⁻¹) and mean air temperature (°C). The average altitude of these four meshes is 1642 m a.s.l., close to the catchment area median altitude (1633 m a.s.l.). Precipitation and temperature data from the in-situ meteorological Saint-Étienne station (belonging to Électricité de France – EDF), located in Saint-Étienne-en-Dévoluy at 1290 m a.s.l., were used to check graphically the validity of the SAFRAN data. This dataset has been supplemented by in situ snow depth measurements from the mountain ski resort of Super Dévoluy, available from early November to late April and provided by Météo-France (Météo-France, 2025). Daily snow depth (in situ snowpack) and fresh snowfall are used as qualitative information on snow cover over the massif at the altitude of 1969 m a.s.l.

The karst groundwater discharging at the Gillardes Spring is calculated by subtraction between the upstream (Pont de la Baume) and downstream Souloise river (Infernet) gauging stations (belonging and managed by EDF). Water levels are recorded using a Panatronic pressure sensor at the upstream station and an OTT Compact Bubble Sensor at the downstream station. These sensors measure hydraulic head at calibrated gauging stations, which is then converted into discharge using rating curves. The rating curves were verified at least three times per year between 2015 and 2019 for the Pont de la Baume station, and once in 2017 for the Infernet station. Data are available at hourly interval from 12 May 2015 to 1 September 2019, and mean daily data are calculated.

2.3 Rain, snow, and discharge variability in the Dévoluy massif over the 4 study years

The analyses in this section are presented by hydrological year (from 1 September to 31 August of the following year) for the period covered by the karst discharge dataset, i.e. from September 2015 to 2019 (Fig. 2). To provide a climatic context, SAFRAN 30-year data (1990–2020) indicate an average annual precipitation of 1242 mm and 89 d yr⁻¹ with air temperature below 0 °C, resulting in a significant snow component over the Dévoluy massif, particularly above the mean SAFRAN reference altitude of 1642 m a.s.l. Over the 4 modelled years, mean annual precipitation is 1170 mm, detail of precipitation for each year is provided in Table S2 in the Supplement.

The interannual average air temperature over the 2015–2019 period was 6 °C, with a maximum monthly average of 15.5 °C in July and a minimum of −2.2 °C in January. The mean interannual snow is 33 % of total precipitation (391 mm, SAFRAN data). Annual snow amounts were highly variable, as illustrated by the in-situ snow depth measurement at the ski resort (Fig. 2b) and the SAFRAN data (Table S1). Most snowfall occurred between October and May, with a peak from December to March. The first two years recorded less than 300 mm of snow (SAFRAN), corresponding to 26 % and 24 % of total precipitation, values close to the minimum observed over the 30-year SAFRAN dataset (23 %). The third year (2017–2018) was particularly snowy, with 571 mm of snow, equivalent to 44 % of annual precipitation (the snowiest of the 30-year period). The fourth year was intermediate, with snow representing 38 % of annual precipitation.

The hydrological regime of the Gillardes Spring is typically alpine, marked by two periods of high water (Fig. 2a). The first, in autumn and early winter, corresponds to the largest rainfall events of the year (Fig. 2a). The second occurs in late spring, when snowmelt adds to the rainfall, maintaining a high discharge (e.g., spring 2018, Fig. 2a). By contrast, low water occurs both in winter, when precipitation falls mainly as snow, and in summer when the snowpack is depleted or melted, and when precipitation is rare and makes little contribution to recharge.

During spring, discharge at Gillardes Spring increases and remains high up to June (Fig. 2a). By contrast, the discharge of the Souloise river upstream the Gillardes Spring has a base flow starting to decrease already in April, reacting only to large rainfall events and causing a short surface runoff before returning to the base flow. Groundwater discharge predominates in the volumes of water measured at the outlet of the massif. On average, groundwater (i.e. Gillardes Spring) accounts for 86 % of the total discharge, while surface flow (i.e. Souloise river upstream the karst springs) accounts for only 14 %. From an annual water budget point of view, the previous percentage gives a rough estimate of the area of the watershed contributing to groundwater discharge. Considering that 86 % of the watershed recharges the groundwater outflowing at the Gillardes Spring, the surface of the karst recharge area is deduced to be around 150 km².

Figure 2 (a) Daily karst groundwater discharge at the Gillardes Spring, daily discharge of Souloise river at Pont de la Baume station, and daily precipitation from SAFRAN reanalysis. (b) In situ snow observations at the Super-Dévoluy snow monitoring station (1976 m a.s.l.): total snow depth and daily fresh snowfall. Grey shaded zones in panel (b) correspond to periods during which the station is closed.

3 Method: Rainfall-snow-discharge modelling

3.1 KarstMod

KarstMod is a lumped rainfall-discharge modelling platform dedicated to karst systems (Mazzilli et al., 2019, 2023, Sivelle et al., 2025), developed by the National KARST Observation Service (Jourde et al., 2018). KarstMod V.3.1.0 has been used in this study. This platform has been used to study the dynamics of karst systems on a variety of themes, e.g. water resources management (Baudement et al., 2017), understanding karst storage dynamics (Frank et al., 2021), modelling tool inter-comparison in karst hydrogeology (Jeannin et al., 2021), impact of climate change on karst flow dynamic and water resources (Fan et al., 2023).

This lumped model transforms daily or hourly precipitation (P) and evapotranspiration (ET) into simulated discharges (Q) via a system of interconnected compartments (also called reservoirs or buckets in the literature) reproducing the hydrological dynamic of the karst. KarstMod provides the option to compute potential evapotranspiration (PET) using Oudin's formulation (Oudin et al., 2005). The platform features an interactive interface that enables modelling to be progressively refined through iterations. The user establishes an initial model structure based on the knowledge of the system (conceptual model), then adjusts and adds complexity to this structure according to the results obtained (Mazzilli et al., 2023).

KarstMod reproduces conceptual models of karst systems by parameterizing different hydrological compartments and then simulating exchanges between them via an outflow equation. The model equations can be found in Mazzilli et al. (2023) and Sivelle et al. (2025). In classical configuration, flow can be linear or nonlinear.

The specific discharge from compartment A to compartment B (Q_AB) is defined by:

$$\begin{array}{}\text{(1)}& {\displaystyle }{Q}_{AB}=k_{AB}{\left({\displaystyle \frac{A}{L_{\mathrm{ref}}}}\right)}^{\propto AB}\mathrm{if}A\mathit{>}\mathrm{0}\text{(2)}& {\displaystyle }{Q}_{AB}=\mathrm{0}\mathrm{otherwise}\end{array}$$

where A is the water level in compartment A [Length] and k_AB is the specific discharge coefficient [L T⁻¹] for the discharge law from compartment A to compartment B (controlling the characteristic transfer rate between the two compartments: higher values leading to faster drainage and shorter residence times), L_ref [L] is a unit length and αAB [–] is a positive exponent (controlling the degree of nonlinearity and thus the sensitivity of discharge to storage variations : with larger values producing stronger flow responses under high-storage conditions). L_ref has been introduced in KarstMod to simplify the units when the law is not linear i.e. when αAB>1 (Mazzilli et al., 2023). For linear discharge laws, the specific discharge coefficient can be converted to a recession coefficient k (dimension [T⁻¹]), for L_ref the reference unit length fixed at 1 mm. In the following, k will be expressed in dimension d⁻¹. Assigning distinct coefficients thus enables to differentiate flow dynamics within the karst system, representing both fast and slow flows. The model's structure is based on two levels (Fig. 3): (i) a mandatory upper compartment (E), representing the karst infiltration zone (soil and epikarst), which receives the liquid precipitation (direct liquid precipitation and melting snowpack Psr) and constitutes the water reserve for evapotranspiration (ET) ; (ii) up to three optional lower compartments (L, M and C), simulating flows through the unsaturated and saturated zones and supplying the outlet (S). In addition, a direct flow from E to S can be activated to reproduce very fast transfers within the system. Added to this is an adjusted parameter E_min in E, which simulates a stored quantity of water only available for evapotranspiration, contained between zero and E_min (E_min<0). For modelling consistency, each compartment is also assigned at an initial water level corresponding to the water it contains at model initialization. Other possibilities exist and are not shown in Fig. 3, e.g., hysteresis, water exchange between M and C, and pumping, but they were not used in this study.

Modelling is divided into three stages, guided by the KarstMod GUI and corresponding to three time periods in the dataset: (1) Warm-up, to avoid any model initialization bias (Mazzilli et al., 2012). (2) Calibration, during which parameters are optimized to obtain the best match between simulated and observed discharge. (3) Validation, during which the set of parameters optimized following calibration is applied to the remaining time series. Parameter optimization is performed via a quasi-Monte-Carlo approach, using Sobol sampling to explore the parameter space. It is based on an objective function (W_obj) calculated for each parameter set. The user defines a minimum value of the objective function (W_{obj_min} value) to be satisfied in order to retain the parameter set, and the number n of parameter sets with W_obj > W_{obj_min} to be saved. The highest result gives one parameter sets for W_{obj_max}. Modeller can then use the best simulation defined by W_{obj_max} or the n best simulations for W_obj>W_{obj_min}. Finally, the platform integrates a set of graphical visualizations to facilitate evaluation of model performance and interpretation of results. Results can be displayed in KarstMod, or exported as CSV files for further work. New release of KarstMod includes a “fixed parameter mode” to upload the n previously calibrated parameter sets and run the model (n times) with a modified input data file to apply scenarios.

3.2 Snow routine

KarstMod integrates a snow routine that simulates the accumulation and melting of snowpack in one or more sub-catchments, on a daily or hourly time scale (Sivelle et al., 2025). This routine is based on a degree-day approach (Bergström, 1992). It requires air temperature (°C) as supplementary input data (precipitation is already a mandatory data in KarstMod). When temperature is below the temperature threshold T_s (°C), precipitation accumulates as snow. When the temperature exceeds T_s, melting occurs according to a melting coefficient MF (mm °C⁻¹ d⁻¹). The meltwater (P_sr) is fed into compartment E as soon as the water retention capacity of the snowpack (CWH) is exceeded. In addition, the refreezing coefficient (CFR) influences the refreezing of liquid water when the temperature falls below T_s. This module does not consider sublimation processes, nor slope orientation. Evaporation and evapotranspiration processes only occur in compartment E, as a lumped process, when water is available (Sivelle et al., 2025).

KarstMod is lumped, not a distributed model. However, the snow routine can be applied to sub-catchments dividing the catchment area into several altitudinal units, similar to the approach used in some hydrological and hydrogeological models such as GR4J with the Cemaneige module (Perrin et al., 2003; Valéry, 2010; Valéry et al., 2014), HBV (Bergström, 1976, 1992) or MORDOR (Garavaglia et al., 2017). This subdivision makes it possible to better represent the temperature gradient and to discriminate rain/snow distribution as a function of altitude. However, the total amount of precipitation applied on each sub-catchment is the same and stays independent of the altitude, i.e., without precipitation lapse rates.

In KarstMod, three parameters define these sub-catchments: the number of sub-catchments, the fraction of the total catchment area they occupy (value between 0 % and 100 %), and a temperature shift. This shift simulates the altitudinal variation in temperature between each sub-catchment.

3.3 Parameter estimation

Model performance is assessed by comparing simulated discharge with observed discharge. Several performance criteria are available in KarstMod (Mazzilli et al., 2023) to calculate the W_obj objective function. Here we use the Non-Parametric Kling-Gupta Efficiency criteria (KGEnp) (Pool et al., 2018) as suggested by Cinkus et al. (2023), which provides a robust evaluation of the agreement between observed and simulated discharge across a wide range of hydrological conditions, combining correlation, bias, and variability within a single composite metric. The calibration is performed using discharge observations only, and KGEnp is therefore applied as a composite discharge-based criterion.

KGEnp can vary between −∞ and 1. A value of 1 indicates a perfect match between observed and simulated discharge. In this study the optimization was performed by exploring the parameters using Sobol routine, until 5000 simulations that met KGEnp > 0.85 were obtained.

KarstMod offers other tools to assess the validity of the results, which are freely interpreted by the modeller: graphical comparison of observed and simulated discharge time series curves, cumulated volume discharged at the spring, internal water levels in each compartment, simple and crossed rainfall-discharge correlograms, distribution of classified discharges. This validation by an expert is part of the modelling procedure. In the specific case of snow-influenced mid-altitude karstic hydro systems, beyond the performance criterion for the entire simulated series (KGEnp), the modeller will seek to validate the consistent consideration of snowpack evolution: (1) cold winter period generates snow stock and liquid precipitation (precipitation from the snow-routine P_sr) becomes lower than total precipitation; (2) in spring, snowmelt is added to precipitation, leading to a sharp increase in groundwater recharge, (3) for the rest of the year, the influence of snow on the hydrology of the system is considered marginal as there are no glacier or snow in mid-altitude areas.

Figure 3 Structure of the KarstMod rainfall-snow-discharge model designed for the Dévoluy case study. Snow routine converts rain into snow and vice versa, and gives input liquid precipitation (P_sr). E, M and C are compartments for flow transfer to the Spring (S). Q_ES is a direct flow from E to the spring.

3.4 Modelling strategy

The modelling strategy in this study comprises three steps: (1) calibration and validation of a simple model without the snow routine (Model 1); (2) calibration and validation of the model including the snow routine (Model 2); and (3) sensitivity analysis of spring discharge to air temperature warming, performed using the previously calibrated Model 2 for two scenarios with modified air temperature (Model 2 T+2 and T+4). These steps and the associated models are described below and summarized in Fig. 4.

The calibration was run on the 2015–2018 period, encompassing contrasted snow cover conditions, while validation was performed over a final hydrological year (1 September 2018 to 31 August 2019), characterized by intermediate snow conditions. Initially, we set W_{obj_min}= 0.5 for the first steps of modelling, and then increased it to W_{obj_min}= 0.85, retaining the 5000 simulations with W_obj>W_{obj_min} to ensure model robustness. For each day, the 5000 results (i.e. 5000 simulated daily discharge) are summarized by their median and 5th–95th percentile in graphs. The specific result for W_{obj_max} is also highlighted. When discussing the results obtained from the 5000 simulations, median absolute deviation (MAD) was calculated to get an interval around the median. To assess the temporal robustness of the 5000 parameter sets, split sampling tests were performed by inverting the calibration and validation periods (validation: September 2015 to September 2016; calibration: September 2016 to September 2019), yielding comparable median W_obj KGENP values of 0.74 vs. 0.75 in both configurations, supporting the temporal robustness of the calibrated parameter sets.

To evaluate model performance, four complementary criteria were selected: Nash-Sutcliffe Efficiency (NSE) for overall fit, Volume Error (VE) for temporal distribution of errors without compensation between positive and negative values, Mean Absolute Error (MAE) for average error magnitude, Percent Bias (PBIAS) for systematic bias direction, and non-parametric Kling-Gupta Efficiency (KGEnp) used for calibration is kept for integrated assessment of correlation, variability, and bias. This combination allows comprehensive evaluation of both global performance and specific aspects to correctly assess snowmelt and spring discharge simulation.

Figure 4 Schematic representation of the modelling and simulation framework. The graphs show observed and simulated discharges of the Gillardes Spring for hydrological year 2017–2018. Model 1 uses a E–M–C–ES KarstMod structure, without snow routine. Model 2 uses a E-M-C-ES KarstMod structure, with snow routine. Two scenarios are then simulated with Model 2 and results displayed by seven-day average: Model 2 T+2 and Model 2 T+4 with a temperature shift applied to input dataset (+2 and +4 °C, respectively).

3.4.1 Model without snow routine (Model 1)

The goal of the first step is to run a basic KarstMod model without snow routine (Fig. 4), to select the compartments and flow paths representing the conceptual hydrogeological functioning of the Dévoluy massif. A classical E-M-C model structure with one upper compartment and two lower compartments is set up. Groundwater is divided in fast flow through C and slow flow through M. A direct flow from E to S (ES) is implemented to better represent very fast flow observed on the discharge data. The model structure is determined on the basis of a preliminary analysis of the hydrograph, which allows dividing the recession curve into successive phases. On a semi-logarithmic plot of discharge versus time, these recession phases appear as distinct linear segments with different slopes, each representing a specific drainage component of the karst system. Three linear segments are recognized and their regression slopes give three recession coefficients used as an initial estimate of the characteristic time for the three conceptual flows through compartments discharging to the spring in KarstMod. For the automatic calibration, these three values define a pre-constrained parameter range based on expert knowledge, within which the optimization is run to search for the best-fit values. This ensures parameters remain within physically plausible bounds while allowing sufficient exploration of the parameter space.

E_min is implemented but left free over a wide parameter range in order to allow the model to adjust actual evapotranspiration. KarstMod potential evapotranspiration routine is activated to calculate PET with Oudin's formula. Input data are total precipitation and daily temperature data from the SAFRAN reanalysis as described in the previous section. For this model 1, snow routine is not activated.

The parameter surface of the catchment area (R_A in KarstMod GUI) requires special attention. It is included in parameters to limit precipitation misestimation in mountain environment. In high elevation areas, SAFRAN is known to provide an underestimate precipitation (Gottardi et al., 2012; Ruelland, 2023). The precipitation underestimation will be corrected in the KarstMod rainfall-snow-discharge modelling by an increase in surface area in order to balance the water budget. Even though the actual karstic catchment area has been estimated at 150 km² based on reliable hydrogeological criteria, the parameter range is set from 150 to 250 km² in the lumped model.

The model requires the calibration of ten parameters (Table 1): five recession coefficients (three for compartment E, one for compartment M, and one for compartment C), parameter E_min, the catchment surface area, and three parameters corresponding to the initial levels of the compartments at the beginning of the warm-up period.

3.4.2 Model with snow routine (Model 2)

Keeping the structure of the initial model (E-M-C-ES) with PET routine, we activated the snow routine with sub-catchments in Model 2 and processed a new calibration. The chosen E-M-C-ES structure was retained as the most appropriate representation of the system based on the previous graphical analysis of the hydrograph. We chose not to add any additional complexity to the model (Hartmann et al., 2012) in order to avoid overparameterization. However, the uniform distribution of snow across the Dévoluy massif being too rough an approximation, sub-catchments are introduced to better discretize snow storage and melting zones. This is justified by the wide elevation range of the catchment, from 870 m a.s.l. at the spring to 2789 m a.s.l. on the summits.

The configuration of sub-catchments for the KarstMod snow routine is based on a hypsometric analysis performed using the 1 m LiDAR DTM (IGN, 2024). Elevation bands were defined with equal area, in line with the recommendation of Valéry et al. (2014). They also recommended using 5 bands, but testing from 1 to 5 bands showed no improvement beyond 4, so four elevation bands were retained for parsimony, each representing 25 % of the total catchment area (Fig. 5). Their boundaries are further supported by observed winter snow cover patterns in the Dévoluy massif.

The SB1 sub-catchment (870–1414 m a.s.l.) corresponds to the least snow-covered zone, with rare snowfalls and a short-lived snow cover. SB2 (1414–1634 m a.s.l.) covers the lower altitudes of the ski resort, with regular snowfalls but little persistent snow cover. SB3 (1634–1980 m a.s.l.) marks the start of the quasi-permanent winter snowpack. Finally, SB4 (1980–2789 m a.s.l.) groups together the highest altitudes, where snow cover is continuous from November to May, with significant snow accumulation. A temperature shift is applied to each sub-catchment to model the altitudinal gradient. The shift is calculated by multiplying the difference in altitude between the SAFRAN measurement point (1642 m a.s.l.) and the median altitude of the sub-catchment by a regional gradient of 0.6 °C per 100 m. This regional gradient is a rough estimate, close to the one proposed by Ruelland (2023), but the model will compensate for this approximation by adjusting the T_s parameter (temperature threshold for melting), to improve the distribution of snow and rain over the catchment. It is negative with altitude higher than the reference and positive with when altitude is lower than the reference. The shifts obtained are +2.0 °C for SB1, +0.8 °C for SB2, −1.0 °C for SB3 and −3.3 °C for SB4. Range of the snow routine parameters are defined from the literature (Ismail et al., 2023; Lafaysse et al., 2017; Sivelle et al., 2025) and are: −2 to +3 °C for T_s, 3 to 8 mm ° d⁻¹ for MF, 0.3 to 0.8 for CFR and 0.05 to 0.25 for CWH (Table 1).

Figure 5 Snow routine sub-catchment parametrization showing four elevation bands representing the altitude distribution over the Dévoluy catchment. Each sub-catchment (SB1–SB4) is defined by its lower and upper altitude limits and the corresponding temperature shift. The blue diamond indicates the SAFRAN input temperatures, and the orange diamonds show the altitude used for the temperature adjustment in each band.

3.4.3 Scenario simulations: altered temperature and snow regimes

We examine the effect of a reduction in snowfall and modification of snowmelt on simulated discharge caused by air temperature increase (Fig. 4). Historical data and global and regional climate models already show warming in France in 2020 between +1.5 and +2 °C in reference to pre-industrial period (e.g. Ribes et al., 2022). RCP 4.5 or 8.5 or equivalent SSP scenarios (Calvin et al., 2023; Ribes et al., 2022) predict a global warming of more than +4 °C before the end of the century. In France, the TRACC provides a reference warming trajectory for the 21st century, including climate indicators such as temperature, precipitation, and extremes, to guide adaptation strategies across multiple sectors (Soubeyroux et al., 2024).

In this study we use the TRACC temperature recommendations to perform a temperature sensitivity analysis, keeping precipitation unchanged. Two temperature perturbations are applied to the calibrated Model 2: observed temperature in 2015-2019 shifted by +2 °C (Model 2 T+2), which is equivalent to +4 °C with respect to 1850–1900, and a shift of +4 °C (equivalent to +6 °C with respect to 1850–1900) to illustrate extreme warming of 2 °C over the mean (Model 2 T+4). This type of warm year is expected in a future climate, taking into account climate internal variability, model uncertainty, and scenario uncertainty (Evin et al., 2021). This approach isolates the impact of warming on snow accumulation, melt and recharge timing without introducing additional uncertainty from precipitation changes.

Table 1 Calibrated parameters and objective function values for Model 1 and Model 2. Values correspond to the median of the 5000 best simulations. The distribution of the parameter sets and the range applied to these parameters are shown in Fig. S1 (Supplement).

4 Results

4.1 Simulation of the Gillardes Spring discharge

4.1.1 Model 1: without snow routine

Calibrated parameters of the 5000 best simulations for Model 1 are given in Table 1, distribution of these parameter is shown in Fig. 11 in supplementary material. This first model, calibrated without snow routine, reproduces the overall discharge dynamics of the Dévoluy karst system but shows important discrepancies with observations in key periods (Fig. 6a). Simulated flows are overestimated during winter months (December to February), when most precipitation falls in the Dévoluy massif as snow and may not immediately contribute to recharge. Conversely, simulated spring discharges (April to early July) are underestimated, with a marked deficit in May and June when observed flows usually reach their seasonal peak. This misfit is particularly visible in May–June 2018 during a very snowy year. These mismatches generate a temporal shift in the simulated hydrograph, with excessive winter discharge and insufficient spring contributions. This pattern is reflected in the model performance, with a median KGEnp value for the 5000 simulations of 0.60 over the calibration period and 0.54 in validation, indicating only moderate ability to reproduce both the timing and magnitude of flows (Table 1). Performance metrics in Table S4 confirm that Model 1, without a snow routine, is inadequate for this catchment. Overall performance is poor (2015–2019: NSE = −0.18 MAE = 3.34 m³ s⁻¹) and this is confirmed during snowmelt periods (April–August: NSE = 0.01, PBIAS =  47.4 %). This systematic 47 % underestimation of spring discharge demonstrates the model's inability to capture snowmelt-driven discharge. Annual metrics show consistent deficiencies across all four years, with snowmelt period PBIAS ranging from −28.6 % to −57.4 %, confirming the necessity of explicit snow representation for this mid-altitude karst system.

4.1.2 Model 2: snow routine activated

Calibrated parameters of the 5000 best simulations for Model 2 are given in Table 1, distribution of these parameters is shown in Fig. 11 in supplementary material. Here the snow routine is activated and the simulation improves substantially as illustrated by the median curve and the narrow envelope of 5th–95th percentile (Fig. 6b). The curve for W_{obj_max} is close to the median result. Precipitation occurring during cold periods are temporarily stored as snow and progressively released during the melting season thanks to the snow routine, which results in a more realistic redistribution of recharge. The spring flow regime is much better reproduced, with the model successfully capturing both the timing and the magnitude of the discharge peak between April and June. The excessive winter recharge previously simulated by the rainfall-only Model 1 is also corrected, bringing simulated flows closer to observed values. The Fig. 6c shows the simulated snowpack for the sub-catchments SB3 and SB4 (the highest elevation bands) and the in situ measured snow depth. The simulated snowpack with snow routine represents the water equivalent of the snowpack rather than its actual thickness, as snow density changes over time. The in-situ measurements, taken at 1969 m (between the elevation ranges of sub-catchments 3 and 4), fall well within the envelope formed by the simulated curves of sub-catchments 3 and 4, indicating that the model captures the overall evolution of the snow cover at this location.

Model efficiency has a good KGEnp median for the 5000 simulations of 0.86 over the calibration period and 0.74 in validation (Table 1). Calculated KGEnp year by year (Table S4) ranges between 0.74 and 0.89 (0.85 for the mean value over the 4 studied years). Other performance metrics in Table S4 demonstrate the effectiveness of the snow routine implementation. The model achieves good overall performance (2015–2019: NSE = 0.70, VE = 0.63, PBIAS = 4.9 %, MAE = 1.65 m³ s⁻¹). During snowmelt periods (April-August), performance is good (NSE = 0.77, VE = 0.71), demonstrating the model's ability to capture spring discharge dynamics. Annual performance varies with snow regime characteristics (Table S1): NSE ranges from 0.30 to 0.87, with the best performance occurring in 2017–2018, the year with highest snowfall contribution (571 mm, 44 % of total precipitation). For each snowmelt periods (April–August), NSE ranges from 0.26 to 0.88 across the four years, confirming that the snow routine enables accurate simulation of discharge dynamics in this mid-altitude karst catchment.

The time series of internal water level in KarstMod compartments, and internal discharge (plotted in Fig. 7 with seven-day averaged values), show how KarstMod distributes precipitation (P_sr) and discharges over time from E to the outlet (S, Spring). Negative water level values in E stand for soil and epikarst water storage available for actual evapotranspiration between E_min and zero. It corresponds to summertime, when ETP is maximum, snowpack has almost completely melted, and precipitation is scarce in this region. Direct flow from E to S (QES) closely follows the water level in E, given quick discharge (high recession coefficient k_ES, Table 1). Internal water levels in compartments M and C vary according to their respective recession coefficient. Compartment M acts for baseflow, with large inertia, discharging to the Spring for months after the recharge period feed by rainfall and snowmelt. Compartment C, by its high storage (high water level) and high dynamic response (high recession coefficient); accounts for discharge during high-water flow in spring season influenced by snowmelt but also during peak events from rainfall or snowmelt during the year.

Overall, these results highlight the importance of explicitly accounting for snow processes in the Dévoluy karst catchment. By integrating snow storage and melting, the model not only aligns more closely with observed hydrological patterns but also provides a higher value of performance criteria, which strengthens its reliability for further scenario analysis.

Figure 6 Observed and simulated discharge for the Dévoluy karst system (2015–2019) for (a) Model 1 and (b) Model 2. Shaded areas show the 5th–95th percentile of the 5000-parameter sets, solid lines the ensemble median, and dashed lines the W_{obj_max} simulation. (c) Model 2 simulated snow water equivalent for subcatchments 3 and 4 (shaded: 5th–95th percentile) and observed snow depth at the Super-Dévoluy station (grey bars); gaps indicate periods with no available measurements.

Figure 7 (a) Median of 5000 simulations of seven-day averaged water levels in compartments E, M, and C (Model 2). (b) Median of 5000 simulations of seven-day averaged internal discharges (QES: discharge from E to S, QMS: discharge from M to S, QCS: discharge from C to S). Seven-day cumulative precipitation is shown on a secondary y-axis (inverted).

4.2 Impact of temperature increase: T+2 and T+4 scenarios

This section addresses the question: what would have been the flow rate of the Gillardes Spring if the rainfall/snowfall proportion from 2015 to 2019 had been affected by an increase in air temperature?

For each discharge value listed below, a range comprising a lower and upper limit (±) is provided. This corresponds to the median absolute deviation (MAD) which allows to better represent the potential variability in simulated discharge for the 5000 parameter sets.

For the +2 °C scenario (Model 2 T+2), simulated discharge starts to exceed that of the reference Model 2 during late autumn and early winter. This difference persists until mid-spring, then the two curves (Fig. 8a) intersect (Model 2 T+2 becomes lower than reference). From July onwards, the curves converge and seem to overlap throughout the summer low-flow period, but the Model 2 T+2 stays below with a lower baseflow. Over the 4 modelled years, the mean simulated discharge decreases from 4.66±0.13 m³ s⁻¹ for Model 2 to 4.41±0.12 m³ s⁻¹ for Model 2 T+2, corresponding to a reduction of approximately 5 % in simulated discharge. For the +4 °C scenario (T+4), the temporal pattern remains broadly similar, but with accentuated differences between the two curves (Fig. 8b): autumn and early winter precipitation gives higher discharge at the Gillardes Spring for increasing air temperature, and lower discharge during spring and summer. For T+4 scenario, spring recession starts earlier. The differences are more pronounced for the hydrological year 2017–2018, which was characterized by heavy snowfall and a substantial seasonal snowpack (Table S1). For the 4 modelled years, the mean simulated discharge decreases from 4.66±0.13 m³ s⁻¹ for Model 2 to 4.14±0.12 m³ s⁻¹ for Model 2 T+4, corresponding to a reduction of approximately 11 % in simulated discharge (Table S2).

Snow pack analysis at the elevation of in situ snow measurements (sub-catchment 3) confirms this shift: the median number of days with snow (> 1 mm snowpack solid component) decreases from 154 d yr⁻¹ in Model 2 to 99 d in T+2 and 50 d in T+4 (−67 %). In 2017–2018, the year with highest snowfall (571 mm, 44 % of precipitation), snow cover duration drops from 186 days to 74 d under T+4 (−60 %), with maximum snowpack solid component reduced by 84 % (from 742 to 118 mm). This indicates that while snow cover duration tends to diminish under warming, it is not entirely eliminated, particularly during mid-winter cold periods.

Moreover, average AET for the 4 years studied, calculated for Models 2, 2 T+2 and 2 T+4 (Table S3) shows a AET increase: respectively 36 % of annual precipitation, 40 % and 44 %. This value found in Model 2 is in line with other studies showing that in mountainous catchment around a third of precipitation goes to atmosphere via evapotranspiration (Fan et al., 2023; Herrnegger et al., 2012).

Figure 8 Seven-day average simulated discharge of the Gillardes Spring over the 2015–2019 period. (a) Comparison of Model 2 and Model 2 T+2 (temperature warming T+2 °C). (b) Comparison of Model 2 and Model 2 T+4 (temperature warming T+4 °C). Shaded areas represent the 5th–95th percentile envelope of the 5000-parameter sets. Solid lines show the ensemble median and dashed lines the simulation with the parameter set corresponding to the maximum objective function value (W_{obj_max}).

5 Discussion

5.1 Insights from the rainfall-snow-discharge model on the mid-altitude Dévoluy karst aquifer

In this study, we proposed a modelling approach based on a model that considers the effects of snow on the recharge of a mountain karst aquifer. Modelling results were greatly improved by considering snow recharge using the KarstMod snow routine, which is an essential component of the Dévoluy karst system (KGEnp median value of the 5000 simulations = 0.85). This kind of routine, widely used in hydrological and hydrogeological models (Bittner et al., 2018; Seibert and Vis, 2012; Valéry et al., 2014), is based solely on precipitation and temperature data and is therefore particularly suited to a parsimonious approach, unlike energy-balance-based models, which are more complex and require additional inputs (Franz et al., 2008).

The optimal calibration of the model was achieved with an E-M-C structure enriched with an additional link between E and S (Q_ES), revealing three contrasting response components. Compartment M gives a slow discharge that supports baseflow. Direct flow from E to the Spring corresponds to very fast discharge reflecting immediate response to recharge events. Compartment C gives an intermediate hydrodynamic response, i.e., fast flow lasting several weeks.

Model 1 without snow routine was not validated because its overall performance criterion (KGEnp) remained below 0.60 and karst discharge was poorly simulated during winter and spring periods, which are typically affected by snowmelt dynamics in mid-altitude mountain catchments. Model 2, with snow routine, was validated and shows during snowy year (Fig. 6b) that observed high-water discharge during spring season is highly supplied by delayed snowmelt recharge. This is confirmed by the years with low snowfalls, during which observed discharge in spring season remains low. Snow thus plays a dual role by shifting the timing of recharge and by modulating its dynamics, as meltwater is progressively released over several weeks. In contrast, liquid rainfall generates short recharge pulses that lead to rapid discharge increases at the Gillardes Spring. The studied area has a high infiltration capacity, characterized by a karstic surface morphology where the rocks are heavily fractured and covered with karren, which fosters infiltration and reduces runoff. The Souloise River watershed (Fig. 1) largely corresponds to the groundwater recharge area. Only 14 % of the water in the Souloise downstream of the hydrosystem comes from surface river. Most of the precipitation does not run off to the stream; instead, it infiltrates and recharges the karst aquifer of the Gillardes Spring. The rapid flows in Model 2 (from E to the spring, and from E to C to the spring) show that the strong karstification observed in the unsaturated zone is used by recharge water for rapid vertical transfer to the aquifer, then through the saturated zone to the spring. These rapid flows were also verified by artificial tracers (Zappelli et al., 2018). The good agreement between simulated and observed summer low flows over 4 hydrological years, including both wet and dry as well as snow-rich and snow-poor years, validates the use of the recession coefficient of compartment M as a robust indicator of baseflow recession ($\mathrm{8.94}\times {\mathrm{10}}^{-\mathrm{3}}$ d⁻¹). This parameter provides a quantitative metric of the aquifer's functioning and allows the simulation of groundwater discharge depletion during snow-free summer periods. This is a typical value of slow flows, but it is high enough to allow groundwater drainage above spring's base level within a few months once recharge ceases. The quantification of rapid and slow flows obtained using Model 2 shows that in the Dévoluy case study Gillardes Spring discharge is highly dependent on recharge distributed over time: either from snowmelt or from successive rainfall events throughout the year.

5.2 Evolution of recharge and discharge under temperature variations: a warning for stakeholders

The snow routine relies on precipitation and temperature input data, with temperature playing a key role as it controls the distribution between rainfall and snowfall, as well as snow accumulation and melting, which in turn affect the timing and magnitude of recharge. By modifying only input temperatures, we can simulate warming scenarios inspired by climate projections for the Alps, which anticipate a rise in temperatures that significantly alters the seasonal distribution of snow (Christensen et al., 2007). This approach isolates the effect of temperature on snow/rain partitioning without introducing additional uncertainties linked to precipitation (Duethmann et al., 2020). These scenarios are not intended as predictive forecasts, but rather as a temperature-sensitivity test to explore how the system responds to temperature forcing in a hypothetical case where precipitation remains as observed between 2015 and 2019.

Mean simulated discharge distribution at seasonal scale and monthly scale are given in Figs. 9 and 10 for the 5000 simulations. When discussing the results for a specific month or season, the median value is given, more or less the median absolute deviation (MAD).

In winter, simulated discharge increase in the warming scenarios by + 43 % (T+2 °C) and +72 % (T+4 °C), with median discharge rising respectively from 3.37±0.25 m³ s⁻¹ in the reference scenario (Model 2) to 4.82±0.22 m³ s⁻¹ and 5.81±0.21 m³ s⁻¹. This is explained by enhanced liquid water availability due to diminished snow accumulation (direct rainfall or short snowfall event and rapid snowmelt). Discharge during the spring season decreases under warming scenarios by 11 % (T+2 °C) and 30 % (T+4 °C), with a median of 6.67±0.24 m³ s⁻¹ in the T+2 °C scenario and 5.24±0.19 m³ s⁻¹ in the T+4 °C, compared with 7.5±0.30 m³ s⁻¹ in the reference scenario. This reduction is consistent with reduced snow accumulation during winter, resulting in lower spring discharge. This pattern, observed in the seasonal mean discharge, is confirmed when using monthly average discharge over the 4 simulated years (Fig. 10), with the T+2 °C and T+4 °C scenarios showing markedly lower discharges than both the observed values and the simulated discharge in the reference Model 2 in April, May and June. During the summer season, median simulated discharge decreases as the temperature scenarios increase: it drops by −48 % (T+2 °C) and −61 % (T+4 °C) in reference to Model 2, decreasing respectively from 3.62±0.35 m³ s⁻¹ to 1.89±0.20 m³ s⁻¹ and further to 1.40±0.16 m³ s⁻¹. This decline illustrates the progressive reduction in water availability during summer under warmer conditions. Focusing to August (Fig. 10), when water demand is high and rainfall and Gillardes Spring discharge is low, mean monthly simulated discharge (over the 4 years studied) decreases from 1.47±0.24 m³ s⁻¹ in Model 2 (and 1.30 m³ s⁻¹ in observed data) to 1.01±0.21 m³ s⁻¹ in Model 2 T+2 (−26 % vs. Model 2) and 0.85±0.19 m³ s⁻¹ (−42 % vs. Model 2) in Model 2 T+4 (Fig. 10). It is important to note that here we compare simulated discharge over the 4 modelled years, rather than focusing on a specific year (e.g., a particularly snowy or dry year). This is because observed discharge depends not only on temperature but also on many other factors, such as total annual precipitation and its seasonal distribution. Here, we focus on the response of the system to temperature changes, and thus on the overall, long-term behavior of the catchment. Furthermore, karst discharge decrease will be even greater in a year when spring precipitation is scarce, leading to an early low baseflow of the karst spring.

In autumn, precipitation observed at the highest elevations of the catchment is partly snowfall, but its timing is highly variable between years (Fig. 2). During early autumn (September to mid/late October) corresponding to the end of the summer low-flow period, Model 2 produces slightly higher flows in September (1.32±0.16 m³ s⁻¹) than Model 2 T+2 and Model 2 T+4 (0.97±0.18 and 0.77±0.18 m³ s⁻¹), whereas from November onward, especially during the first snow events, flows with T+2 °C and T+4 °C become higher (monthly median jumping from 6.98±0.31 m³ s⁻¹ in Model 2 to 8±0.36 m³ s⁻¹ in both Models T+2 and T+4). These differences are not sufficient to substantially modify the distribution of values shown in the violin plot, but they reflect the loss of the current tipping point in November from a rainfall-dominated to a snow-dominated regime. Overall, these results illustrate a clear shift in seasonal flow patterns under warming scenarios, with a transition toward rainfall-driven recharge: earlier and higher winter discharges, reduced spring high-water period, and mean summer discharge drastically reduced. Moreover, temperature increase impacts evapotranspiration (ET). Potential ET is currently correlated with temperature in Oudin's formula, generating PET increase for scenarios T+2 and T+4 °C. Median Actual ET calculated by KarstMod water budget in compartment E for the 5000 simulations also increases in average from 427±16 to 468±16 mm yr⁻¹ (+9 %) and 512±18 mm yr⁻¹ (+16 %), respectively for Model 2, Model 2 T+2, Model 2 T+4, indicating increased water stress and reduced water availability for discharge (Table S3).

The uncertainty quantified through the distribution of 5000 calibrated parameter sets, varies across seasons and scenarios. Expressed as percentage of the median discharge, the percentage of seasonal MAD ranges across the three models from 3.6 % to 7.4 % in winter, 3.6 % to 4.0 % in spring, and 9.1 % to 13.6 % in summer. Despite this variability, the 5th–95th percentile ranges (Fig. 9) show consistent trends across all 5000 simulations: winter discharge systematically increases under warming, spring discharge decreases, and summer discharge shows the most important reduction. This demonstrates that while parametric uncertainty affects the magnitude of change, the direction of seasonal redistribution under warming is consistent across the entire parameter ensemble.

These changes have direct implications for water resource management, as they alter the availability of water during critical periods for agriculture, industry, and drinking water supply.

While these scenarios isolate temperature effects, a fully prospective study using climate projection data could extend this analysis by incorporating projected precipitation changes and simulating longer periods to assess multi-year drought impacts, such as consecutive dry years predicted by climate models (Calvin et al., 2023). Testing these projections with multiple hydrological models (e.g., ensemble approaches like EXPLORE2; (Marson et al., 2024) would further reduce uncertainties and provide a wide range of possible future scenarios to be considered for water resource management.

Figure 9 Mean seasonal discharge distributions for the Dévoluy karst system over the 2015–2019 period according to scenario simulation described in Fig. 4. For each season, the star symbol indicates the observed mean seasonal discharge. Violin plots show the distribution of mean seasonal discharge across the 5000 parameter sets for Model 2 (gold), Model 2 T+2 °C (blue), and Model 2 T+4 °C (red). The vertical bar and white dot indicate the 5th–95th percentile range and the median of the 5000 parameter sets.

Figure 10 Mean monthly discharge distributions for the Dévoluy karst system over the 2015–2019 period, shown for each calendar month (September to August) according to scenario simulation described in Fig. 4. For each month, the star symbol indicates the observed mean monthly discharge. Violin plots show the distribution of mean monthly discharge across the 5000 parameters set for Model 2 (gold), Model 2 T+2 °C (blue), and Model 2 T+4 °C (red). The vertical bar and white dot indicate the 5th–95th percentile range and the median of the 5000 parameter sets.

5.3 Addressing water balance uncertainties in mountain catchments

As discussed above, snow plays a central role in the hydrological cycle, acting both as a temporary storage of winter precipitation and as the main source of spring and early summer discharge. An accurate representation of snow accumulation and melting is therefore essential to adjust the water balance. However, uncertainties in snowfall measurement and spatial distribution often lead to an imbalance between simulated recharge and observed discharge.

It should be noted that a good hydrological model will correctly reproduce variations of karst spring discharge or streamflow over time, but it must also give a balanced water budget. A balanced water budget is based on a good estimate of actual evapotranspiration and input precipitation data, representative of the catchment. Some hydrological models have an adjustment parameter to balance the water budget, such as the catchment water exchange coefficient in GR models (for instance in GR4J model – Perrin et al., 2003), or the precipitation correction factor in MORDOR model (e.g. Evin et al., 2024). Ruelland (2023) proposed a snowfall adjustment parameter for mountain catchments. In KarstMod, water budget adjustment can be done by two methods: adjusting the recharge area (RA) or using Q_loss, which simulates a discharge that leaves compartment E but does not reach S, the outlet of the modelled basin.

In a study case where PET and precipitation are well constrained, the parameter RA can be set by the user to the known hydrogeological catchment area. In the opposite case, this parameter is left free, within a range which may exceed the actual size of the studied catchment. The model will adjust the RA value automatically, and typically there are two possible scenarios: (1) the value obtained is indicative of the likely surface area of the actual catchment, (2) the value obtained is greater than the geologically acceptable value. In this second case, and if there are no human disturbances on the water cycle in the study catchment, this means that the precipitation given as input is insufficient to balance the water budget, i.e., there is more water discharging at the karst spring than is being recharged in the known catchment. Additionally, this imbalance may be true if there is a transfer of water between hydrological catchments. But if there is no water transfer between catchments, this imbalance may also illustrate underestimated precipitation input data. This point has been systematically highlighted in hydrological studies of mountain catchments in Europe (Evin et al., 2024; Gottardi et al., 2012; Pulka et al., 2024; Ruelland, 2024). Field rain gauge data and reanalyzed SAFRAN data underestimate precipitation. As highlighted by Ruelland (2023), precipitation may be insufficiently represented within the mountain catchments because snowfall is under-captured (depending on the precipitation gauges and their wind shields), and because precipitation is subject to orographic gradients (depending on the atmospheric conditions). Whatever the underlying reason, Ruelland (2023) concluded that the snowfall adjustment parameter allows total input precipitation to be increased, which clearly helps to fulfill the water balance equation and to more realistically simulate both snow and runoff.

Optimal calibration of the KarstMod model yielded a median catchment area (RA) over the 5000 simulations equal to 198±7 km², whereas the actual area is estimated at around 150 km². The RA parameter was used here as an adjustment parameter. The Dévoluy case is geologically interesting because the catchment is perched, with valleys that isolate the aquifer all around the massif. We can therefore deduce that the actual size of the catchment is well constrained, and moreover, underground exchanges between watersheds are non-existent. The RA calibrated value is almost 1,3 times higher than the expected value, which would correspond to an increase of 32 % in total precipitation over the Dévoluy massif and confirm the underestimation of SAFRAN precipitation in this mountainous catchment.

Another potential source of imbalance in mountain water budgets is snow sublimation. Although sublimation can represent a significant loss in high-elevation environments, recent studies suggest its relative contribution decreases under climate warming. For instance, Sexstone et al. (2018), in a Colorado Rocky Mountain case study, show that snow sublimation rates corresponding to climate warming simulations remained unchanged or slightly increased, but total sublimation losses decreased by up to 6 % because of a reduction in snow covered area and duration. In the case of Dévoluy, sublimation is not explicitly taken into account in the lumped rain-snow-discharge model, and therefore has no impact on the water budget. Ruelland (2023) had also shown that a lumped hydrological modelling on a dataset covering 17 mountainous catchments in the French Alps and Pyrenees gives good results without including a sublimation term in the snow routine. However, the reduction in snow-covered area and duration increases the actual evapotranspiration of non-snow-covered areas, which is taken into account by KarstMod in compartment E. However substantial uncertainties remain about future ET evolution. Lemaitre-Basset et al. (2022a) indicate that elevated atmospheric CO₂ may reduce ET via stomatal effects, potentially counteracting temperature-driven increases in PET. In addition, climate warming is expected to induce significant phenological changes, such as earlier leaf-out and a lengthening of the growing season (Menzel and Fabian, 1999; Richardson et al., 2013), which can alter both the magnitude and seasonal distribution of evapotranspiration. An earlier onset and extended duration of vegetation activity may enhance transpiration over longer periods, thereby modifying the seasonal partitioning between evapotranspiration and recharge. These ecohydrological processes are not explicitly represented in the present modelling framework, and their effects may either amplify or offset temperature-driven changes in evapotranspiration. Thus, future evapotranspiration trajectories are likely to depend not only on temperature and atmospheric CO₂, but also on vegetation dynamics and phenological shifts. In a complementary study, Lemaitre-Basset et al. (2022b) shows that the choice of PET formulation (Oudin vs. Penman-Monteith) does not substantially affect the long-term trend, although absolute values differ slightly. This supports the use of Oudin's formulation in KarstMod for assessing relative warming impacts. Moreover, as noted by Oudin et al., (2005), hydrological model parameter optimization against observed streamflow can partially compensate for PET uncertainties. Overall, while our results are primarily temperature-driven, they remain robust for discussing the impacts of warming on water resources in the Dévoluy catchment, acknowledging the residual uncertainties related to CO₂ effects, PET formulation, and precipitation input.

6 Conclusion

This study highlights the central role of snow in controlling recharge and discharge dynamics in the mid-mountain karst aquifer of the Dévoluy massif. Using a simple daily rainfall-snow-discharge model, we demonstrated that explicitly accounting for snow accumulation and melting is essential to reproduce both the timing and magnitude of karst spring discharge. The three-component structure of the model reveals contrasting response pathways, from rapid groundwater flow generated by rainfall events to delayed contributions from snowmelt and slow drainage sustaining groundwater baseflow. These results confirm that snow not only shifts recharge in time but also regulates its dynamics, thereby extending groundwater availability at karst spring during the summer period.

Warming scenarios have been tested for +2 and +4 °C temperature increases as advised by the French TRACC storyline (Soubeyroux et al., 2024), to help water managers incorporating climate change into environmental and water resources studies. Higher air temperature further shows that temperature rise leads to earlier and more variable snowmelt peaks, enhanced winter discharge, and reduced August discharge, with significant consequences for seasonal water availability, increasing water stress and potentially threatening water resources during critical low-flow periods. Mean August discharge has been modelled for the 2015–2019 period, and applied scenarios show that discharge is reduced by 26 and 42 %, respectively, under +2 °C and +4 °C warming, from 1.47±0.24 m³ s⁻¹ in Model 2 to 1.09±0.21 and 0.85±0.19 m³ s⁻¹.

These results provide valuable insights into the sensitivity of mountain karst aquifers to temperature forcing and claim toward adaptation strategies for water management under climate change.

Finally, this work illustrates both the strengths and limitations of rainfall-snow-discharge hydrological models in mountain karst environments. While they can capture key processes with limited data, uncertainties in precipitation inputs and water budget balance remain major challenges, especially in snow-dominated systems. Future research should combine improved precipitation datasets, including better accounting for snowfall, with climate projection scenarios to refine assessments of recharge and discharge dynamics. Such approaches are essential to anticipate the resilience and vulnerability of mid-mountain water resources in a warming climate.