Aquifer recharge in the Piedmont Alpine zone: Historical trends and future scenarios

and future scenarios Elisa Brussolo1, Elisa Palazzi2, Jost von Hardenberg3,2, Giulio Masetti4, Gianna Vivaldo4, Maurizio Previati5, Davide Canone5, Davide Gisolo5, Ivan Bevilacqua5, Antonello Provenzale4, and Stefano Ferraris5 1Research Center, Società Metropolitana Acque Torino S.p.A., Torino, Italy 2Institute of Atmospheric Sciences and Climate, National Research Council of Italy (CNR), Torino, Italy 3Department of Environment, Land and Infrastructure Engineering (DIATI), Politecnico di Torino, Torino, Italy 4Institute of Geosciences and Earth Resources, National Research Council of Italy (CNR), Pisa, Italy 5Interuniversity Department of Regional and Urban Studies and Planning (DIST), Politecnico di Torino and Università di Torino, Torino, Italy Correspondence: Elisa Palazzi (e.palazzi@isac.cnr.it)

Owing to its hydrogeological features (see also the rivers system shown in Fig. 1), the foothill aquifer systems generally 125 sustain the infiltration of both local rainfall and streamwater originating in mountain catchments. These systems are highly sensitive to variation and changes in meteoclimatic variables, first of all in precipitation regimes.

