Nossana Spring (Northern Italy) under Climate Change: Projections of Future Discharge Rates and Water Availability

Abstract

Nossana represents an important pre-Alpine karst spring for drinking supply, sustaining a water distribution system serving 300,000 people. The goal of this study was to project Nossana discharges and evaluate potential supply limits for four future periods (2021–2040, 2041–2060, 2061–2080, 2081–2100). Bias-corrected Regional Climate Models (RCMs), part of the EURO-CORDEX experiment and forced by three emission scenarios (RCP2.6, RCP4.5, RCP8.5), were evaluated, statistically downscaled, and used as input in a calibrated rainfall-runoff model ensemble. For each emission scenario, the calibrated model ensemble considered three RCMs and ten rainfall-runoff parameterizations. Projected ensemble mean discharges are lower than observations (3%–23%) for all RCPs, though they do not show a clear trend between the four time periods. Days characterized by discharges lower than actual water demand are projected to decrease, except for the RCP8.5 emission scenarios and the period 2081–2100. Conversely, the same consecutive days are expected to increase after 2060 for all emission scenarios. These results reflect the projected precipitation trend, characterized by longer, drier summer periods and wetter autumns in comparison to today’s climate. Also, they indicate a possible need for alternative drinking water resources. The proposed methodology was demonstrated to deliver useful quantitative information for water management in the mid- long-term period.

1. Introduction

Karst springs play a fundamental role in large-scale human water supply, both from a strategic and socio-economic point of view [1]. Quantifying actual and future project spring discharges is extremely important to manage water resources in karstic areas, especially in view of the effects related to climate change. According to the Intergovernamental Panel on Climate Change (IPCC), the mean temperature is expected to increase globally throughout the 21st century, with a consequent increase in the frequency and magnitude of heat waves, drought periods and intense precipitation events in many regions of the world [2]. The prediction of the response to climate change in hydrogeological systems becomes vital, especially if exploited for drinking water supply [3] or the sustainability of groundwater dependent ecosystems [4].

The application of traditional groundwater flow models for the prediction of fractured media is complicated due to the duality of the flow systems [5]; this is particularly true for karst systems [6,7]. A typical karst system is formed by a fractured rock matrix, which includes both diagenetic rock micro pores and small fractures of tectonic origin, as well as a network of widely articulated karst conduits [1,8]. This dual nature leads to a hybrid flow behavior, Darcian within the rock matrix, and turbulence in the conduits. In addition, flow in the conduits can be assimilated to full pipe (pressurized) flow or to free surface flow, based on the channel saturation. Attempts to define physically based, distributed models incorporating some or all of these characteristics have been made by many authors [9,10,11,12,13,14,15]. However, the definition of a conceptual model that could depict the system heterogeneity reliably, in terms of both 3D geometry and parameterization, can be very data demanding, time consuming, and costly [16]. Speleological exploration (e.g., [17]) and artificial tracers (e.g., [18]) can furnish information on spring catchment areas and the main geometric properties of the karst system. Additional insights regarding different water sources and water residence time in the aquifer can be derived from natural tracer analyses, such as water isotopes, major ions, trace elements, and dissolved organic carbon (e.g., [19,20,21,22,23]). In-situ hydraulic methods, such as pumping tests, can provide knowledge on the degree of confinement of the aquifer and its quantitative parameters (e.g., [24,25]). Non-invasive geophysical methods can also be used to gather information on aquifer geometry and boundaries (e.g., electric and seismic methods, [26,27]), to identify major voids and conduits (e.g., gravimetric and geomagnetical methods, [28]), and to recognize preferential infiltration pathways (e.g., self-potential methods, [29]). In addition, to allow model calibration, a continuous monitoring of at least precipitation and spring discharge is essential [16].

When the main goal of the study concerns the response of a karst spring to meteorological events, the continuous monitoring of precipitation and discharge could suffice for a modeling setup. Input-output models can be used to relate the precipitation contribution to the outflow through empirical equations based on lumped parameters (e.g., [30,31]), neural networks (e.g., [32]), or input-response functions (e.g., [33,34]). An increasing number of studies involving such modelling approaches have been documented for karst environments in recent years (e.g., [35,36,37,38,39]). Due to the lack of physical processes representation, a limitation of these models concerns their loss of reliability if applied in conditions that largely differ from calibration [40,41]. Parameter transferability has been largely tested (e.g., [42,43,44]), suggesting an increase in simulation errors and uncertainties with increasing differences in mean rainfall. Also, performance losses appear larger when transferring parameters from wet to dry periods than vice versa. Vaze et al. [42] concluded that calibrated models could be used for climate change impact studies if mean rainfall changes between future and calibration periods fall within a −15% to +20% range.

