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

Nathan Rispal, Bruno Arfib, Philippe Audra, Pierre Henry, Benoît Viguier, Alexandre Zappelli, Ludovic Mocochain, Marianna Jagercikova, Christine Vallet-Coulomb, Hélène Miche, et al.

► To cite this version:

Nathan Rispal, Bruno Arfib, Philippe Audra, Pierre Henry, Benoît Viguier, et al. Impact of seasonal snow on the recharge of a mountain karst aquifer under climate change: the Dévoluy case study (Southern Alps, France). 2026. hal-05465997

HAL Id: hal-05465997

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

Preprint submitted on 19 Jan 2026

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 CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

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

Nathan Rispal<fnref n="1" />, Bruno Arfib<fnref n="2" />, Philippe Audra<fnref n="1" />, Pierre Henry<fnref n="2" />, Benoît Viguier<fnref n="3" />, Alexandre Zappelli<fnref n="2" />, Ludovic Mocochain<fnref n="4" />, Marianna Jagercikova<fnref n="4" />, Christine Vallet-Coulomb<fnref n="2" />, Hélène Miche<fnref n="2" />, Laurent Cadilhac<fnref n="5" />

*Correspondence to*: Nathan Rispal (nathan.rispal@univ-cote-dazur.fr) and Bruno Arfib (arfib@cerege.fr).

**Abstract.** 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 melt in controlling both flow timing and magnitude. Under +2 °C 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 28% to 44%, 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.

## 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 m 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

<footnote n="1">¹ Univ. Côte d'Azur, Polytech'Lab- UPR 7498, Sophia- Antipolis, France</footnote>

<footnote n="2">²Aix Marseille Univ., CNRS, IRD, INRAE, CEREGE, Aix-en-Provence, France.</footnote>

<footnote n="3">³Géoazur, Univ. Côte d'Azur, CNRS, Observatoire de la Côte d'Azur, IRD, Valbonne, France.</footnote>

<footnote n="4">⁴Independent researchers.</footnote>

<footnote n="5">⁵Agence de l'Eau Rhône Méditerranée et Corse, Lyon, France.</footnote>
<figure><img src="image_3.png" /></figure>

numerous factors such as deposition thickness, weather conditions, slope exposure and soil infiltration capacity (Bayard et al., 2005; Eeckman et al., 2024; 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 & 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., 2014a). 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.

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., 2014a), 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 modeling 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., 2014a).

The impact of snow on the hydrodynamics of karst hydrosystems 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 modeling, Ç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 anticipate the effects of climate change on karst hydrological regimes (Doummar et al., 2018). These approaches make it possible to project the evolution of discharge rates under different warming scenarios using complex climatic projections. Finally Fan et al. (2023) have studied a mountainous karst system and the impact of snow on recharge, and show the relation between period of snowmelt and maximum spring discharge.
<figure><img src="image_4.png" /></figure>

Despite these advances, the understanding of hydrological functioning of mid-altitude karst systems influenced by snow remains limited. The strong spatio-temporal variability of snow accumulation and melt, 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 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 snow, yet it has never been investigated through hydrological modelling before. This limestone massif, composed of Mesozoic and Cenozoic formations folded and faulted by the Alpine orogeny, hosts karst developed mainly in Upper Cretaceous and Eocene limestones (Bonhomme, 1972; Luparini, 1975; Parizot et al., 2025). Altitudes range from 870 m asl. at the Gillardes Spring – the main outlet – to 2789 m asl. at the Obiou summit. Snowpack is generally present above 1700 m asl. from December to April.

The aim of this study is to (1) evaluate the ability of the KarstMod V3 platform (Sivelle et al., 2025) to simulate daily discharge of a karstic spring in a mid-altitude mountain environment, (2) quantify the contribution of seasonal snow to spring discharge, and (3) assess the impacts of temperature increase scenarios. The rainfall-snow-discharge model incorporates a degree-day snow routine applied to sub-catchments. It is implemented on the Dévoluy basin using four years of daily observations of precipitation, air temperature, and discharge. Sensitivity to snow dynamics is then explored by simulating modified air temperature scenarios. This approach highlights how temperature-driven changes in snow accumulation and melt influence the timing and magnitude of karst discharge providing first insights into the potential hydrological response of Alpine karst systems to future warming.

## 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 asl.), and displays the characteristic morphology of a perched syncline, with a central depression located around 900 m asl. surrounded by steep ridges over 2000 m asl. The basin opens out to the north through the Souloise gorges, the lowest altitude of the massif at 870 m asl., 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
<figure><img src="image_5.png" /></figure>

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 m and 800 m (Parizot et al., 2025).

These Senonian deposits are themselves unconformably overlain by Eocene conglomerates, Eocene-Oligocene nummulitic limestones and marls, 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 compartmentalize 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 asl.) 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 likely allows 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 asl., 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 asl.) 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 asl., the karstic Gillardes Spring is composed of two mains 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, and excludes the Bure Plateau drained southward.
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 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.</figcaption></figure>

## 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 x 8 km long, with 4 meshes covering the study area (# 7481, 7482, 7565, and 7566). The average value of the 4 selected meshes were used as input data for KarstMod model at daily time step from 2015 to 2019: precipitation (mm/day) and mean temperature (°C). The average altitude of these four meshes is 1642 m asl., close to the catchment area median altitude (1633 m asl.). Precipitation and

130

135
<figure><img src="image_8.png" /></figure>

temperature data from the in-situ meteorological Saint-Étienne station (belonging to Electricité de France - EDF), located in Saint-Étienne-en-Dévoluy at 1290 m asl., were contributed 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 snowfall are used as qualitative information on snow cover over the massif at the altitude of 1969 m asl.

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 a 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 four study years

The analyses in this section are presented by hydrological year (from September 1st to August 31st 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 climatological context, SAFRAN 30-year data (1990-2020) indicate an average annual precipitation of 1242 mm and 89 days per year 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 asl.

From 2015 to 2019, the proportion of snow extracted from SAFRAN data varied between 24% and 44% of total precipitation and the mean interannual value is 33% (391 mm). Most snowfall occurred between October and May, with a peak from December to March. 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.

Over the four years 2015 to 2019, snow amounts were highly variable, as illustrated by the in situ snow depth measurement at the ski resort (Fig. 2d) and the SAFRAN data (Fig. 2c and Table S1 in the supplementary material). 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. 2 b). The first, in autumn and early winter, corresponds to the largest rainfall events of the year (Fig. 2a.b). The second occurs in late spring, when snowmelt adds to the rainfall, maintaining a high discharge (e.g., spring 2018, Fig. 2 a.b). By contrast, low water occurs
<figure><img src="image_9.png" /></figure>

both in winter, when precipitation falls mainly as snow, and in summer when the snowpack is depleted or melted, and when precipitations are rare and makes little contribution to recharge.

During spring, discharge at Gillardes Spring increases and remains high up to June (Fig. 2b). 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, it accounts for 86% of the total discharge, while surface flow accounts for only 14%. From an annual water budget point of view, the previous percentage gives a rough estimate of the surface 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><img src="image_10.png" /></figure>

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

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

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

<figure><img src="image_14.png" /><figcaption>Fig. 2. (a) Daily precipitation and air temperature from SAFRAN reanalysis; (b) Observed discharge at the Gillardes Spring and Souloise river at Pont de la Baume station; (c) annual snow quantity from SAFRAN data, (d) measured in situ snow data at Super Dévoluy ski resort (1976 m asl): snow depth (daily height of snowpack), and daily observed snowfall. In (d), “No data” shaded areas correspond to times during which the in situ station is closed.</figcaption></figure>
<figure><img src="image_15.png" /></figure>