Water balance terms
In this paper we refer to water balance as the balance between water inputs and outputs at the catchment scale, as schematically shown in Figure 2. 130 Being precipitation much more episodic than evapotranspiration, daily time steps are recommended with respect to longer ones, in order to avoid the underestimation of drainage and recharge (being the distinction between them one of lag time (e.g. Healy, 2010). Therefore in this work the computation is done at the daily scale in the soil column, and results are then aggregated at the yearly and quarterly time scale and, spatially, at the catchment scale. The daily soil model will be described later in the text. The yearly catchment water balance can be simplified as follows 135 (Healy, 2010): where P q represents the sum of liquid precipitation (rainfall) and snowmelt, AET the Actual Evapotranspiration and Q out the drainage. Input surface or subsurface flow can be neglected because all catchments are bordered by the mountain divide. Storage To account for the different ground characteristics and variations, together with the physical description of the hydrological processes, all the water balance variables (except for output discharge) must be first evaluated at a fine spatial resolution and then aggregated (i.e., upscaled) at watershed level. To this end, a horizontal spatial resolution of 250 m was considered as a good compromise between the necessity of high-resolution and the computational resources required: the study area was 145 discretised into a grid of 652 x 521 pixels, with a 250 m × 250 m spatial resolution.
The model takes into account the effects of irrigation on the value of AET . However, input and output irrigation terms in surface and groundwater balance are assumed to balance out at the catchment level. Industrial water withdrawals are mainly for hydroelectric production and therefore they are not changing the water balance at the watershed scale, because they give back the captured water at short distance.

150
The contribution of glacier and permanent snow melting is disregardable with respect to the catchment scale. However the evaluation of snowmelt from the snowpack is quite important as it represents the amount of daily snow equivalent contributing to the water balance. To this end, as a first step, air temperature was used to distinguish between snowfall (P s ) and rainfall (P r ) (DeWalle and A.Rango, 2008): for temperatures lower or equal to -2 • C all the precipitation amount can be considered as snowfall, while for temperature higher or equal to 5 • C all the precipitation amount can be considered as rainfall. For 155 temperatures between -2 • C and 5 • C, we linearly interpolated rainfall and snowfall amounts to breakdown the different percentages of precipitation type. The resulting snowfall P s is used as an input to a a simple bucket model in each pixel, used to integrate the storage of snow on ground, S. The daily snowmelt, N f , is subtracted from the same bucket and it is evaluated using the method by Zeinivand and Smedt (2009), where the melted snow is a function of the difference between the daily mean temperature T mean and melting temperature T 0 (0 • C) and of the rainfall amount (P r evaluated in mm/day), as follows: were k rain (rainfall melt-rate factor) and k snow (melt-rate factor) are fixed parameters. We use the same values as in Zeinivand and Smedt (2009): k rain = 0.0757 • C −1 and k snow =3 mm day −1 • C −1 .
A correct evaluation of the water balance allows a reliable estimation of water availability. To this end we applied the following chain in the study area: 1. identification and retrieval of the daily meteoclimatic (i.e. temperature and precipitation) data; 165 2. regridding of the meteoclimatic and hydrological data at the spatial resolution required by the hydrological model, namely a square 250 m grid; 3. setting-up of a mathematical model that accounts for soil-vegetation-atmosphere interactions and providing actual evapotranspiration estimates (see Section 4).
Observational and reanalysis datasets have been used to obtain meteoclimatic data for the period 1959-2017, as described 170 in Section 3, while for the assessment of future water balance scenarios we relied on climate projections of temperature and precipitation (see Section 5). Table 1. Average climatological monthly adiabatic lapse-rate from Rolland (2003) for Alpine region (expressed in K/km) Jan Feb Mar Apr May June July Aug Sep Oct Nov Dec 4.5 5.0 5.8 6.2 6.5 6.5 6.5 6.5 6.0 5.5 5.0 4.5   (Allen et al., 1998). A simplified scheme is shown in Figure 4, which includes only one layer of the bucket model. The water inputs for the model are precipitation (sum of rainfall and melted snow) and irrigation, if any. The outputs are AET and drainage (sum of runoff and deep percolation). The distinction between runoff and deep percolation is not included in the model, as well as the modelling of capillary rise.

215
For each layer, the bucket equation is as follows: where P q is rainfall plus snowmelt; I * is irrigation; Q out is drainage; AET is actual evapotranspiration; ∆S is the variation of water storage in the soil per day. In a bucket model (the one by Baudena et al. (2012) was applied in the same area in northwestern Italy) water can flow upwards within each layer via evapotranspiration or evaporation, and downwards via deep percolation, as a function of soil water content of each layer: if it exceeds field capacity, the water surplus percolates into the 220 layer immediately below.

Input data of the soil water model
The soil hydraulic properties have been estimated via pedotransfer (PTF) following Schaap et al. (2001) where r is the fraction of volume occupied by stones, and L d is the layer depth, which is obtained by dividing by seven the root zone depth z. The root zone depth has been obtained using the land cover classes from the BDTRE database (Regione Piemonte, 2018), listed in Table 2, converted into root zone depth using different root zone depths for each class, as listed in Table 3. The resulting root zone depth was compared with that also provided by the Regione Piemonte Soil Map (IPLA, 2007), choosing the minimum between the two. For trees the water from the whole soil depth can be depleted and the IPLA 230 depths have been used. Because of the high spatial resolution of BDTRE, the raster has been produced with a 5 m × 5 m grid, Plain meadows (< 1000m asl) 3 Orchards 4 Horticulturals 5 Plain broadleaves (< 800m asl) 6 Mountain broadleaves (> 800 m asl) 7 Coniferous 8 Vineyards 9 Mountain grassland (> 1000m asl) 10 Bare soil 11 Bare rocks 12

Impervious surfaces 13
Water 14 Others 15 Glaciers and permanent snow 16 then resampled at 250 m by the mode criterium. Homogenising data as a function of the evapotranspiration behavior of each class, the different land cover classes have been aggregated into 16, as reported in Table 2.
Finally, considering that irrigation (here considered as water input) significantly modifies actual evapotranspiration, for quantifying the actually irrigated fields, we use the Regional Irrigation Information system (Regione Piemonte, 2016), the agri- collected in experiments performed in real world farms both for surface irrigation (Canone et al., 2015) and for sprinklers (Canone et al., 2016).

Actual Evapotranspiration computation 240
The Actual Evapotranspiration AET was calculated starting from potential evapotranspiration P ET , then reduced by considering the actual soil water content, obtained from the bucket soil model for each layer, with AET of each day being the sum of the depletions from the single layers. P ET was calculated with the model of Hargreaves and Samani, as in Allen et al. (1998), using daily maximum and minimum air temperature data from the OI dataset and extraterrestrial radiation modeled as in Aguilar et al. (2010). Also, the reduction of evapotranspiration (from PET to AET) due to actual soil water content was 245 modeled by using the coefficients Kc and Ks related to crop and soil respectively according to Allen et al. (1998).

Future projections of precipitation and temperature
Reliable estimates of the hydrological response and of water availability in the coming decades generally require the implementation of a modelling chain consisting of global climate models (GCMs), which provide climate scenarios for the entire planet, regional climate models (RCMs) nested into global models providing lateral and boundary conditions for the regional 250 simulation and, depending on the resolution which needs to be achieved, further downscaling procedures. At the end of the modelling chain, a hydrological model is thus forced with a high-resolution climatic input to simulate the hydrological response at the scale of interest. In this study, we used a small multi-member ensemble of the RCA4 regional climate model forced by five GCM simulations, able to provide the climatic variables of interest -minimum and maximum temperature and precipitation -at a spatial resolution of ∼12.5 km. More detailed analyses of the interplay between several GCMs and RCMs 255 (Sorland et al., 2018) are out of the operational water management scope of this paper. In this work we proceeded similarly to another recharge study (Allen et al., 2010)

RCA4 Regional Climate Model
We analyzed the available simulations of the RCA4 RCM (Strandberg et al., 2014) driven by five different GCMs, namely EC-

265
Earth, CNRM-CM5, IPSL-CM5A-MR, HadGEM2-ES, MPI-ESM-LR, which provide lateral and boundary conditions for the regional simulation. Using one single RCM allowed us to obtain an ensemble of reasonably homogeneous simulations at the regional level, but representing at the same time model uncertainties in future projections captured by the different large-scale GCMs. RCA4 is one state-of-the-art RCM which participated to CORDEX, the Coordinated Regional Climate Downscaling Experiment (http://wcrp-cordex.ipsl.jussieu.fr/, Giorgi et al., 2009) Strandberg et al. (2014). at about twice the pre-industrial level by the end of the century. RCP8.5 is an extreme business-as-usual scenario in which greenhouse gas emissions are not expected to stabilise and carbon dioxide concentrations will be more than tripled at the end of the century compared to preindustrial levels. properties of the climatic input. Bias correction methods are usually applied to correct the differences between climate model output and observed climatologies and are different depending on the variable and specific application which is considered 300 (Hempel et al., 2013;Maraun, 2013;Maraun et al., 2010). In this study we used an additive correction factor for adjusting the temperature and a multiplicative correction factor for precipitation (a standard procedure for positive-defined fields) applied pixel by pixel and constant in time, in order to correct the differences in the long-term climatology calculated over a common time period between the simulated and observed fields. To this end, we calculated the long term mean of the historical Euro-CORDEX simulations and of the OI dataset, already interpolated on the UTM grid, in the 30-year long period from 1986 to 305 2015. In order to maintain the physical consistency between the minimum temperature and the maximum temperature, and thus avoid possible inversions, the same correction factor was used for daily minimum and maximum temperature data, calculated from the average between the daily minimum temperature and the maximum temperatures. The correction factor calculated to correct the model bias in the historical reference period is then applied to the future simulations. In addition to the bias adjustment, temperatures have been further corrected for the lapse rate, as already done for the OI data, based on Rolland 310 (2003).

Historical trends
Historical data from 1959 to 2017 were analyzed in each river catchment and sub-catchment, aggregating all data previously evaluated at the spatial resolution of 250 m, to find out significant trends both of the water balance terms and the meteorological 315 variables, for both the hydrological year and the quarterly analysis.
A linear regression was performed to calculate the temporal trends. The goodness of line fit was estimated with the coefficient of determination R 2 and its statistical significance was evaluated considering the p − value at 5% level (i.e., 95% significance) (Wilks, 2011).

Hydrological year analysis 320
During the hydrological year, a complete snow melting occurs in the considered area. Thus, in Eq. 1 P q refers to the whole precipitation amount and P q − AET to drainage. Figure  Figure 5 reveals that, for this catchment, the daily maximum temperature has increased during the study period, with a statistically significant trend. The total precipitation, P q , has decreased (p=0.059) while AET has increased (p=0.054), and drainage has significantly decreased (p=0.032). The catchment considered here is in the driest part of the study region, which is also the part where most of the significant trends of precipitation and AET were detected in the data analysis. The same 330 approach is adopted for all the catchments in the study area. The results are summarised in Table 4     While yearly precipitation trends are overall negative, the quarterly analysis shows that a significant decrease of precipitation in the first and second quarters occurred from 1959 to 2017, confirming that precipitation trends in NW Italy depend on the considered temporal aggregation (Brunetti et al., 2006). The drainage trend analysis shows an overall reduction in the first and   Table 5 shows the drainage trend slopes from 2018 up to 2050 for each catchment in the study area using climate projections under the RCP 4.5 and the RCP 8.5 scenarios. Bold font cells indicate trend values that are statistically significant. Drainage trends in the RCP 4.5 scenario are often negative but not statistically significant. These negative trends become stronger and occasionally statistically significant in the more extreme RCP 8.5 scenario.

390
As a general statement, taking into account also the variability within the projection ensemble, the rather small values of drainage trends and their limited significance suggest that the Greater Turin metropolitan area will probably experience a yearly steady-state drainage situation till 2050, with total precipitation that shows a slight decrease. We can also notice a slight drainage increase with higher precipitation amounts, a decreasing trend with higher daily maximum temperature and, above all, a strong interannual variability (see, as an example, the standard deviation evaluated for the detrended timeseries shown in 395 Figure 8).
The projections of P q − AET (drainage) for two different cross-sections of the study area, the Dora Riparia station in Turin (catchment ID 66) and the Orco station in S. Benigno (catchment ID 109), are shown in Figure 8. Each row corresponds to each of the five GCMs which drive the RCM; the results obtained with the RCP 4.5 (RCP 8.5) scenario are in red (green). The slope of the regression line, the coefficient of determination R 2 , the p-value and the standard deviation are also indicated. We 400 choose to show these two cases because the two catchments are among the largest ones (as from Table 4, 1243 km 2 and 829 km 2 , respectively) and they represent, respectively, the drier (western) and the wetter (northern) parts of the region.
The overall results at yearly time scale, not shown for the sake of brevity, can be summarized as follows: -T max shows positive trends, that are almost always significant for all model realizations and both scenarios. Also T min has similar trends;

405
-P q has either positive or negative trends, rarely significant; -AET has either positive and negative trends, rarely significant; slight positive drainage trends with higher precipitation amounts, negative trends with higher T max and high interannual variability.
A large variability between the different projections was obtained, as in Crosbie et al. (2013) and in Persaud et al. (2020), 410 with a much less clear pattern in spatial variability if compared with the past data trends evaluations.

Quarterly analysis
Projections of drainage evaluated at the seasonal timescale show again a large variability within the ensemble of projections.
However, a general tendency to drainage increase in the first quarter and to drainage decrease in the third and fourth quarters, JAS and ON D, emerges. Tables 6-9 show the drainage trend (expressed in mm/quarter/year) in the four Quarters. All models 415 and scenarios are indicated. Bold font cells indicate trend values that are statistically significant.
In the first quarter (JF M ) the models mostly agree to indicate a drainage increase till 2050 over the whole study area, with a wider agreement for RCP4.5 (4 out of 5 models concur in all catchments). In the second quarter (AM J) in all river basins at least 3 out of 5 models estimate a drainage decrease for the RCP 8.5 scenario, not for RCP4.5. In the third (3 out of 5 models) and fourth (4 out of 5 models) quarters there is an overall tendency of drainage decrease in the whole study area for RCP4.5.

420
For the other meteoclimatic and water balance variables, at quarterly timescale we can report that (we omit the details for brevity): -T max and T min show positive trends in all quarters, almost always significant for all models and scenarios; in all quarters both rainfall plus snowmelt P q (including snowmelt) and rainfall (P r ) have both positive and negative trends, rarely significant. Rainfall has almost always a weak positive trend in the first quarter, rarely significant. Snowfall 425 shows broad negative trends, almost always significant; -AET shows broad positive trends in the first and second quarters, despite some differences between the models. In the third quarter an overall decreasing tendency can be observed, while in the fourth quarter trend slopes have values close to zero. tion time variability plays the major role. Konapala et al. (2020) indicate changes in long-term drainage, but they recall the limitations in the ability of current generation coupled climate models to capture the key drivers of persistent weather extreme.

Conclusions
Assessing the impacts of climate change on groundwater resources represents a priority in water management, besides being an important scientific challenge. In this study we focused on an area where the aquifer providing water to about 2.3 million 435 people is catched in watersheds characterized by a very large spatial variability of precipitation, due to the proximity of high mountains and of the sea. In this work we assessed the trends of precipitation, temperature and actual evapotranspiration in order to estimate the trend of drainage as a proxy of the water available for drinking purposes. The analyses have been performed both at the hydrological year and at the quarterly timescales. We analysed past and future conditions, using a historical climate dataset providing temperature and precipitation data for the period 1959-2017 and future projections from 440 a multi-member ensemble of a regional climate model up to 2050, in order to be compliant with the time frame of water management objectives. Our analysis revealed a very strong interannual variability of drainage in the historical period, as well as remarkable spatial differences, i.e., across the different watersheds. We found that the driest part of the region (the Dora Riparia Valley) shows significant negative drainage trends both yearly and in the spring quarter. The increasing trend of actual evapotranspiration is typical of colder and moister areas, like England (Blyth et al., 2018). The combination of it with a very 445 variable precipitation input has to be monitored in the future for its potential effects on drainage. As for the future projections, there is a large inter-model variability, as in Crosbie et al. (2013) and in Persaud et al. (2020), of all hydrological variables, which is evident for all catchments, and almost no significant trends. This last result is in line with Gudmundsson et al. (2017) for this part of Europe. Runoff was not considered in this study, but following Kumar et al. (2016) it could be useful to quantify its variability with climate change.

450
This study constitutes a knowledge basis which helps for a better informed management, infrastructural and supply decisions for the study area considered, and our methodology could be extended also to other regions. This research is the result of a shared work developed together with a drinking water authority and academic/research institutions, in a general context which stimulated and supported a participatory planning that is scientifically grounded. Our results will be integrated in the definition and refinement of study-area long-term guidelines and strategic development for groundwater resources protection 455 and infrastructural provision and planning as a part of a wider effort to strengthen good water resource management and governance.
Data availability. The datasets presented in this study can be obtained upon request to the corresponding author.
Responsible Resilience Interdepartmental Centre (R3C) DIST-PoliTO for valuable collaboration and ARPA Piemonte for kind support. The observational datasets used in this study are all freely available. Post-processed data are reproducible following the detailed descriptions of Section 2 and 4 and available upon request to the authors.