Climate projections are derived from Global Circulation Models (GCMs) forced according to the representative concentration pathways (RCP, [45]) proposed by the IPCC [2]. Usually, these models are run at large horizontal resolutions (≥ 80 km) and are therefore not able to reproduce climate states at regional or local scales reliably [46]. To obtain climate projections at a finer resolution by reproducing the physics of the processes, dynamical downscaling through Regional Climate Models (RCMs) is performed (e.g., [47,48]). With the ongoing climate modeling experiments—for example, the Coordinated Regional Climate Downscaling Experiment (CORDEX, http://www.cordex.org/), WorldClim dataset project (http://www.worldclim.org/), and the Climate Change, Agriculture, and Food Security research program (CCAFS, http://ccafs-climate.org/downscaling/)—future climate projections are made available over several domains worldwide with a horizontal resolution of up to 12.5 km. Some studies (e.g., [49,50]) demonstrated that dynamical downscaling can be applied to achieve an even larger detail (≤4 km). However, these methods are computationally very intensive and often show performance losses over complex terrains [46,51]. The computationally inexpensive alternative is statistical downscaling, which provides a predictive framework of local conditions (with almost no bias with historical records) through statistical relationships with regional climatic aspects. The hypothesis is that the relationship between large-scale atmospheric circulation and local scale dynamics is constant [52,53,54]. Therefore, statistical downscaling allows for producing large ensembles of climate change projections for extended time periods [46,55].

The overarching goal of this study is to propose a methodological approach to evaluate possible variations in the discharge regime of a karst spring in basins affected by climate change. The selected study area is the hydrogeological basin of the Nossana Spring (northern Italy), which is heavily exploited for drinking water. Using a control data period covering 1998–2017, the specific objectives of this study are: (i) quantification of the expected changes in precipitation and temperature in the study area, for twenty-year periods till 2100, as resulting from the statistical downscaling of RCMs run under three RCPs (2.6, 4.5, 8.5); (ii) calibration and validation of a lumped-parameter model (GR4J with CemaNeige [56,57,58]) based on observed data; (iii) recognition of possible limits in the future utilization of the spring as a drinking supply based on simulated discharge and actual warning thresholds.

2. Study Area

The Nossana Spring is located in the central Prealps of Seriana Valley, within the Lombardy Region of the Province of Bergamo (Northern Italy) (Figure 1). Its hydrogeological basin covers an area of about 80 km² and is characterized by large differences in altitude, ranging from 474 m a.s.l. (Nossana Spring) to 2512 m a.s.l. (Pizzo Arera Mountain).

The spring is managed by the public company UniAcque S.p.A. and feeds the main local water distribution system, serving over 300,000 people. The observed mean annual discharge (1998–2017) is 3.77 m³ s⁻¹. For the same period, the observed daily variability is extremely high, with values ranging from 0.55 m³ s⁻¹ to 18.00 m³ s⁻¹. In a single day, discharge varied up to 12 m³ s⁻¹, suggesting high flow velocities characterizing the system.

To manage the water resource, UniAcque S.p.A. set two warning thresholds based on average water demand according to the period of the year (one threshold for the cold-wet season and one threshold for the hot-dry months). Below these discharge thresholds, there is a need to integrate the spring resources pumping water from deep wells in the Serio River Valley. This entails different issues, from costs increase to quality decrease. Moreover, nowadays the additional discharges are limited (0.50 m³ s⁻¹). During wet, cold periods, the warning threshold is set at 1.32 m³ s⁻¹; during dry, hot periods, the threshold is set at 1.52 m³ s⁻¹. These two discharge rates include a share for UniAcque S.p.A. (0.80 and 1.00 m³ s⁻¹, respectively) and a share for environmental flow and use of the downstream municipalities (0.52 m³ s⁻¹). Water that is not used by UniAcque S.p.A. flows into the Serio River, the main watercourse of the Seriana Valley.

Many studies have defined a hypothetical hydrogeological catchment for the spring [59,60], considering the geological, geomorphological, and tectonic evidences of the area [61,62,63,64]. The karst massif that characterizes the area consists of calcareous and calcareous-silicoclastic series (medium-late Triassic). The main calcareous formation is represented by “Calcare di Esino” (middle-late Trias), constituted by deeply karstified massive or stratified limestone, with a maximum thickness of about 800 m. Karst evidences (e.g., dolines and tunnels) are mainly recognized in the western and north-eastern area of the basin. The late-Triassic marl formations with low permeability (Breno Formation, Metamorphic Bergamasco Limestone and Gorno Formation) link the massif to its southern, western, and eastern borders along the tectonic contacts with the Calcare di Esino formation. Three main fragile structures mark the limits of the massif along the northern side (“Valtorta-Valcanale thrust”), west side (“Grem fault”), and southeast side (“Clusone fault”). More specific geological information about the study area can be found in literature [63,65,66].

According to Vigna and Banzato [67], based on the geological and hydrographic characteristics described above, Nossana can be classified as a karst spring fed by a dominant drainage system. These systems are usually characterized by a high permeability, which is linked to a considerable karstification process. The outflow network is well-organized with a series of main natural tunnels and secondary conduits that rapidly discharge incoming infiltration water. The most important features are the very low or complete absence of a classic phreatic zone (indeed, the presence of a phreatic zone has not been verified for Nossana yet) and the very high flow velocity.

The climate of the study area is wet and temperate. According to data recorded at the meteorological station of Clusone (591 m a.s.l.), located at around 4 km east of Nossana Spring, rainfall is distributed throughout the year. Mean annual rainfall calculated for the period 1998–2017 is 1384 mm, with November as the wettest month (average of 160 mm) and January as the driest (average of 65 mm). May, June, October, and November yield similar amounts of rainfall. Conversely, the winter lowest values indicate a change from liquid (rainfall) to solid precipitation (snow/ice). Average monthly temperatures vary between 1.4 °C (January) and 21.7 °C (July), with an overall mean annual temperature of 11.5 °C. According to Köppen–Geiger climate classification [68], the climate of the area can be classified as Cfb; namely, temperate with a wet summer, the mean temperature of the warmest month below 22 °C, and at least four months above an average of 10 °C.

3. Data

Data used to perform the study were derived from three main sources. Daily spring discharge values (m³ s⁻¹) were provided by UniAcque S.p.A. for the period 1998–2017. Observed daily precipitation (mm), minimum and maximum temperature (°C) time series (1998–2017) were downloaded for the Clusone meteorological station from the web portal of the Regional Agency for the Environment (ARPA Lombardia, https://www.arpalombardia.it/Pages/Meteorologia/Richiesta-dati-misurati.aspx). In addition, daily time series of precipitation, minimum temperature, and maximum temperature resulting from climate simulations were extracted for the period 1970 to 2100 from selected RCM model runs of the EURO-CORDEX database. The latter were accessed through the Earth System Grid Federation (ESGF https://esgf-data.dkrz.de). In detail, bias corrected RCM model runs available at 12.5 km resolution for all emission pathways (RCP2.6, RCP4.5, RCP8.5) were selected (

Table 1. Summary of the climate model data sources used for this study.

ID	Global Circulation Model	Ensemble Member	Regional Climate Model	Institution
Mod_1	ICHEC-EC-EARTH	r12i1p1	RCA4_v1	Swedish Meteorological and Hydrological Institute (SMHI)
Mod_2	MPI-ESM-LR	r2i1p1	REMO2009_v1	Max Planck Institute for Meteorology, Climate Service Center
Mod_3	MPI-ESM-LR	r1i1p1	REMO2009_v1	Max Planck Institute for Meteorology, Climate Service Center

). This selection was done to allow a robust and coherent comparison between the projected future emission pathways, which depend on the political and economic actions that will be undertaken to reduce greenhouse gas emissions. The RCP2.6 scenario foresees that emissions will be halved from the year 2050; for RCP 4.5 this target will be reached in 2080; while according to RCP8.5 emissions will continue to increase with today’s rate [2]. The constraint to have the same models for all the RCPs limited their number; however, they include two GCMs, three different ensemble runs, and two RCMs. All models were corrected for bias by the Swedish Meteorological and Hydrological Institute (SMHI) using a distribution-based scaling method [69] based on MESAN reanalysis data (1989–2010).

4. Methods

The study is based on four main steps, which are summarized in Figure 2. First, climate models are evaluated to verify their ability to reproduce the general climatology of the study area. Second, RCMs future projections are downscaled through statistical methods and local climate changes are evaluated. Third, observed time series of precipitation and temperature are used as input in a rainfall-runoff hydrologic model, which is calibrated and validated against observed discharge. Finally, the downscaled time series are used to feed the calibrated hydrologic model and possible changes in the spring discharge rate characteristics are analyzed.

4.1. Climate Models Evaluation

Daily time series of precipitation, minimum temperature, and maximum temperature were extracted for the three RCMs realizations and the three RCPs (total of nine combinations) from the grid cell where the Clusone meteorological station is located. The RCMs were evaluated by comparing their results for the control period (1998–2017) with the available observational data. Representative Concentration Pathways begin to diverge after 2005; therefore, between this date and the end of the control period, the same model, under different RCPs, could produce data with small differences. For this reason, the evaluation included all nine modeling combinations.

The evaluation consisted of the comparison of monthly rainfall and temperature climatology, and of the verification of their similarity through two performance indices. The selected indices were the Nash-Sutcliffe Efficiency (NSE—[70]) and the Mean Absolute Error (MAE). The RCMs were considered to be generally able to catch the local climatology for positive (>0.8) NSE values and MAE values lower than 20% of the mean monthly rainfall (temperature). The indices are calculated as:

$$\mathrm{NSE}=1.0-\frac{{{\displaystyle \sum }}_{\mathrm{t}=1}^{\mathrm{N}}{\left({\mathrm{OBS}}_{\mathrm{t}}-{\mathrm{MOD}}_{\mathrm{t}}\right)}^2}{{{\displaystyle \sum }}_{\mathrm{t}=1}^{\mathrm{N}}{\left({\mathrm{OBS}}_{\mathrm{t}}-\overline{\mathrm{OBS}}\right)}^2},$$

(1)

$$\mathrm{MAE}=\frac{1}{\mathrm{N}}{{\displaystyle \sum }}_{\mathrm{t}=1}^{\mathrm{N}}\left|{\mathrm{OBS}}_{\mathrm{t}}{-\mathrm{MOD}}_{\mathrm{t}}\right|,$$

(2)

where OBS refers to observations, MOD to modelled data, t indicates the time step (in this case the month), and N is the total number of time steps.

4.2. Statistical Downscaling

Precipitation and temperature data for the nine modeling combinations were statistically downscaled using a two-step approach, based on change factors and weather simulators (e.g., [71]). Assuming a stable climatology for twenty-year periods, the downscaling was performed for four time intervals (2021–2040; 2040–2060; 2061–2080; 2081–2100). Change factors [72] were calculated from daily modeled data, on a monthly (half-monthly) basis, for a selected set of rainfall (temperature) statistics, for each modeling combination and time period. For details on change factors calculation of different rainfall and temperature statistics, the reader can refer to Kilsby et al. [73]. Once calculated, change factors were applied to the same set of observed statistics to obtain values describing the local scale future behavior for each month of the year. These local scale, future statistics were used as inputs to fit a stochastic rainfall and temperature generator returning daily time series. Rainfall, minimum temperature, and maximum temperature time series were generated for each modeling combination and time interval.

4.2.1. Rainfall Generator

For future precipitation simulation, RainSim V3.1.1 [74] was used. RainSim is a rainfall generator software able to work in both single site and spatial-temporal mode. It is based on a Neymann-Scott Rectangular Pulses (NSRP) model [75,76], which manages rainfall occurrence and amount as a single process. Storms arrive as Poisson processes and each storm is characterized by a random number of rain cells, each having a certain origin delay time from the beginning of the storm, a certain duration, and a certain intensity. In single site mode—the one used in this study—to generate rainfall time series, RainSim needs to fit five internal parameters for each calendar month (

Table 2. RainSim parameters for time series generation in single site mode [74].

Parameter	Description	Units
λ−1	Mean time between adjacent storm origins	[h]
β−1	Mean waiting time for raincell origins after storm origin	[h]
η−1	Mean duration of raincell	[h]
Ν	Mean number of raincells per storm	[-]
ξ−1	Mean intensity of a raincell	[mm·h−1]

). These parameters are fitted through an optimization function embedded in the software, which works based on a set of input statistics (e.g., for future periods those obtained from the change factor approach described in Section 4.2). For each month, seven statistics were used to fit the model: mean daily rainfall, daily variance, daily skew, lag-1 autocorrelation, probability of dry days (threshold 0.2 mm), probability of dry day to dry day transition, and probability of wet day to wet day transition.

4.2.2. Temperature Generator

To generate maximum and minimum temperature time series, the approach presented by Richardson [77] was adopted and implemented in the R software (www.r-project.org). The R-scripts were first developed by Camera et al. [78] and adjusted for this study. Temperature was modeled as a stochastic process with the daily mean and standard deviation conditioned on the wet and dry state of the day (threshold 0.2 mm). Richardson [77] proposed an approach in which residual time series are generated, so the first step consisted of computing the residuals for the minimum and maximum temperature time series recorded in Clusone. To do so, means and standard deviations of the 20-year time series were calculated for dry and wet days for each half-month. Going forward, a first-order Fourier series was fitted through these values, using the least squares method to obtain mean and standard deviation values for each calendar day (1–365) for both dry and wet states.

As in [77], a weak stationary generating process [79] was used to simulate time series. Only one residual time series was generated for each variable (minimum and maximum temperature) because in working with residuals the effect of the wet or dry state of the day was removed. Once the residuals time series were obtained, the daily values of minimum and maximum temperature were calculated multiplying the residual by the standard deviation and adding the mean, according to the state of day. Before back-transforming the generated residual time series into actual temperature values, standard deviation and mean were modified by applying the calculated change factors. In this way, the simulated future temperature time series were linked to the rainfall series previously generated.

4.2.3. Evaluation of the Local Climate Change Signature

A visual comparison between the monthly climatology calculated from observations and those derived from the generated time series was carried out for each time interval (2021–2040, 2041–2060, 2061–2080, 2081–2100) and RCP, to verify possible shifts in the general climate trend of the area. In addition, changes in mean annual precipitation and temperature were calculated. It was verified that mean annual rainfall changes projected by the statistical downscaling of every modelling combination fell between −15% and +20% in comparison to the control period. This is the range suggested by Vaze et al. [42] as a limit for a solid and reliable temporal transferability of lumped model parameter sets.

4.3. Hydrologic Modelling

To model the spring response to precipitation events, the daily rainfall-runoff (here applied as rainfall-discharge) four-parameter GR4J model (Génie Rural Journalier with 4 parameters, [56]), extended with the CemaNeige snow accounting routine [57,58,80], was applied. Model calibration, validation, and future discharge simulations were performed using the open-source AirGR 1.4.3.52 package [81,82] in the R 3.5.3 software (www.r-project.org). GR4J was selected since it is under continuous development and widely known (e.g., [83,84]), as well as being characterized by a parsimonious but flexible structure that allows for the simulation of complex catchments [85]. Also, it was already applied and verified in karst environments (e.g., [38,85,86]).

4.3.1. Model Description

A schematic representation of the model is presented in Figure 3. Model inputs are time series of daily precipitation (P) and daily potential evapotranspiration (E). The model structure consists of two non-linear reservoirs, the production store and the routing store, a delay function element, and a groundwater exchange component. From the inputs, a net rainfall (Pn) or a net evapotranspiration (En) is calculated. In case of net evapotranspiration, a quantity Es (actual evapotranspiration) is subtracted from the production store. In case of net rainfall (Pn), this is split into two components. The first component (Ps) is added to the production store, while the second component is added to the release from the production store (Perc) to give a total amount Pr. Ps, Es, and Perc are defined as functions of the production store level (S) and of its maximum capacity (X1). Pr is then split into two fixed amounts; 0.9Pr is routed by the unit hydrograph HU1 so that a quantity Q9 enters the routing store, while 0.1Pr is routed through the unit hydrograph HU2 so that a quantity Q1 is available for fast drainage towards the basin outlet. The unit hydrographs are defined based on the model parameter X4. Before reaching the outlet, Q1 could lose or gain water from interactions with groundwater, resulting in Qd. Water is also released towards the basin outlet from the routing store (Qr) based on the store level (R), its maximum capacity (X3), and interactions with groundwater. The latter are controlled by the parameter X2 (exchange coefficient). The model output discharge (Q) is given by the sum of Qr and Qd. In a karst environment the roles played by the production store and the routing store could be assimilated to those of the rock matrix and the conduits, respectively.

The CemaNeige module allowed us to account for delays in the GR4J precipitation input, transforming part of the precipitation into snowfall and simulating snow melting according to temperature data and two parameters. Model inputs are daily precipitation, minimum temperature, and maximum temperature time series, which need to be referred to a specific elevation. In addition, the distribution of the elevations within the study area is necessary. Based on these data, with an internal routine, the model divides the study area in elevation layers of equal areas and the user can optionally define a mean annual snowfall value for each. For layers with median elevations lower than 1500 m a.s.l., CemaNeige partitions precipitation (P_RS) into rainfall (RAIN) and snow (SNOW) as follows:

$$\mathrm{RAIN}=\{\begin{array}{c}{\mathrm{P}}_{\mathrm{RS}} \mathrm{if} \mathrm{Tmax}>0 °\mathrm{C}; \mathrm{Tmin}>0 °\mathrm{C}\\ 0 \mathrm{if} \mathrm{Tmax}\le 0 °\mathrm{C}; \mathrm{Tmin}\le 0 °\mathrm{C}\\ {\mathrm{P}}_{\mathrm{RS}}·\frac{\mathrm{Tmax}}{\mathrm{Tmin}-\mathrm{Tmax}} \mathrm{if} \mathrm{Tmax}>0 °\mathrm{C}; \mathrm{Tmin}<0 °\mathrm{C}\end{array}$$

(3)

$$\mathrm{SNOW}={\mathrm{P}}_{\mathrm{RS}}-\mathrm{RAIN}$$

(4)

where Tmax and Tmin are daily maximum and minimum temperature, respectively. For layers with median elevations equal or higher to 1500 m a.s.l., the model calculates the rainfall fraction as follows:

$$\mathrm{RAIN}=\{\begin{array}{c}{\mathrm{P}}_{\mathrm{RS}} \mathrm{if} \mathrm{Tmax}>3 °\mathrm{C}; \mathrm{Tmin}>3 °\mathrm{C} \\ 0 \mathrm{if} \mathrm{Tmax}\le -1 °\mathrm{C}; \mathrm{Tmin}\le -1 °\mathrm{C} \\ {\mathrm{P}}_{\mathrm{RS}}·\frac{\mathrm{Tmax}}{\mathrm{Tmin}-\mathrm{Tmax}} \mathrm{if}-1 °\mathrm{C}<\mathrm{Tmax} \mathrm{and} \mathrm{Tmin}\le 3 °\mathrm{C} \end{array}$$

(5)

The snow fraction is calculated according to Equation (4) in this case also. Melting is then controlled by temperature and the two model parameters: Kf [-], which is the degree-day melting factor, and Ctg [mm° C⁻¹ day⁻¹], namely the weighting coefficient for snowpack thermal state, which describes the thermal inertia of the accumulated snow. The higher is Kf, the larger is the snowmelt; the higher is Ctg, the later snow melts. Recently, Riboust et al. [80] presented a modified and extended version of the CemaNeige module to account for the melting heterogeneity of the snow. To do so, they enriched the model with two additional parameters to simulate hysteresis cycles with a linear model. The two extra parameters define the accumulation threshold at which the Snow Cover Area (SCA) equals 100% and the melting threshold at which the SCA drops below 100%. In this study, both the standard and the hysteretic CemaNeige modules were tested.

4.3.2. Model Calibration and Validation

Regarding the inputs, potential evapotranspiration was calculated using minimum and maximum temperature records applying the Hargreaves equation [87]. For the CemaNeige module, the elevation of the Clusone meteorological station (591 m a.s.l.) was specified as the reference one. Based on the elevation distribution within the study area, five equal area layers were created (458–882, 882–1222, 1222–1518, 1518–1789, 1789–2507 m a.s.l.). For each layer, the mean annual snowfall value (45 mm, 60 mm, 100 mm, 155 mm, 200 mm, respectively) was derived, in water equivalent, from the available national and regional snowfall maps [88,89]. For the conversion between snowfall and water equivalent, the conventional 10:1 water density to snow density proportion (1 cm of snow equals 1 mm of rainfall) was applied [90].

The observation period (1998–2017) was split into three intervals: two years for model spin-up (1998–1999) and nine years each for calibration and validation (2000–2008 and 2009–2017, respectively). Ten thousand random parameter sets were generated within the ranges specified in

Table 3. Parameter ranges used to generate the random samples for model calibration and validation.

GR4J
Parameter	Description	Calibration Range
X1	Production store capacity [mm]	100–1200
X2	Intercatchment exchange coefficient [mm/day]	−5–3
X3	Routing store capacity (mm)	20–300
X4	Time constant of unit hydrograph (day)	1.1–2.9
CemaNeige
X5-Kf	Weighting coefficient of the snowpack thermal state [-]	0–1
X6-Ctg	Day-degree rate of melting [mm °C−1 day−1]	0–10
X7-Tacc	Accumulation threshold [mm]	0–45
X8-Tmelt	Fraction of annual snowfall defining the melt threshold [-]	0–1

. The ranges of the GR4J parameters correspond to the 80% confidence intervals obtained by Perrin et al. [56] calibrating the model over 429 catchments worldwide. The ranges of the CemaNeige model parameters are those suggested by Ayzel et al. [91] for Kf and Ctg. For the accumulation threshold (X7–Tacc), the mean annual snowfall of the lower elevation layer was considered to be the maximum value of the range (minimum zero). For the melting threshold (X8–Tmelt), the maximum flexibility was given to the model. The model was run with each parameter set for both 2000–2008 and 2009–2017 periods. Only parameter sets returning, for both periods, a Kling-Gupta Efficiency (KGE, [92,93]) larger than 0.7 and a Nash-Sutcliffe Efficiency [70] calculated on logarithmic values (lNSE) larger than 0.5 were retained for future simulations. KGE was selected as one of the two performance indices since it is rather complete, addressing linear correlation (r), bias (β), and variability (γ) occurring between simulated and observed discharge values. Similarly, lNSE was selected because it is the most common evaluation index for hydrologic models and in its logarithmic formulation it gives particular importance to low flow values, which are of particular interest in this study for the exploitation of the spring for drinking supply. The formulation of KGE is presented below, while for lNSE the reader can refer to Equation (1).

$$\mathrm{KGE}=1-\sqrt{{\left(\mathrm{r}-1\right)}^2+{\left(\mathsf{\beta }-1\right)}^2+{\left(\mathsf{\gamma }-1\right)}^2},$$

(6)

with:

$$\mathrm{r}=\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{N}}{\mathrm{OBS}}_{\mathrm{i}}·{\mathrm{MOD}}_{\mathrm{i}}-\frac{1}{\mathrm{N}}·\left({{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{N}}{\mathrm{OBS}}_{\mathrm{i}}\right)\left({{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{N}}{\mathrm{MOD}}_{\mathrm{i}}\right)}{\sqrt{\left[{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{N}}{\mathrm{OBS}}_{\mathrm{i}}{}^2-\frac{1}{\mathrm{N}}·{\left({{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{N}}{\mathrm{OBS}}_{\mathrm{i}}\right)}^2\right]·\left[{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{N}}{\mathrm{MOD}}_{\mathrm{i}}{}^2-\frac{1}{\mathrm{N}}·{\left({{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{N}}{\mathrm{MOD}}_{\mathrm{i}}\right)}^2\right]}} ,$$

(7)

$$\mathsf{\beta }=\frac{{\mathsf{\mu }}_{\mathrm{MOD}}}{{\mathsf{\mu }}_{\mathrm{OBS}}},$$

(8)

$$\mathsf{\gamma }=\frac{{\mathrm{CV}}_{\mathrm{MOD}}}{{\mathrm{CV}}_{\mathrm{OBS}}},$$

(9)

where μ is the mean flow rate [m³ s⁻¹], and CV is the coefficient of variation [-].

In addition, the outputs of the model obtained with the calibrated and validated parameter sets were checked against two emergency discharge thresholds, counting the days, and the maximum number of consecutive days below those discharges. Modeled values were compared with observations.

4.4. Changes in Discharge Regimes at Nossana Karst Spring

To evaluate the changes in discharge regime at Nossana karst spring, the model was run in ensemble mode with all the calibrated and validated parameter sets, and the downscaled future time series of precipitation and temperature as input. The last two years of observations (2016–2017) were used as the spin-up period. Changes in discharge regimes were analyzed mainly in terms of low flow. Variations in the number of days and consecutive days below the warning thresholds were calculated. In addition, an analysis of the general sensitivity of the spring to climate change was carried out, analyzing variations in average discharges in relation to modifications in precipitation amounts, precipitation partitioning into rain and snow, and temperature.

5. Results

5.1. Climate Models Evaluation

The results of the RCMs evaluation are presented in

Table 4. Performance statistics for the considered RCMs calculated in comparison with observed precipitation, minimum temperature, and maximum temperature monthly climatology. NSE is the Nash –Sutcliffe Efficiency coefficient, MAE the Mean Absolute Error and %MAE the normalized Mean Absolute Error. For RCMs name attribution, the reader can refer to Table 1.

	Precipitation			Tmin			Tmax
Model	NSE	MAE (mm)	%MAE	NSE	MAE (°C)	%MAE	NSE	MAE (°C)	%MAE
Mod_1 RCP 2.6	0.79	12.35	10.71	0.96	1.10	17.69	0.97	1.18	7.09
Mod_1 RCP 4.5	0.66	14.97	12.97	0.97	1.00	17.19	0.97	1.08	6.48
Mod_1 RCP 8.5	0.76	12.71	11.20	0.97	1.00	15.96	0.98	0.90	5.43
Mod_2 RCP 2.6	0.31	19.81	17.16	0.97	0.90	14.20	0.97	0.99	5.95
Mod_2 RCP 4.5	0.52	17.44	15.12	0.97	1.00	16.13	0.96	1.27	7.65
Mod_2 RCP 8.5	0.31	22.90	19.85	0.97	0.90	15.42	0.96	1.28	7.66
Mod_3 RCP 2.6	0.23	20.53	17.79	0.97	1.00	15.96	0.96	1.19	7.14
Mod_3 RCP 4.5	0.37	17.17	14.88	0.98	0.80	12.79	0.97	0.99	5.95
Mod_3 RCP 8.5	0.52	17.78	15.41	0.98	0.80	12.69	0.97	0.98	5.88

. Regarding precipitation, Mod_1 (see Table 1) is the model that best responds under all three RCP scenarios (divergence after 2005). Its NSE values result higher than 0.65 and the normalized MAE values lower than 15%. The other two models show larger errors and a lower fit with the observed climatology but they are generally able to represent the precipitation annual cycle, too. In addition, all models can satisfactorily catch the minimum and maximum temperature monthly climatology. NSE values are all very high (>0.95), while normalized MAE values for mean temperature range between 4% and 7%. All models were considered suitable for statistical downscaling.

5.2. Statistical Downscaling—Local Climate Change Signature

The results obtained from the statistical downscaling of the RCMs are presented by RCP scenario in twenty-year periods, to be consistent with observations (1998–2017). For convenience, the twenty-year periods 2021–2040, 2041–2060, 2061–2080, and 2081–2100 are defined as p1, p2, p3, and p4, respectively.

Projected changes in mean annual precipitation, for all models and RCPs, fall within the range suggested by Vaze et al. [42] for the application of rainfall-runoff models outside their calibration-validation conditions (

Table 5. Relative changes (%) in mean annual precipitation, in comparison to observations (1998–2017), for the three RCMs (Mod_1, Mod_2, Mod_3, see Table 1) under different RCPs, for four future periods (p1, 2021–2040; p2, 2041–2060; p3, 2061–2080; p4, 2081–2100).

	RCP2.6			RCP4.5			RCP8.5
Period	Mod_1 [%]	Mod_2 [%]	Mod_3 [%]	Mod_1 [%]	Mod_2 [%]	Mod_3 [%]	Mod_1 [%]	Mod_2 [%]	Mod_3 [%]
p1	−3.8	2.7	2.4	1.6	9.4	7.1	1.7	6.3	−14.1
p2	6.6	−5.2	0.7	−18.5	3.4	6.8	15.1	5.9	−6.9
p3	−4.4	5.4	−13.4	−3.0	0.9	−0.7	2.5	11.0	−2.4
p4	−4.4	8.3	−2.0	2.2	9.7	−8.4	−10.4	4.8	−7.9

). The only exception is the decrease (−18.5%) projected by Mod_1 for p2 and RCP4.5. In terms of temperature, mean annual changes range from 0.7 °C (RCP2.6, p1) to 5.8 °C (RCP8.5, p4). The actual occurrence of these changes would lead to a modification in the climatic regime of the area, according to the Köppen and Geiger classification [68]. The shift would be from Cfb to Cfa (i.e., a temperate climate characterized by wet summer and temperatures of the warmer month always higher than 22 °C). Except that for RCP8.5 and p4, the different RCMs usually agree in terms of projected temperature, while they show some differences in precipitation.

For the low emission scenario (RCP 2.6) (Figure 4), the monthly precipitation trend usually remains unchanged. Exceptions are January after 2040 and the spring season, both characterized by a slight increase in precipitation (15-20 mm). Also, Mod_1 projects a large increase of precipitation in May, for p1 and p2 (>25 mm), and an even larger decrease of precipitation in September (>50 mm), starting from p3. Temperatures in p1, p2, and p3 are projected to remain constant in the summer period, while showing a slight increase (around 1.0 °C) in the winter months. A strong increase in temperature is projected by Mod_1 in p4 for all the months of the year, up to a maximum of about 3 °C.

For the mid emission scenario (RCP4.5, Figure 5), in comparison to observations, the models (ensemble mean) project a slight increase in precipitation in the spring (around 20 mm, early spring in p1, shifting to May moving to p4), a marked decrease in the summer (up to 40 mm in August and 30 mm in September during p3), and another increase in autumn with peaks in October (around 40 mm in p3). The periods p2 and p3 are the most affected. Temperature trends are generally preserved, with a slightly higher increase observed in the winter months rather than in the summer period for p1 and p2. In the other two periods, changes appear quite homogeneous throughout the year, reaching up to approximately 2.5 °C.

The precipitation regime of the high emission scenario (RCP8.5, Figure 6) shows large differences in comparison to observations and the other two scenarios, exhibiting an extreme autumn peak (up to 322 mm/month according to Mod_2 for November in p2), a precipitation increase in the winter months (particularly in January), and a marked decrease in the summer period (25 mm or more in July, August and September for all models in p4). Temperatures constantly increase from p1 to p4, when they show monthly values in the summer of up to 4.7 °C (ensemble mean) larger than in the control period.

5.3. Hydrologic Model Calibration and Validation

The calibration-validation procedure led to a recognition of ensembles of ten 6-parameter sets (GR4J and CemaNeige model) and eight 8-parameter sets (GR4J and CemaNeige with hysteresis) able to satisfy the established performance criteria. The ensemble of calibrated-validated sets for the 6-parameter model with the corresponding performance indices are presented in

Table 6. Calibrated and validated parameter sets for the GR4J-CemaNeige model. KGE is the Kling–Gupta Efficiency, lNSE is the logarithmic Nash—Sutcliffe Efficiency, Cal indicates the calibration period (2000–2008), and Val indicates the validation period (2009–2017). For model parameters (X1–X6) see Table 3.

Parameter Set	X1	X2	X3	X4	X5	X6	lNSE Cal	KGE Cal	lNSE Val	KGE Val
Set 1	1091.94	2.07	171.26	1.35	0.01	3.22	0.67	0.73	0.51	0.71
Set 2	1078.07	2.81	226.15	1.32	0.18	4.38	0.66	0.72	0.51	0.72
Set 3	1196.85	2.98	219.46	1.31	0.08	2.91	0.65	0.72	0.53	0.72
Set 4	1187.37	2.53	188.70	1.34	0.18	2.06	0.64	0.72	0.52	0.71
Set 5	1164.21	2.97	182.71	1.18	0.01	6.33	0.67	0.76	0.53	0.72
Set 6	1075.06	2.50	207.87	1.30	0.05	3.10	0.66	0.72	0.52	0.71
Set 7	1118.65	2.72	173.46	1.30	0.06	6.74	0.67	0.75	0.52	0.71
Set 8	1146.42	2.38	144.62	1.30	0.04	3.74	0.68	0.76	0.54	0.71
Set 9	1050.21	2.55	181.51	1.12	0.10	5.07	0.68	0.76	0.52	0.71
Set 10	1158.85	2.74	176.41	1.46	0.15	5.89	0.67	0.75	0.52	0.71

. Since no significant differences in the simulated discharges were found between the two snow modules (for conciseness the results obtained with hysteresis are not shown), the model with six parameters was preferred.

Simulated discharge time series for the calibration and validation periods are shown in Figure 7, together with the non-exceedance probabilities. From the hydrograph, simulated discharges seem to approximate low-flow observations well, while peak flows are generally overestimated. Low-flow non-exceedance probabilities are also somewhat overestimated, especially for the calibration period. This model behavior is underlined by the mean discharge values and number of days below the warning thresholds. During calibration, the simulated mean flow is on average 7.7% (9.2% for validation) higher than the observed flow (3.48 m³ s⁻¹ and 4.00 m³ s⁻¹ for calibration and validation, respectively). In terms of the number of days below the thresholds, the model ensemble underestimates observations in both calibration and validation period. For the calibration period, the days below the 1.32 m³ s⁻¹ (1.52 m³ s⁻¹) threshold are 758 (992), while for the model ensemble they range from 141 to 310 (307-513). Similarly, for the validation period there are 335 (522) observed days below the 1.32 m³ s⁻¹ (1.52 m³ s⁻¹) threshold, while the model ensemble simulates a number of days between 122 and 207 (216–391). The model ensemble returns better results in terms of the maximum number of consecutive days below the thresholds, which are slightly underestimated with the exception of the lower threshold during the calibration period. Specifically, for the calibration period and the 1.32 m³ s⁻¹ (1.52 m³ s⁻¹) threshold, this index is calculated as 10 (46) days for observations and 58–63 (35–44) days for the model ensemble. The consecutive number of days below the threshold can therefore be considered the most reliable parameter to evaluate the future behavior of the Nossana Spring and discuss its future management. It is an indicator of a critical period length, during which the service company could need additional water resources to meet the demand.

5.4. Changes in Discharge Regimes at Nossana Karst Spring

The distribution of simulated mean discharge rates for the rainfall-runoff model ensemble, presented by RCP scenario, RCM, and twenty-year period, is shown through the boxplots in Figure 8. Simulated thirty-member (RCM and rainfall-runoff) ensemble mean discharges are lower than the mean observed flow (3.77 m³ s⁻¹) for every twenty-year period and RCP scenario, ranging from −3% (for p1 and RCP4.5) to −23% (for p4 and RCP8.5). However, a monotonous trend with time cannot be recognized. Also, the high emission scenario (RCP8.5) results in noticeably lower mean discharges than the other two RCPs only for p1 and p4.

Considering single RCPs, it is noted that the variability within the rainfall-runoff ensemble appears to be limited, while the variability between RCMs for the same period of time is much larger. Indeed, mean discharge values are strictly connected with mean annual rainfall amounts, which are extremely variable (see Table 5). For example, the mean annual rainfall decrease of 18.5% projected by Mod_2 for RCP4.5 and p2 corresponds to a 39% decrease in mean discharge. Similarly, the mean annual rainfall increase of 15.1% projected by Mod_1 for RCP8.5 and p2 corresponds to a 6% increase in mean discharge. However, although a decrease in mean annual rainfall is always closely related to a decrease in mean annual discharge, an increase in rainfall does not necessary lead to an increase in discharge. This lacking increase is related to a simultaneous increase in temperature and therefore evapotranspiration rates. As an example, for Mod_2 and RCP8.5, rainfall increases of about 6% and 5% (in p2 and p4) correspond to decreases in mean discharge of about 5% and 9%, respectively.

Figure 9 shows the relationships between the average snowpack thickness (expressed in snow water equivalent) and the 30-day average discharge in the Nossana basin, for the observation modelling period (2000–2017) and the four future periods under RCP 4.5. Under RCP 2.6 and RCP 8.5, the variables behave according to the same principles and thus they are not shown. First, Figure 9 shows that the model can reproduce the observed annual cycle of the 30-day discharge, with a slight underestimation of the peak occurring in May. The snowpack thickness is an output of the CemaNeige routine; no observed data are available, and so therefore it is possible that the snowpack thickness and the spring snowmelt are slightly underestimated. However, considering that April and May observed average precipitation is similar (Figure 5), the increase in the 30-day average discharge occurring during these two months is due to snowmelt. This dynamic is well represented by the model outputs. In future periods, the average thickness of the snowpack decreases with time, with similar relative paths in all months. Total precipitation in January, February, and March is not projected to decrease much (Figure 5), meaning that the relative future partitioning into rainfall and snow will favor the first. This reflects in 30-day discharges generally higher in all future intervals than during the observation period. Conversely, the 30-day discharge peak that used to occur in late May will be less pronounced and anticipated. In summer, the projected decreases in 30-day discharge are associated to both increasing temperature and decreasing precipitation (Figure 5). The 30-day discharge peaks projected for the autumn period (clearly evident in p1, p2, and p4) are mainly linked to increasing precipitation and only secondarily to increasing temperature, which affects precipitation partitioning. RCP 8.5 shows a similar annual trend in comparison to RCP 4.5 but average 30-day discharge during late spring and summer (May-September) are lower, up to an average of around 1.2 m³ s⁻¹ (Mod_1 p4). In addition, in p4 the spring peak for Mod_1 disappears and similar average discharges occur from late January to mid-June. Conversely, the major difference between RCP 4.5 and RCP 2.6 is the absence, for the latter, of the discharge peak in autumn (related to an absence of precipitation peak, see Figure 4).

The boxplots in Figure 10 (Figure 11) show the distribution of the number of days and maximum number of consecutive days below the 1.32 m³ s⁻¹ (1.52 m³ s⁻¹) warning threshold for the rainfall-runoff model ensemble, presented by RCP scenario, RCM, and the twenty-year period. Variability within the rainfall-runoff ensemble appears larger for these indices rather than for the mean discharge. Especially for the consecutive number of days below the warning thresholds and dry period, single parameterizations lead to values with differences larger than 70% between each other.

In terms of ensemble mean, the number of days below the warning thresholds (both values) are projected to decrease in all cases, except for p4 under RCP8.5. Looking at single RCMs, there are few cases of simulations counting an increasing number of days below the warning thresholds, in comparison to observations, mainly related to Mod_1 and Mod_3 under RCP4.5 and RCP8.5. Under RCP2.6, only Mod_3 projects an increasing number of days below thresholds limited to p3.

Conversely, in terms of ensemble mean, the maximum number of consecutive days below the thresholds is projected to increase after 2060 for all emission scenarios (except p3, RCP2.6), although not all RCMs agree. The length of the longest period below the 1.32 m³ s⁻¹ (1.52 m³ s⁻¹) warning threshold can be as long as 36 (64) extra days. This trend is particularly evident for the 1.52 m³ s⁻¹ warning threshold (Figure 11). However, this trend is particularly alarming for p4 since temperature is projected to increase from 1.0 °C to 5.8 °C (all RCPs and RCMs considered singularly), which will probably cause an increase in the water demand.

6. Discussion

Groppelli et al. [94,95] produced five future climate scenarios by statistical downscaling of three GCMs (2010–2060), forced under the A2 IPCC emission scenario [96], over the Oglio basin, the main valley east of the Seriana Valley. The A2 emission scenario, derived in the framework of 4th IPCC report [96], is somewhat comparable to RCP8.5, derived in the framework of the 5th IPCC report [2]. For the ten years around 2050 (in comparison to 1990–1999), they found an increase in mean annual temperature at 2000 m a.s.l. in the range 1.6–4.8 °C. At a lower elevation (around 600 m a.s.l.), for a comparable period (2041–2060) and scenario (RCP8.5), in this study a temperature increase in the range 1.8–2.8 °C was found. Larger increases could be expected only toward the end of the century (up to 5.8 °C). In terms of mean annual precipitation, Groppelli et al. [94,95] found variations between −15% and 40%, while in this study results suggest more limited changes (from −7% to 15%). The 18.5% decrease was obtained for the period 2041–2060, but under RCP4.5.

It is recognized that the selection of a rainfall-runoff modeling approach implies the assumption that the karst system will not undergo profound changes in the next 80 years, as pointed out by Hartmann et al. [16]. Specifically on the utilized rainfall-runoff model (GR4J with CemaNeige as snow accounting routing), it was found that the CemaNeige snow routine including hysteresis did not increase the performance of the model in terms of simulated discharge. This is in agreement with Riboust et al. [80], who found that accounting for hysteresis improved model simulations limited to snow cover area.

Gattinoni and Francani [60] calibrated a distributed equivalent porous media model (MODFLOW) for the simulation of Nossana Spring discharges and applied it to build spring depletion curves under variable recharge. The calibrated model and the depletion curves were then used to assess changes in spring behavior according to expected regional changes in temperature and precipitation, for the period 2080–2099 in comparison to 1980–1999, according to GCM forced under the A1B emission scenario [96]. All the models considered in their study suggested a decrease in mean annual precipitation (from 4% to 27%) and an increase in temperature (from 2.2 to 5.1 °C), which is in contrast with the results of this study obtained for the comparable RCP4.5 scenario (precipitation changes between −8% and 10%; temperature changes between 2.0 and 2.6 °C). However, due to the maximum expected precipitation decrease and evapotranspiration increase, they suggested a possible decrease in the spring discharge up to 40%, which is similar to the 39% decrease calculated in this study for the period 2041–2060 (Mod_2, see Figure 8b) for a precipitation decrease of 18.5% (Table 5) and a temperature increase of 2.1 °C. Gattinoni and Francani [60] also suggested that changes of the precipitation regime during the recharge season (March to November) are expected to have the largest influence on spring discharge. This is proven by the results of this study, which showed how the changes in precipitation partitioning between January and May will change the 30-day average discharge, anticipating and decreasing the observed late spring high flows. In addition, the projected rainfall decreases affecting the summer period (June to September) will reduce Nossana Spring discharges during this season. Higher recharge, in comparison to the recent past, can be expected in October, November, and December.

7. Conclusions

Climate change directly impacts the water cycle and consequently groundwater resources. Therefore, planning their management is essential for their protection and responsible use. In this study, a comprehensive methodological approach for the evaluation of future climate and discharge variations (up to 2100, in four twenty-year intervals) at Nossana karst spring was presented. The method is based on the statistical downscaling of three bias-corrected EURO-CORDEX RCMs forced by three different emission scenarios (RCP2.6, RCP4.5, RCP8.5) and the implementation of a rainfall-runoff model (GR4J) extended with a snow accounting routine (CemaNeige). Today, Nossana Spring is exploited by UniAcque S.p.A. for the supply of drinking water to 300,000 people (including the town of Bergamo). Simulated discharges were evaluated in terms of mean flow and in comparison to actual water demands (expressed as discharge warning thresholds), to provide projections that can be operationally useful to the service company.

Based on the presented results, the following can be concluded:

• The considered bias-corrected EURO-CORDEX RCMs have very good skills in reproducing observed temperature climatology (NSE > 0.95 and relative MAE < 10%) over the study area, while larger errors persist regarding precipitation (NSE between 0.25 and 0.65, relative MAE between 10% and 20%);
• According to the downscaled RCMs data, in comparison to 1998–2017, mean temperature will likely increase throughout the rest of the XXI century, from 0.7 °C in 2021–2040 (RCP4.5, Mod_2) to 5.8 °C in 2081–2100 (RCP8.5, Mod_1);
• Downscaled RCMs data do not show a clear trend in precipitation. For all twenty-year periods and RCP scenarios, there are single RCMs projecting increasing and decreasing rainfall (except 2021–2040, RCP2.6, all increasing). Variations in mean annual rainfall varies between −18.5% (2041–2060, RCP4.5, Mod_2) and 15.1% (2041–2060, RCP8.5, Mod_2);
• A pronounced decrease of precipitation is expected in the summer period after 2060, as most RCM-RCP combinations show;
• Mean discharges are generally projected to decrease in comparison to observed flow (3.77 m³ s⁻¹) since changes in mean annual precipitation usually do not balance increases in evapotranspiration rates due to higher temperatures;
• Variability in the projected mean discharges is mainly linked to the meteorological input rather than the rainfall-runoff model parameterization;
• The maximum number of consecutive days below the warning thresholds was recognized as the best index to evaluate the spring low flow conditions;
• After 2060, the length of the periods with discharge lower than the warning thresholds is expected to increase. These periods could last up to 64 days (86%) longer than in 1998–2017.

The present study indicates that additional water resources might be needed to satisfy the population water demand in the Nossana Spring area, especially after 2060. Therefore, the presented approach can be considered useful to provide indications regarding the management of the resource in the mid- long-term period. In fact, UniAcque S.p.A. will soon start to investigate possible alternative sources for water supply. However, it is evident that the results include large uncertainties derived from the climate input variability, especially precipitation, and as new emission scenarios are produced, the study should be updated.

Furthermore, the study focused on the quantitative aspects of the spring behavior; however, it is possible that a change in the discharge regime could also influence the chemical-physical characteristics of the water and its quality. A chemical-isotopic monitoring of the spring could improve the conceptual model of the aquifer system and the knowledge of the system response. This would allow a refinement of the method, raising it to a qualitative and quantitative prevision of the resource behavior.