# 3. Method : Rainfall-snow-discharge modelling

## 3.1. KarstMod

KarstMod is a rainfall-discharge modelling platform dedicated to karst systems (Mazzilli et al., 2019, 2023, Sivelle et al., 2025), developed by the National KARST Observation Service (SNO KARST, France). KarstMod V.3.0.13 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), model 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 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 the Oudin formulation (Oudin, 2006). 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, 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:

$$Q_{AB} = k_{AB} \left(\frac{A}{L_{ref}}\right)^{\alpha_{AB}} \quad \text{if } A>0$$

$$Q_{AB} = 0 \quad \text{otherwise}$$

where A is the water level in compartment A [L] and $k_{AB}$ is the specific discharge coefficient [L/T] for the discharge law from compartment A to compartment B, $L_{ref}$ [L] is a unit length and $\alpha_{AB}$ [-] is a positive exponent. $L_{ref}$ has been introduced in KarstMod to simplify the units when the law is not linear (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 day⁻¹. Assigning distinct coefficients thus enable 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 Emin in E, which simulates a stored quantity of water only available for evapotranspiration, contained between zero and Emin (Emin<0). For modelling consistency, each compartment is also
<figure><img src="image_16.png" /></figure>

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. Modeling is divided into three stages, guided by the KarstMod GUI and corresponding to successive time interval in the dataset : 1) Warm-up, to avoid any model initialization bias. 2) Calibration, during which parameters are optimized to obtain the best match between simulated discharge 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 (Wobj) calculated for each parameter set. The user defines a minimum value of the objective function (Wobj_min value) to be satisfied in order to retain the parameter set. The set with the highest Wobj value is then used for validation, which is based on the same objective function, applied to a different time interval. Finally, the platform integrates a set of graphical visualizations to facilitate evaluation of model performance and interpretation of results.

**3.2. Model optimization**

Model performance is assessed by comparing simulated discharge with observed discharge. For this purpose, KarstMod offers a set of tools based on statistical criteria and visual exports in the form of graphs.

Several performance criteria are available in KarstMod (Mazzilli et al., 2023) to calculate the Wobj 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 good evaluation of the fit between observed discharge and simulated discharge in a wide variety of contexts. It is a multi-component criterion based on the diagonal decomposition of the Nash-Sutcliffe Efficiency (NSE) to separate correlation, bias, and variability that is comparable to multi-objective criteria for calibration purpose.

KGEnp can vary between -∞ and 1. A value of 1 indicates a perfect match between observed and simulated discharge. In this study, we set a minimum performance of 0.8.

KarstMod offers other tools to assess the validity of the results, which are freely interpreted by the modeler: 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 “expert” validation by the modeler is part of the modeling procedure. In the specific case of snow-influenced mid-altitude karstic hydro systems, beyond the performance criterion for the entire simulated series (KGEnp), the modeler will seek to validate the consistent consideration of snowpack evolution: (1) cold winter period generate snow stock, and Psr (meltwater) becomes lower than P; (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><img src="image_17.png" /></figure>

### 3.3. Snow routine

Version 3 of 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 (Psr) 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).

The KarstMod main model is not a distributed model, but the snow routine can be applied to sub-catchments dividing the catchment area into several altitudinal units (sub-catchments), 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., 2014b), 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.

<figure><img src="image_18.png" /></figure>
<figure><img src="image_19.png" /></figure>

Fig. 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 (Psr). E, M and C are compartments for flow transfer to the Spring (S). QES is a direct flow from E to the spring.

**3.4. Modelling strategy**

This section presents the step-by-step procedure used to model rainfall-snow-discharge processes at the Gillardes karst Spring using the KarstMod model over the period 2015 to 2019. The modelling approach is conducted using a limited number of calibration parameters and input time series. As a conceptual model, it is based on the modeler’s understanding of key atmosphere-surface-groundwater interactions, with a focus on two main processes: recharge together with fast and slow groundwater flows.

The Dévoluy massif, characterized by seasonal snow, is a relevant case study: the Gillardes Spring show typical karstic behavior, with rapid responses during high-flow events and sustained baseflow during dry periods in winter and summer.

The modeling 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 performed using the previously calibrated Model 2, including a scenario with the snow routine deactivated (Model 2 SRoff) and 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 (September 1, 2018, to August 31, 2019), characterized by intermediate snow conditions. Initially, we set Wobj min=0.5 for the first steps of modeling, and then increased it to Wobj min=0.8, retaining the 10,000 simulations with Wobj>Wobj min to ensure model robustness. Among these 10,000 simulations, the one that maximizes Wobj was selected and used to generate graphical outputs. Cross calibration-validation tests did not improve performance, with similar KGENp values remaining above 0.8.

**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. Emin 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
<figure><img src="image_20.png" /></figure>

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², the parameter range is set from 150 to 250 km².

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 Emin, 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. Snow routine with sub-catchments is used because the simplification of a uniform distribution of snow or rain across the Dévoluy massif is too rough in a first approximation. Indeed, actual elevations range from 870 m asl. at the Spring to 2789 m asl. on the summits. Every rainy day, therefore, there are areas of the catchment that receive rain while at higher altitudes it snows, and the evolution of the snowpack over time differs between high-, medium- and low-altitude zones. Modeling can therefore be further refined by dividing into sub-catchments by elevation band to better discretize snow storage and melting zones.

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 we tested from 1 to 5 elevation bands and did not notice any improvement between 4 and 5 bands. So, for the sake of parsimony, four elevation bands were set up, each representing 25% of the total catchment area (Fig. 5). The boundaries of these classes are supported by observations of the extent of winter snow cover in the Dévoluy massif.

The SB1 sub-catchment (870-1414 m asl.) corresponds to the least snow-covered zone, with rare snowfalls and a short-lived snow cover. SB2 (1414-1634 m asl.) covers the lower altitudes of the ski resort, with regular snowfalls but little persistent snow cover. SB3 (1634-1980 m asl.) marks the start of the quasi-permanent winter snowpack. Finally, SB4 (1980-2789 m asl.) 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 negative as the altitude rises, and positive as it falls. This shift is calculated by multiplying the difference in altitude between the SAFRAN measurement point (1642 m asl.) and the median altitude of the sub-catchment, by a regional gradient of 0.6°C/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. 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.°C⁻¹.day⁻¹ for MF, 0.3 to 0.8 for CFR and 0.05 to 0.25 for CWH.
<figure><img src="image_21.png" /></figure>

However, one parameter caught our attention, the temperature threshold $T_s$. It has been defined as follows: when the air temperature is below $T_s$, precipitation accumulates as snow, and when the temperature exceeds $T_s$, melting occurs according to a melting coefficient MF. We highlight that $T_s$ does not have to be zero. Its value is calibrated and depends on the altitude at which the air temperature of the input file is measured, and on the median altitude of the catchment area. For example, if the air temperature is measured at low altitude, then temperatures will be higher than the representative catchment average; in this case, $T_s$ will have a positive value (up to several degrees Celsius) to enable the model to accumulate snow at catchment scale. Conversely, $T_s$ will be negative if the air temperature input is measured at high altitude, above the median altitude.
We have also tested other structures (number of compartments and flows) that did not bring any significant improvement. We mention them here, but we will not show the results. Such structures can also be found in other case studies in the literature (i.e. Hartmann et al., 2012; Sivelle et al., 2025), such as adding a flow through an additional parallel compartment (L), or adding a hysteresis threshold, or an exchange between the M and C flows.

<figure><img src="image_22.png" /><figcaption>Fig. 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.</figcaption></figure>

### 3.4.3. Scenario simulations : altered temperature and snow regimes

We explore the sensitivity of simulated discharge to snow dynamics by modifying air temperature, which affects rain/snow distribution and snowmelt. The previously calibrated Model 2 is used to test these changes, without further calibration or structural modifications. The model is run in “run mode” with fixed parameters, and total precipitation input remains identical to the original simulation.

First, we consider a scenario in which the snow routine is deactivated (Model 2 SOff) using previously calibrated Model 2 with SAFRAN input precipitation and temperature. In this configuration, all precipitation is treated as rainfall, regardless of temperature. The goal is to simulate a fully rain-dominated regime, in which no snow is accumulated or melted throughout the year. This approach allows isolating the role of the snow module and assessing its specific contribution to discharge generation in a karst system influenced by seasonal snow dynamics.
<figure><img src="image_23.png" /></figure>

Then we examine the effect of a reduction in snowfall and modification of snowmelt on simulated discharges caused by air temperature increase (Fig. 4). Historical data and global and regional climate models show an already warming in France in 2020 between +1.5°C 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, reference warming trajectory for climate change adaptation have been proposed in 2024 as storyline named TRACC (Soubeyroux et al., 2024), to help water managers to incorporate climate change in environmental or water resources studies. We thus use the French “TRACC” recommendation to fix the temperature increase to model. We tested two scenarios of temperature change with the KarstMod 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).

Importantly, in Model 2 SRoff, no changes are made to the input temperature or to the potential evapotranspiration (PET), which remains identical to the reference simulation (Model 2). This distinguishes Model 2 SRoff from the temperature-shifted models (T+2 and T+4); these latter models incorporate a PET calculation influenced by warmer temperatures. In Model 2 SRoff, any change in simulated discharge arises solely from the removal of snow-related processes and the immediate routing of precipitation into the hydrological system.

While KarstMod is not specifically designed for climate change projections, it works like any other lumped hydrological models (e.g., GR, MORDOR), which have been successfully used in climate sensitivity studies such as the Explore2 project (Sauquet et al., 2025). These simulations allow assessment of the potential consequences of warming on karst aquifer recharge and discharge dynamics under simplified but plausible future conditions.
<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /><figcaption>Fig. 4. Schematic representation of the modeling 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. Three scenarios are thensimulated 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°C and +4°C, respectively), Model 2 SRoff with snow routine deactivated.</figcaption></figure>

380
<figure><img src="image_26.png" /></figure>

<table><caption>Table 1. Calibrated parameters and objective function values for Model 1 and Model 2.</caption><thead><tr><th>Parameter</th><th></th><th>Model 1</th><th>Model 2</th></tr></thead><tbody><tr><td>k<sub>ES</sub> [day<sup>-1</sup>]</td><td>Recession coefficient from compartment E to the Spring</td><td>2.00 x10<sup>-1</sup></td><td>2.42x10<sup>-1</sup></td></tr><tr><td>k<sub>EM</sub> [day<sup>-1</sup>]</td><td>Recession coefficient from compartment E to M</td><td>1.42x10<sup>-1</sup></td><td>6.22x10<sup>-2</sup></td></tr><tr><td>k<sub>EC</sub> [day<sup>-1</sup>]</td><td>Recession coefficient from compartment E to C</td><td>8.46x10<sup>-2</sup></td><td>2.10x10<sup>-1</sup></td></tr><tr><td>k<sub>MS</sub> [day<sup>-1</sup>]</td><td>Recession coefficient from compartment M to the Spring</td><td>8.91x10<sup>-3</sup></td><td>9.35x10<sup>-3</sup></td></tr><tr><td>k<sub>CS</sub> [day<sup>-1</sup>]</td><td>Recession coefficient from compartment C to the Spring</td><td>7.95x10<sup>-2</sup></td><td>6.45x10<sup>-2</sup></td></tr><tr><td>E<sub>0</sub> [mm]</td><td>Initial water level in E</td><td>5.84</td><td>37.33</td></tr><tr><td>M<sub>0</sub> [mm]</td><td>Initial water level in M</td><td>32.59</td><td>23.30</td></tr><tr><td>C<sub>0</sub> [mm]</td><td>Initial water level in C</td><td>44.88</td><td>6.26</td></tr><tr><td>E<sub>min</sub> [mm]</td><td>Minimum water level in E</td><td>-60.20</td><td>-32.26</td></tr><tr><td>RA [km<sup>2</sup>]</td><td>Recharge area</td><td>207.95</td><td>191.00</td></tr><tr><td>Latitude [°]</td><td>Latitude of the area for PET routine</td><td>44.6</td><td>44.6</td></tr><tr><td>K1 [°C]</td><td>Fixed parameter for PET routine</td><td>100</td><td>100</td></tr><tr><td>K2 [°C]</td><td>Fixed parameter for PET routine</td><td>5</td><td>5</td></tr><tr><td>T<sub>S</sub> [°C]</td><td>Snowmelt temperature threshold</td><td>/</td><td>1.74</td></tr><tr><td>MF [mm. °C<sup>-1</sup> .d<sup>-1</sup>]</td><td>Melt factor</td><td>/</td><td>4.09</td></tr><tr><td>CFR [ - ]</td><td>Refreezing coefficient</td><td>/</td><td>4.44x10<sup>-1</sup></td></tr><tr><td>CWH [ - ]</td><td>Water capacity holding of snow</td><td>/</td><td>5.49x10<sup>-2</sup></td></tr><tr><td>KGEnp (calibration)</td><td>Wobj value for calibration</td><td>0.65</td><td>0.89</td></tr><tr><td>KGEnp (validation)</td><td>Wobj value for validation</td><td>0.56</td><td>0.80</td></tr></tbody></table>

# 4. Results

## 4.1. Simulation of the Gillardes Spring discharge

### 4.1.1. Model 1: without snow routine

Calibrated parameters for Model 1 are given in Table 1. 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 (
<figure><img src="image_27.png" /></figure>

Fig. 4a). 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. 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 KGEnp value of 0.65 over the calibration period and 0.56 in validation, indicating only moderate ability to reproduce both the timing and magnitude of flows (Table 1).

### 4.1.2. Model 2: snow routine activated

Calibrated parameters for Model 2 are given in Table 1. Here the snow routine is activated and the simulation improves substantially (

Fig. 4b). Precipitations 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 four sub-catchments 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 value of 0.89 over the calibration period and 0.80 in validation (Table 1). Calculated KGEnp year by year (Table S2 in supplementary material) ranges between 0.80 and 0.90 (0.83 for the mean value over the 4 studied years). This highlights the model's ability to reproduce discharge in different hydrological conditions, from low to high influence of seasonal snow over the 4 years of modelling. Actual evapotranspiration (AET) calculated with KarstMod, as a water budget in compartment E, is in average 36% of total precipitation.

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 (Psr) 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 Emin 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
<figure><img src="image_28.png" /></figure>

recession coefficient) accounts for discharge during spring season high-water flow influenced by snowmelt and during peak events from rainfall or snowmelt.

425

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 more robust performance criteria, which strengthens its reliability for further scenario analysis.

<figure><img src="image_29.png" /><figcaption>Fig. 4. Simulated discharge of the Gillardes Spring compared to the observed discharge, between May 2015 to September 2019, (a) with Model 1 (without snow routine), (b) Model 2 (with snow routine) and (c) snow comparison between in situ observed snow depth at the Super-Dévoluy ski resort versus simulated water stored in snowpack in the sub-catchments of the snow routine.</figcaption></figure>

430
<figure><img src="image_30.png" /></figure>

<figure><img src="image_31.png" /><figcaption>Fig. 5. (a) Model 2 seven-day averaged simulated water levels of KarstMod compartments E, M, C. (b) Model 2 seven-day averaged internal KarstMod discharges (QES: from E to S, QMS: from M to S, QCS: from C to S). Seven-day precipitation is shown on a secondary y-axis (inverted).</figcaption></figure>

## 4.2. Scenario simulations

### 4.2.1. All precipitation forced as rainfall

What would have been the discharge rate of the Gillardes Spring if the precipitation from 2015 to 2019 had fallen only as rain? To address this question, we simulated a scenario (Model 2 SRoff) in which the catchment no longer received snowfall by deactivating the snow routine in KarstMod, without temperature change. This model differs from Model 1 in that it uses the calibrated Model 2 that gave validated results, whereas Model 1 was not validated.

Annual mean simulated discharges for the Model 2 SRoff range from 3.61 to 5.43 m³·s⁻¹, equivalent to Model 2 ranging between 3.59–5.38 m³·s⁻¹ (Table S2). Simulated discharge with Model 2 and Model 2 SRoff is plotted on seven-day average in Fig. 8b. The main differences appear in two periods: (1) during late autumn and winter, (2) during spring and early summer.

In autumn and winter, simulated discharge without the snow routine is higher: total precipitation falls as rain and immediately contributes to recharge into the compartments. At Gillardes Spring, both base flow and flood peaks increase. Conversely, in spring and early summer, the flow rate without snow routine does not show the spring high-water period. The discrepancy is all the more noticeable in spring 2018 and 2019, as these were the two snowiest years (Table S1). The difference between the two curves illustrates the drop in flow in spring if the snowmelt no longer supports the flow (Model 2 SRoff). During the low
<figure><img src="image_32.png" /></figure>

water period at the end of summer and the high discharge in early autumn, the two curves overlap, as there is no remaining
snow in this mid-mountain catchment.

4.2.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 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. 8c) 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 modeled years, the mean simulated discharge decreases from 4.58 m³·s⁻¹ for Model 2 to 4.34 m³·s⁻¹ for Model 2 T+2, corresponding to a reduction of approximately 5.2% 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. 8d): 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 modeled years, the mean simulated discharge decreases from 4.58 m³·s⁻¹ for Model 2 to 4.09 m³·s⁻¹ for Model 2 T+4, corresponding to a reduction of approximately 10.7% in simulated discharge (Table S2).

Finally, the comparison between simulated discharge from Model 2 T+4 and Model 2 SRoff shows very few differences (Fig. 8e). Both models produce nearly identical discharge dynamics, with similar timing and magnitude of spring peaks as well as comparable baseflow levels during summer. The main difference is a lower discharge for Model 2 T+4 because the rise of temperature increases the calculated PET by Oudin formula and the AET calculated by water balance in compartment E. 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, 39% and 43%. 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><img src="image_33.png" /></figure>

<figure><img src="image_34.png" /><figcaption>Fig. 6. Seven-day average of simulated discharges ($m^3.s^{-1}$) under different temperature scenarios. (a) snow quantity for each hydrological year; (b) Model 2 SR Off and the reference Model 2; (c) Model 2 T+2 and Model 2; (d) Model 2 T+4 and Model 2; (e) Model 2 T+4 and Model 2 SR Off.</figcaption></figure>

480
<figure><img src="image_35.png" /></figure>

# 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 modeling 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. This kind of routine, widely used in hydrological and hydrogeological models (Bittner et al., 2018; Seibert and Vis, 2012; Valéry et al., 2014a) 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 (QES), 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.

The comparison of Model 2 and Model 2 SRoff during snowy year (Fig. 8b) shows that current high-water discharge during spring season is highly supplied by delayed snowmelt recharge (-26% of discharge during spring when snow routine is deactivated). 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.

Comparison with a rainfall-only version of the model quantitatively highlights the effect of snow on both the timing and the magnitude of Spring discharge. This difference is particularly relevant for water managers, as snow storage and gradual melting extend the availability of water in summer, supporting human uses such as agriculture, tourism and drinking water supply. 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 (9.35x10⁻³ day⁻¹). This parameter provides a quantitative metric of the aquifer’s functioning and allows the simulation of groundwater discharge depletion during snow-free summer periods.

## 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 melt, which in turn affect the timing and magnitude of recharge. By modifying 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
<figure><img src="image_36.png" /></figure>

uncertainties linked to precipitation (Duethmann et al., 2020). These scenarios are not intended as predictive forecasts but rather as a way to explore how the system responds to temperature forcing.

Seasonal discharge distribution is given in Fig. 7. In winter, mean flows increase in the warming scenarios (T+2°C and T+4°C), with mean discharge rising from 3.11 m³·s⁻¹ in the baseline to 4.66 m³·s⁻¹ (+50%) and 5.70 m³·s⁻¹ (+83%), respectively, suggesting 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, with a mean of 6.42 m³·s⁻¹ in the T+2 °C scenario, compared with 7.24 m³·s⁻¹ in the reference scenario (-11%), and further decreasing to 5.05 m³·s⁻¹ in the T+4 °C scenario (-30%). This reduction is consistent with reduced snow accumulation during winter, resulting in lower spring discharge. This pattern, observed in the seasonal seven-day 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. During the summer season, mean discharge decreases as the temperature scenarios increase: the mean discharge drops from 3.54 m³·s⁻¹ in the reference scenario (Model 2) to 1.86 m³·s⁻¹ in the T+2 °C scenario (-48%), and further to 1.42 m³·s⁻¹ (-60%) in the T+4 °C scenario. 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, monthly mean discharge (over the 4 years studied) decreases from 1.19 m³·s⁻¹ in Model 2 (and 1.30 m³·s⁻¹ in observed data) to 0.86 m³·s⁻¹ in Model 2 T+2 (-28% vs Model 2) and 0.67 m³·s¹ (-44% 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. But it is important to note that karst discharge decrease will be even greater in a year when spring precipitation are 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). While median simulated discharges for the three models shown in Fig. 9 are quite similar, 7-day averaged hydrographs (Fig. 8) and monthly average (Fig. 10) reveal different behaviors than seasonal statistics: 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.12 m³·s⁻¹) than Model 2 T+2 and Model 2 T+4 (0.85 m³·s⁻¹ and 0.69 m³·s⁻¹), whereas from November onward, especially during the first snow events, flows with +2°C and +4°C become higher (monthly average jumping from 7.02 m³·s⁻¹ in Model 2 to 8.20 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 boxplots, 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 reduced by half. Moreover, temperature increase impacts evapotranspiration (ET). Potential ET is currently correlated with temperature in Oudin formula, generating PET increase for
<figure><img src="image_37.png" /></figure>

scenarios T+2°C and T+4°C. Actual ET calculated by KarstMod water budget in compartment E also increases in average from 417 mm/year to 457 mm/year (+9.5%) and 498 mm/year (+19.4%), respectively for Model 2, Model 2 T+2, Model 2 T+4, indicating increased water stress and reduced water availability for discharge (Table S3). 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. A fully prospective study using climate projection data could be developed in future work to extend this analysis (Marson et al., 2024).

<figure><img src="image_38.png" /><figcaption>Fig. 7. Seasonal distribution of simulated discharge (7-day mean) for winter (Dec, Jan., Feb.), spring (Mar., Apr. May), summer (Jun., Jul., Aug.), and autumn (Sep., Oct., Nov.). Boxplots summarize the variability of each scenario within each season (the whiskers are drawn within the 1.5 interquartile range, the box gives 25th percentile, 50th percentile and 75th percentile). Black circles give the mean.</figcaption></figure>

<figure><img src="image_39.png" /><figcaption>Fig. 8. Mean monthly discharge at Gillardes Spring over the 4 hydrological years with observed discharge and modelled discharge from Model 2 and Model 2 T+2 and Model 2 T+4.</figcaption></figure>
<figure><img src="image_40.png" /></figure>

560 **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. KarstMod has no water budget adjustment parameter, but the “RA: catchment recharge area” parameter acts as an adjustment parameter.

In a study case where PET and precipitation are well constrained, this 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 catchment area (RA) equal to 191 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

565

570

575

580

585

590
<figure><img src="image_41.png" /></figure>

are non-existent. The RA calibrated value is almost 1.3 times higher than the expected value, which would correspond to an increase of 27% in total precipitation over the Dévoluy massif.

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 modeling 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. Thus, future ET evolution may be more complex than a simple temperature–PET relationship. In a complementary study, Lemaitre-Basset (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 the Oudin formulation in KarstMod for assessing relative warming impacts. Moreover, as noted by Oudin (2006), 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.

## 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°C 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
<figure><img src="image_42.png" /></figure>

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 modeled for the 2015–2019 period, and applied scenarios show that discharge is reduced by 28 and 44%, respectively, under +2°C and +4°C warming, from 1.19 m³·s⁻¹ in Model 2 to 0.86 m³·s⁻¹ and 0.67 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.

Data availability

The KarstMod modelling software used in this study is freely available online (https://hal.science/hal-02071006). The discharge data from the Gillardes Spring for the study period are available at: https://data.oreme.org/snokarst/snokarst_map#zoom=9&lat=44.954&lon=5.934&layer=OpenStreetMap&overlays=TTT&selected Layer=snokarst station&selected station=44.76032.

SAFRAN reanalysed meteorological data are accessible via the GeoSAS platform (https://geosas.fr/edr-viewer/), and snow measurements (ID 07595 station Super Dévoluy) are available from Météo-France public data (https://donneespubliques.meteofrance.fr/?fond=produit&id produit=94&id rubrique=32).

Author contribution

Nathan Rispal: Conceptualization, Investigation, Methodology, Writing - original draft. Bruno Arfib: Conceptualization, Investigation, Methodology, Writing - original draft. Philippe Audra: Funding acquisition, Supervision, Project administration, Writing - original draft. Pierre Henry: Writing - original draft. Benoît Viguier: Writing - review and editing. Alexandre Zappelli: Writing - review and editing, Ludovic Mocochain: Writing - review and editing. Marianna Jagercikova: Writing - review and editing. Christine Vallet-Coulomb: Writing - review and editing. Hélène Miche: Writing - review and editing. Laurent Cadilhac: review and editing.
<figure><img src="image_43.png" /></figure>

Competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared
to influence the work reported in this paper.

Acknowledgement

We thank EDF for providing stream flow data used in this study. Financial support was provided by the Agence de l'Eau Rhône Méditerranée Corse (AERMC project n°2023 5381) and Comité de spéléologie Régional Sud (CSR Sud). We are also grateful to Guillaume Cinkus, Vianney Sivelle, and Axel Flinck for valuable discussions that improved this work. This work was carried out within the framework of the French Karst observatory network (SNO Karst) initiative of the INSU/CNRS, which aims to strengthen knowledge-sharing and promote cross-disciplinary research on karst systems.

References

Andrew, J. T. and Sauquet, E.: Climate Change Impacts and Water Management Adaptation in Two Mediterranean-Climate Watersheds: Learning from the Durance and Sacramento Rivers, Water, 9, 126, https://doi.org/10.3390/w9020126, 2017.

Données Publiques de Météo-France - Observations du réseau nivo-météorologique: https://donneespubliques.meteofrance.fr/?fond=produit&id_produit=94&id_rubrique=32, last access: 16 October 2025.

Audra, P. and Palmer, A.: The Vertical Dimension of Karst: Controls of Vertical Cave Pattern, Treatise Geomorphol., 6, 186–206, https://doi.org/10.1016/B978-0-12-374739-6.00098-1, 2013.

Barnett, T., Adam, J., and Lettenmaier, D.: Potential Impacts of a Warming Climate on Water Availability in Snow-Dominated Regions, Nature, 438, 303–9, https://doi.org/10.1038/nature04141, 2005.

Baudement, C., Arfib, B., Mazzilli, N., Jouves, J., Lamarque, T., and Guglielmi, Y.: Groundwater management of a highly dynamic karst by assessing baseflow and quickflow with a rainfall-discharge model (Dardennes springs, SE France), Bull. Société Géologique Fr., 188, 40, https://doi.org/10.1051/bsgf/2017203, 2017.

Bayard, D., Stähli, M., Parriaux, A., and Flühler, H.: The influence of seasonally frozen soil on the snowmelt runoff at two Alpine sites in southern Switzerland, J. Hydrol., 309, 66–84, https://doi.org/10.1016/j.jhydrol.2004.11.012, 2005.

Bergström, S.: Development and Application of a Conceptual Runoff Model for Scandinavian Catchments, 1976.

Bergström, S.: The HBV model - its structure and applications., SMHI RH, 1–44, 1992.

Bittner, D., Narany, T. S., Kohl, B., Disse, M., and Chiogna, G.: 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, 2018.

Bonhomme, J. L.: Etude hydrogéologique et hydrodynamique du karst des calcaires sénoniens. Massif du Dévoluy, Hautes-Alpes, Conservatoire national des arts et métiers CNAM, France, 312 pp., 1972.
<figure><img src="image_44.png" /></figure>

680 Branche, E.: The Durance-Verdon River Basin in France: The Role of Infrastructures and Governance for Adaptation to Climate Change, in: Increasing Resilience to Climate Variability and Change: The Roles of Infrastructure and Governance in the Context of Adaptation, edited by: Tortajada, C., Springer, Singapore, 129–155, https://doi.org/10.1007/978-981-10-1914-2_7, 2016.

685 Çallı, S. S., Çallı, K. Ö., Tuğrul Yılmaz, M., and Çelik, M.: Contribution of the satellite-data driven snow routine to a karst hydrological model, J. Hydrol., 607, 127511, https://doi.org/10.1016/j.jhydrol.2022.127511, 2022.

690 Calvin, K., Dasgupta, D., Krinner, G., Mukherji, A., Thorne, P. W., Trisos, C., Romero, J., Aldunce, P., Barrett, K., Blanco, G., Cheung, W. W. L., Connors, S., Denton, F., Diongue-Niang, A., Dodman, D., Garschagen, M., Geden, O., Hayward, B., Jones, C., Jotzo, F., Krug, T., Lasco, R., Lee, Y.-Y., Masson-Delmotte, V., Meinshausen, M., Mintenbeck, K., Mokssit, A., Otto, F. E. L., Pathak, M., Pirani, A., Poloczanska, E., Pörtner, H.-O., Revi, A., Roberts, D. C., Roy, J., Ruane, A. C., Skea, J., Shukla, P. R., Slade, R., Slangen, A., Sokona, Y., Sörensson, A. A., Tignor, M., Van Vuuren, D., Wei, Y.-M., Winkler, H., Zhai, P., Zommers, Z., Hourcade, J.-C., Johnson, F. X., Pachauri, S., Simpson, N. P., Singh, C., Thomas, A., Totin, E., Alegría, A., Armour, K., Bednar-Friedl, B., Blok, K., Cissé, G., Dentener, F., Eriksen, S., Fischer, E., Garner, G., Guivarch, C., Haasnoot, M., Hansen, G., Hauser, M., Hawkins, E., Hermans, T., Kopp, R., Leprince-Ringuet, N., Lewis, J., Ley, D., Ludden, C., Niamir, L., Nicholls, Z., Some, S., Szopa, S., Trewin, B., Van Der Wijst, K.-I., Winter, G., Witting, M., Birt, A., Ha, M., Arias, P., Bustamante, M., Elgizouli, I., Flato, G., Howden, M., Méndez-Vallejo, C., Pereira, J. J., Pichs-Madruga, R., Rose, S. K., Saheb, Y., Sánchez Rodríguez, R., Ürge-Vorsatz, D., Xiao, C., Yassaa, N., et al.: IPCC, 2023: Climate Change 2023: Synthesis Report. Contribution of Working Groups I, II and III to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change [Core Writing Team, H. Lee and J. Romero (eds.)]. IPCC, Geneva, Switzerland., Intergovernmental Panel on Climate Change, https://doi.org/10.59327/ipcc/ar6-9789291691647, 2023.

700 Chen, Z., Auler, A. S., Bakalowicz, M., Drew, D., Griger, F., Hartmann, J., Jiang, G., Moosdorf, N., Richts, A., Stevanovic, Z., Veni, G., and Goldscheider, N.: The World Karst Aquifer Mapping project: concept, mapping procedure and map of Europe, Hydrogeol. J., 25, 771–785, https://doi.org/10.1007/s10040-016-1519-3, 2017.

Christensen, J. H., Carter, T. R., Rummukainen, M., and Amanatidis, G.: Evaluating the performance and utility of regional climate models: the PRUDENCE project, Clim. Change, 81, 1–6, https://doi.org/10.1007/s10584-006-9211-6, 2007.

705 Cinkus, G., Mazzilli, N., Jourde, H., Wunsch, A., Liesch, T., Ravbar, N., Chen, Z., and Goldscheider, N.: When best is the enemy of good – critical evaluation of performance criteria in hydrological models, Hydrol. Earth Syst. Sci., 27, 2397–2411, https://doi.org/10.5194/hess-27-2397-2023, 2023.

De Waele, J. and Gutierrez, F.: Karst Hydrogeology, Geomorphology and Caves, https://doi.org/10.1002/9781119605379, 2022.

710 Dettinger, M., Redmond, K., and Cayan, D.: Winter orographic precipitation ratios in the Sierra Nevada - Large-scale atmospheric circulations and hydrologic consequences, J. Hydrometeorol., 5, 1102–1116, https://doi.org/10.1175/JHM-390.1, 2004.

715 Doummar, J., Hassan Kassem, A., and Gurdak, J. J.: Impact of historic and future climate on spring recharge and discharge based on an integrated numerical modelling approach: Application on a snow-governed semi-arid karst catchment area, J. Hydrol., 565, 636–649, https://doi.org/10.1016/j.jhydrol.2018.08.062, 2018.

Duethmann, D., Blöschl, G., and Parajka, J.: Why does a conceptual hydrological model fail to correctly predict discharge changes in response to climate change?, Hydrol. Earth Syst. Sci., 24, 3493–3511, https://doi.org/10.5194/hess-24-3493-2020, 2020.
<figure><img src="image_45.png" /></figure>

Eckert, N., Corona, C., Giacona, F., Gaume, J., Mayer, S., van Herwijnen, A., and Hagenmuller, P.: Climate change impacts on snow avalanche activity and related risks, Nat. Rev. Earth Environ., 5, https://doi.org/10.1038/s43017-024-00540-2, 2024.

Eeckman, J., De Grenus, B., Miesen, F., Thornton, J., Brunner, P., and Peleg, N.: Dynamics of snow melt infiltration into mountain soils: an instrumental approach in the Nant Valley, Swiss Alps, EGUSphere, 1–23, https://doi.org/10.5194/egusphere-2024-1832, 2024.

Evin, G., Somot, S., and Hingray, B.: Balanced estimate and uncertainty assessment of European climate change using the large EURO-CORDEX regional climate model ensemble, Earth Syst. Dyn., 12, 1543–1569, https://doi.org/10.5194/esd-12-1543-2021, 2021.

Evin, G., Le Lay, M., Fouchier, C., Penot, D., Colleoni, F., Mas, A., Garambois, P.-A., and Laurantin, O.: Evaluation of hydrological models on small mountainous catchments: impact of the meteorological forcings, Hydrol. Earth Syst. Sci., 28, 261–281, https://doi.org/10.5194/hess-28-261-2024, 2024.

Fan, X., Goeppert, N., and Goldscheider, N.: Quantifying the historic and future response of karst spring discharge to climate variability and change at a snow-influenced temperate catchment in central Europe, Hydrogeol. J., 31, 2213–2229, https://doi.org/10.1007/s10040-023-02703-9, 2023.

Fayad, A., Gascoin, S., Faour, G., López-Moreno, J. I., Drapeau, L., Page, M. L., and Escadafal, R.: Snow hydrology in Mediterranean mountain regions: A review, J. Hydrol., 551, 374–396, https://doi.org/10.1016/j.jhydrol.2017.05.063, 2017.

Ferry, S. and Grosheny, D.: Excursion dans le sud du bassin subalpin. Turbidites et brèches carbonatées (Jurassique et Crétacé du bassin subalpin). Excursion du Groupe Français du Crétacé., Groupe Français du Crétacé, France, 2013.

Ford, D. and Williams, P. D.: Karst Hydrogeology and Geomorphology, John Wiley & Sons, 581 pp., 2007.

Frank, S., Goeppert, N., and Goldscheider, N.: Improved understanding of dynamic water and mass budgets of high-alpine karst systems obtained from studying a well-defined catchment area, Hydrol. Process., 35, e14033, https://doi.org/10.1002/hyp.14033, 2021.

Franz, K. J., Hogue, T. S., and Sorooshian, S.: Operational snow modeling: Addressing the challenges of an energy balance model for National Weather Service forecasts, J. Hydrol., 360, 48–66, https://doi.org/10.1016/j.jhydrol.2008.07.013, 2008.

Garavaglia, F., Le Lay, M., Gottardi, F., Garçon, R., Gailhard, J., Paquet, E., and Mathevet, T.: Impact of model structure on flow simulation and hydrological realism: from a lumped to a semi-distributed approach, Hydrol. Earth Syst. Sci., 21, 3937–3952, https://doi.org/10.5194/hess-21-3937-2017, 2017.

Garçon, R.: Prévision opérationnelle des apports de la Durance à Serre-Ponçon à l'aide du modèle MORDOR. Bilan de l'année 1994-1995, Houille Blanche, 71–76, https://doi.org/10.1051/lhb/1996056, 1996.

Gidon, M. and Pairis, J.-L.: Le rôle des mouvements tectoniques éocènes dans la genèse des structures de retrémité NE du Dévoluy et dans celle du chevauchement de Digne, Géologie Alp., 73–83, 1976.

Goldscheider, N., Chen, Z., Auler, A. S., Bakalowicz, M., Broda, S., Drew, D., Hartmann, J., Jiang, G., Moosdorf, N., Stevanovic, Z., and Veni, G.: Global distribution of carbonate rocks and karst water resources, Hydrogeol. J., 28, 1661–1677, https://doi.org/10.1007/s10040-020-02139-5, 2020.
<figure><img src="image_46.png" /></figure>

Gottardi, F., Obled, C., Gailhard, J., and Paquet, E.: Statistical reanalysis of precipitation fields based on ground network data and weather patterns: Application over French mountains, J. Hydrol., 432–433, 154–167, https://doi.org/10.1016/j.jhydrol.2012.02.014, 2012.

755

Grosjean, A., Pittet, B., Gardien, V., Leloup, P., Mahéo, G., and Barraza Garcia, J.: Tectonic heritage in drainage pattern and dynamics: the case of the French South Alpine Foreland Basin (ca.45–20 Ma), Basin Res., 29, 26–50, https://doi.org/10.1111/bre.12134, 2017.

Hartmann, A., Lange, J., Vivó Aguado, À., Mizyed, N., Smiatek, G., and Kunstmann, H.: A multi-model approach for improved simulations of future water availability at a large Eastern Mediterranean karst spring, J. Hydrol., 468–469, 130–138, https://doi.org/10.1016/j.jhydrol.2012.08.024, 2012.

760

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

Herrnegger, M., Nachtnebel, H.-P., and Haiden, T.: Evapotranspiration in high alpine catchments – an important part of the water balance!, Hydrol. Res., 43, 460–475, https://doi.org/10.2166/nh.2012.132, 2012.

765

Huet, B., Lasseur, E., Briais, J., Bellahsen, N., Loget, N., Rubino, J.-L., and Suc, J.-P.: Palaeoenvironmental and drainage network evolution of the Oligocene Western Alpine Foreland Basin, https://doi.org/10.5194/egusphere-egu23-7337, 2023.

Institut national de l'information géographique et forestière (IGN), 2024: Lidar HD data – Classified point cloud download, available at: https://geoservices.ign.fr/lidarhd, last access: 15 October 2025.

770

Ismail, M. F., Bogacki, W., Disse, M., Schäfer, M., and Kirschbauer, L.: Estimating degree-day factors of snow based on energy flux components, The Cryosphere, 17, 211–231, https://doi.org/10.5194/tc-17-211-2023, 2023.

Jagercikova, M., Mocochain, L., Lebatard, A.-E., Bourles, D. L., Leanni, L., Sartégou, A., and Zappelli, A.: Découverte et étude de remplissages karstiques allochtones d’âge miocène dans l’Obiou (Dévoluy, Alpes françaises) Implications géomorphologiques et paléogéographiques, Karstologia, 2021.

775

Jeannin, P.-Y., Artigue, G., Butscher, C., Chang, Y., Charlier, J.-B., Duran, L., Gill, L., Hartmann, A., Johannet, A., Jourde, H., Kavousi, A., Liesch, T., Liu, Y., Lüthi, M., Malard, A., Mazzilli, N., Pardo-Igúzquiza, E., Thiéry, D., Reimann, T., Schuler, P., Wöhling, T., and Wunsch, A.: Karst modelling challenge 1: Results of hydrological modelling, J. Hydrol., 600, 126508, https://doi.org/10.1016/j.jhydrol.2021.126508, 2021.

Jódar, J., González-Ramón, A., Martos-Rosillo, S., Heredia, J., Herrera, C., Urrutia, J., Caballero, Y., Zabaleta, A., Antigüedad, I., Custodio, E., and Lambán, L. J.: Snowmelt as a determinant factor in the hydrogeological behaviour of high mountain karst aquifers: The Garcés karst system, Central Pyrenees (Spain), Sci. Total Environ., 748, 141363, https://doi.org/10.1016/j.scitotenv.2020.141363, 2020.

780

Lafaysse, M., Cluzet, B., Dumont, M., Lejeune, Y., Vionnet, V., and Morin, S.: A multiphysical ensemble system of numerical snow modelling, The Cryosphere, 11, 1173–1198, https://doi.org/10.5194/tc-11-1173-2017, 2017.

Lemaitre-Basset, T., Oudin, L., and Thirel, G.: Evapotranspiration in hydrological models under rising CO2: a jump into the unknown, Clim. Change, 172, 36, https://doi.org/10.1007/s10584-022-03384-1, 2022a.

785

Lemaitre-Basset, T., Oudin, L., Thirel, G., and Collet, L.: Unraveling the contribution of potential evaporation formulation to uncertainty under climate change, Hydrol. Earth Syst. Sci., 26, 2147–2159, https://doi.org/10.5194/hess-26-2147-2022, 2022b.
<figure><img src="image_47.png" /></figure>

790 Lemot, F., Valla, P. G., van der Beek, P., Jagercikova, M., Niedermann, S., Carcaillet, J., Sobel, E. R., Andò, S., Garzanti, E., Robert, X., Balvay, M., Bernet, M., Glodny, J., and Mocochain, L.: Miocene cave sediments record topographic, erosional and drainage development in the Western European Alps, Earth Planet. Sci. Lett., 621, 118344, https://doi.org/10.1016/j.epsl.2023.118344, 2023.

795 Luparini, V.: Etude hydrogéologique du massif du Dévoluy, Hautes-Alpes - Alpes françaises, Faculté des Sciences de l'Université de Grenoble, France, 142 pp., 1975.

Magand, C., Ducharne, A., Le Moine, N., and Brigode, P.: Parameter transferability under changing climate: case study with a land surface model in the Durance watershed, France, Hydrol. Sci. J., 60, 1408–1423, https://doi.org/10.1080/02626667.2014.993643, 2015.

Maire, R.: Hydrogéologie du plateau de Bure, Voconcies, 9–16, 1980.

800 Marson, P., Corre, L., Soubeyroux, J.-M., Sauquet, E., Robin, Y., Vrac, M., and Dubois, C.: Explore2 – Rapport de synthèse sur les projections climatiques régionalisées, report, METEO FRANCE ; INRAE ; Institut Pierre-Simon Laplace, 2024.

Martinec, J. and Rango, A.: Parameter values for snowmelt runoff modelling, J. Hydrol., 84, 197–219, https://doi.org/10.1016/0022-1694(86)90123-X, 1986.

805 Mazzilli, N., Guinot, V., Jourde, H., Lecoq, N., Labat, D., Arfib, B., Baudement, C., Danquigny, C., Dal Soglio, L., and Bertin, D.: 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, 2019.

Mazzilli, N., Sivelle, V., Cinkus, G., Jourde, H., and Bertin, D.: KarstMod User Guide - version 3.0, https://doi.org/10.1016/j.envsoft.2017.03.015, 2023.

810 Meckel, L. D., Ford, M., and Bernoulli, D.: Tectonic and sedimentary evolution of the Dévoluy Basin, a remanent of the Tertiary western Alpine foreland basin, SE France, Géologie Fr., 2, 3–26, 1996.

Meeks, J. and Hunkeler, D.: Snowmelt infiltration and storage within a karstic environment, Vers Chez le Brandt, Switzerland, J. Hydrol., 529, 11–21, https://doi.org/10.1016/j.jhydrol.2015.06.040, 2015.

Monjuvent, G.: Le Drac, morphologie, stratigraphie et chronologie quaternaires d'un bassin alpin, phdthesis, Université Paris-Diderot - Paris VII, 1979.

815 Oudin, L.: Une formule simple d'évapotranspiration potentielle pour la modélisation pluie-débit à l'échelle du bassin versant, Houille Blanche, 92, 113–120, https://doi.org/10.1051/lhb:2006109, 2006.

Parizot, O., Espurt, N., Viguier, B., Duclaux, G., Zappelli, A., Lahfid, A., Rispal, N., Mocochain, L., Henry, P., Jagercikova, M., Audra, P., Arfib, B., and Cadilhac, L.: 3D Structural Modelling and Restoration of a Deformed Alpine Karst Reservoir: Insights into the Groundwater Flows of the Dévoluy Massif (French Alpine Foreland), Basin Res., 37, e70065, https://doi.org/10.1111/bre.70065, 2025.

820 Perrin, C., Michel, C., and Andréassian, V.: Improvement of a parsimonious model for streamflow simulation, J. Hydrol., 279, 275–289, https://doi.org/10.1016/S0022-1694(03)00225-7, 2003.

Pool, S., Vis, M., and Seibert, J.: Evaluating model performance: towards a non-parametric variant of the Kling-Gupta efficiency, Hydrol. Sci. J., 63, https://doi.org/10.1080/02626667.2018.1552002, 2018.
<figure><img src="image_48.png" /></figure>

825 Porthault, B.: Le Crétacé supérieur de la Fosse vocontienne et des régions limitrophes: (France sud-est), micropaléontologie, stratigraphie, paléogéographie, Université Claude Bernard - Lyon I, Lyon, 265 pp., 1974.

Porthault, B.: Excursion dans la "Fosse vocontienne". Le Crétacé supérieur de la "Fosse vocontienne" et des régions voisines. Excursion du Groupe Français du Crétacé., France, 1976.

830 Pulka, T., Herrnegger, M., Ehrendorfer, C., Lücking, S., Avanzi, F., Formayer, H., Schulz, K., and Koch, F.: Evaluating precipitation corrections to enhance high-alpine hydrological modeling, J. Hydrol., 645, 132202, https://doi.org/10.1016/j.jhydrol.2024.132202, 2024.

Reisch, C. E. and Toran, L.: Characterizing snowmelt anomalies in hydrochemographs of a karst spring, Cumberland Valley, Pennsylvania (USA): evidence for multiple recharge pathways, Environ. Earth Sci., 72, 47–58, https://doi.org/10.1007/s12665-013-2935-5, 2014.

835 Ribes, A., Boé, J., Qasmi, S., Dubuisson, B., Douville, H., and Terray, L.: An updated assessment of past and future warming over France based on a regional observational constraint, Earth Syst. Dyn., 13, 1397–1415, https://doi.org/10.5194/esd-13-1397-2022, 2022.

Rispal, N., Audra, P., Mocochain, L., Jagercikova, M., Henry, P., and Demory, F.: Understanding the sedimentary dynamics of alpine caves: the example of the Dévoluy Massif (France), 2025.

840 Ruelland, D.: Development of the snow- and ice-accounting routine (SIAR), J. Hydrol., 624, 129867, https://doi.org/10.1016/j.jhydrol.2023.129867, 2023.

Ruelland, D.: Potential of snow data to improve the consistency and robustness of a semi-distributed hydrological model using the SAFRAN input dataset, J. Hydrol., 631, 130820, https://doi.org/10.1016/j.jhydrol.2024.130820, 2024.

845 Sauquet, E., Evin, G., Siauve, S., Aissat, R., Arnaud, P., Bérel, M., Bonneau, J., Branger, F., Caballero, Y., Colléoni, F., Ducharne, A., Gailhard, J., Habets, F., Hendrickx, F., Héraut, L., Hingray, B., Huang, P., Jaouen, T., Jeantet, A., Lanini, S., Le Lay, M., Magand, C., Mimeau, L., Monteil, C., Munier, S., Perrin, C., Robelin, O., Rousset, F., Soubeyroux, J.-M., Strohmenger, L., Thirel, G., Tocquer, F., Tramblay, Y., Vergnes, J.-P., and Vidal, J.-P.: A large transient multi-scenario multi-model ensemble of future streamflow and groundwater projections in France, EGUSphere, 1–41, https://doi.org/10.5194/egusphere-2025-1788, 2025.

850 Seibert, J. and Vis, M. J. P.: Teaching hydrological modeling with a user-friendly catchment-runoff-model software package, Hydrol. Earth Syst. Sci., 16, 3315–3325, https://doi.org/10.5194/hess-16-3315-2012, 2012.

Sexstone, G. A., Clow, D. W., Fassnacht, S. R., Liston, G. E., Hiemstra, C. A., Knowles, J. F., and Penn, C. A.: Snow Sublimation in Mountain Environments and Its Sensitivity to Forest Disturbance and Climate Warming, Water Resour. Res., 54, 1191–1211, https://doi.org/10.1002/2017WR021172, 2018.

855 Sivelle, V., Cinkus, G., Mazzilli, N., Labat, D., Arfib, B., Massei, N., Cousquer, Y., Bertin, D., and Jourde, H.: Improvement of the KarstMod modelling platform for a better assessment of karst groundwater resources, Hydrol. Earth Syst. Sci., 29, 1259–1276, https://doi.org/10.5194/hess-29-1259-2025, 2025.

Soubeyroux, J.-M., Dubuisson, B., Bernus, S., Samacoïts, R., Rousset, F., Schneider, M., Drouin, A., Madec, T., Tardy, M., and Corre, L.: A quel climat s’adapter en France selon la TRACC ?, Meteo-France, 2024.
<figure><img src="image_49.png" /></figure>

860 Taucare, M., Daniele, L., Viguier, B., Vallejos, A., and Arancibia, G.: Groundwater resources and recharge processes in the Western Andean Front of Central Chile, Sci. Total Environ., 722, 137824, https://doi.org/10.1016/j.scitotenv.2020.137824, 2020.

Valéry, A.: Modélisation précipitations débit sous influence nivale: Elaboration d'un module neige et évaluation sur 380 bassins versants, Institut des Sciences et Industries du Vivant et de l'Environnement, AgroParisTech, 418 pp., 2010.

865 Valéry, A., Andréassian, V., and Perrin, C.: 'As simple as possible but not simpler': What is useful in a temperature-based snow-accounting routine? Part 2 – Sensitivity analysis of the Cemaneige snow accounting routine on 380 catchments, J. Hydrol., 517, 1176–1187, https://doi.org/10.1016/j.jhydrol.2014.04.058, 2014a.

Valéry, A., Andréassian, V., and Perrin, C.: 'As simple as possible but not simpler': What is useful in a temperature-based snow-accounting routine? Part 2 – Sensitivity analysis of the Cemaneige snow accounting routine on 380 catchments, J. Hydrol., 517, 1176–1187, https://doi.org/10.1016/j.jhydrol.2014.04.058, 2014b.

Vidal, J.-P., Martin, E., Franchistéguy, L., Baillon, M., and Soubeyroux, J.-M.: A 50-year high-resolution atmospheric reanalysis over France with the Safran system, Int. J. Climatol., 30, 1627–1644, https://doi.org/10.1002/joc.2003, 2010.

Viviroli, D., Dürr, H. H., Messerli, B., Meybeck, M., and Weingartner, R.: Mountains of the world, water towers for humanity: Typology, mapping, and global significance, Water Resour. Res., 43, https://doi.org/10.1029/2006WR005653, 2007.

875 Zappelli, A., Belleville, A., and Jagercikova, M.: Hydrogeoly of Dévoluy karstic system: New insights with dye tracing and rainfall-discharge model, 2018.

Zierl, B. and Bugmann, H.: Global Change Impacts on Hydrological Processes in Alpine Catchments, Water Resour. Res. - WATER RESOUR RES, 410, https://doi.org/10.1029/2004WR003447, 2005.