Introduction

Recharge of mountainous aquifers in Mediterranean areas, and in the rest of the world, is dependent on precipitation amount and distribution (Viviroli et al. 2007; Penna et al. 2016; Schreiner-McGraw and Ajami 2022), with effects also related to the elevation of the ridges constituting the recharge areas. Infiltration towards the water table is strongly influenced by the persistence of snowpack and related temperature changes (Wu and Xu 2005; Luce et al. 2014; Sospedra-Alfonso et al. 2015; IPCC 2019; Lucianetti et al. 2020; Somers and McKenzie 2020). In a framework of increasing awareness of climate-change effects, such as infiltration reduction and aquifer resource depletion (the latter magnified by human consumption), it is important to update groundwater recharge assessment methods based on actual data in order to guarantee sustainable management of groundwater resources, as required also to attain achievement of the United Nation’s Sustainable Development Goals (SDGs; IAH 2017; Ferguson et al. 2020; Gleeson et al. 2020; Wu et al. 2020; IPCC 2021).

Fractured carbonate aquifers are frequently characterized by a fragile equilibrium, because of their relationship with karst processes (e.g. dissolution, surface morphology, underground conduits development, dual-flow aquifer response), which modulate the input/output of groundwater flow. Thus, recharge of such aquifers is strongly correlated with topographic effects, and these become an essential piece of basic information for different management purposes, for example establishing protection zones for drinking water supply, managing climate-change effects such as water shortage during drought periods, and more generally planning a rational use of groundwater resources (Hartmann et al. 2014; Fiorillo et al. 2015; Vigna and Banzato 2015; Goldscheider et al. 2020; Iacurto et al. 2020).

A very well-known tool to identify the mean elevation of recharge areas is based on the properties of water stable isotopes; these remain unchanged after water infiltration and during groundwater flow in the saturated zone in a classical recharge/discharge cycle not affected by deep fluid contribution (Kendall and McDonnell 1998; Liebminger et al. 2007; Carucci et al. 2012; Beria et al. 2018; Sbarbati et al. 2021). Consequently, the use of stable isotopes for calculating the mean recharge elevation is still a basic tool applicable to small, medium and large hydrogeological basins, with particular reference to mountainous aquifers where recharge areas are distributed in a wide elevation range (e.g. up to thousands of meters). The mean elevation of the recharge area calculated by isotope values of spring water is also representative in large basins, because longer flowpaths facilitate the mixing of recharge contributions, leading to steady average values of stable isotopes (Jasechko et al. 2017; Jódar et al. 2021).

In hydrogeology, the computed isotope recharge elevation (CIRE; Petitta et al. 2010) is a common concept to express the mean elevation of the recharge area by measuring the isotope composition in discharge areas (from springs and wells), to be compared with rainfall isotope values recorded at different elevations in the related recharge area of the same aquifer feeding the springs/wells (Tallini et al. 2014; Meng et al. 2015; Mussi et al. 2017; Sappa et al. 2018; Matheswaran et al. 2019; Tazioli et al. 2019; Allerano et al. 2020; Lucianetti et al. 2020; Lugodisha et al. 2020; Pfeiffer et al. 2021). Of course, this approach has limitations, whereby the major constraint is the knowledge of an “average” isotopic value of the discharge water, which has been mixed along the flowpath, collecting different inputs from different recharge zones and recharge episodes. Consequently, results are also time-dependent, because of the likely change in isotope composition not only with the season but also during a single precipitation event (Longinelli and Selmo 2003; Minissale and Vaselli 2011; Giustini et al. 2016). Nevertheless, in fractured carbonate aquifers where fast flow in karst conduits is limited, as happens in central Italy, groundwater flow is mainly governed by the fracture network, leading to steady discharge at the base-flow springs; in these cases, stable isotopes show steady values (Petitta et al. 2015, 2018; Mastrorillo et al. 2020), limiting the seasonal effect of recharge-water isotope composition on the groundwater discharge.

In this framework, the calculation of CIRE after measuring rainfall values and building elevation/stable isotopes relationships for recharge is possible and related correlation lines are available in the literature at the national, regional and sometimes local (single aquifer) scale (Desiderio et al. 2005; Barbieri et al. 2005; Tarragoni 2006; Tazioli et al. 2007; Spadoni et al. 2009; Petitta et al. 2010; Minissale and Vaselli 2011; Nanni et al. 2013). In terms of hydrogeological budget, knowing the extent of the recharge area and the total amount of spring discharge, a simple equation allows one to obtain the infiltration rate at a yearly scale for hydrogeological basins (Boni et al. 1986). The subsequent challenge is to assess the distribution of this effective infiltration with elevation, starting from isotopic values useful only for an “average” calculation of the CIRE.

In fact, especially in mountainous carbonate fractured aquifers, recharge areas close to the discharge point are delivering fast groundwater flow, compared with high-elevation recharge zones whose contribution has usually longer residence time. Spring discharge hydrographs are able to reveal response time to infiltration, but chemiographs and tracer tests are necessary to assess groundwater residence time, which frequently does not correspond to the response time. The spatial heterogeneity and the temporal distribution of recharge rate on these aquifers not only depend on the tectono-stratigraphic setting and on water table depth, but surficial karst features (e.g. endorheic basins, fracture density) and snow cover persistence over time also have a significant role. In fact, recent studies suggest that snowpack produces a higher rate of recharge (especially in fractured aquifers) with respect to an equivalent amount of rain (Wu and Xu 2005; Kløve et al. 2017; Cummings and Eibert 2018; Nygren et al. 2020).

Consequently, it would be generally inferred that the average infiltration rate has a different contribution, in terms of rainfall amount, with elevation (higher areas offer higher percentage of recharge), as enhanced by snow cover presence and persistence with time (Beniston et al. 2003; Wu and Xu 2005; Nanni et al. 2013; Chen et al. 2018; Lucianetti et al. 2020; Notarnicola 2020; Di Marco et al. 2021). Conversely, the stable isotopes seem to offer the sole “average” value of the CIRE. Comparing the CIRE with the mean elevation of the recharge area in different carbonate aquifers, it becomes possible to verify the correspondence or the discrepancy between topographic mean elevation and CIRE. This approach independently highlights possible distributed and enriched contribution from higher-elevation zones of the recharge areas.

The aim of this work is: (1) to find a realistic distribution with the elevation of infiltration rates, and (2) to assess the influence of high-elevation snowpack as the preferential input for aquifer recharge in mountainous carbonate fractured aquifers. These challenges have been tackled by analyzing the correspondence and/or the discrepancy between the mean topographic elevation of the recharge areas and the related CIRE obtained by stable isotope data taken from aquifers of central Italy. A validation of the findings has been conducted by modern satellite images, to investigate apparent anomalies found during conducted recharge elevation comparisons. By this approach, it is here intended to offer not only a contribution to fill the existing gaps in aquifer recharge assessment, but mainly, to elucidate the potential effects of changes in precipitation nature (snow/rain) on future availability of groundwater resources.

Materials and Methods

Hydrogeological recharge area datasets

The investigated area is located in central Italy, along the Apennines carbonate ridge; in this context carbonate reliefs represent wide recharge areas feeding the main regional fractured aquifers (Petitta 2009). These hydrogeological systems are highly permeable thanks to fracturing and karst processes, and therefore they are characterized by high effective infiltration rates (reaching up to 800–1,000 mm/year, according to Boni et al. 1986). The outcropping of continental deposits at the boundaries of carbonate regional aquifers associated with recent tectonic activities have limited mature karst development in the saturated zone of the aquifer (Cavinato and De Celles 1999; Petitta 2009), affecting groundwater flow. In the studied carbonate aquifers, the Meso-Cenozoic sequence of the carbonate platform has a thickness greater than 2,000 m; this is because the upper Miocene in the Apennine chain was uplifted by a compressive tectonic phase, thrusting the sequence towards the E–NE, followed during the Quaternary by extensional tectonic movement, which has strongly contributed to the permeability due to fracturing. Recent karst development improves the rock permeability, with the main effects on the infiltration phase due to the development of surface karst morphology. Nevertheless, the continuous Quaternary shift of the permeability limits at the base of the aquifers by contact with detrital/alluvial deposits (due to changes of base level by erosion and eustatic changes) inhibits a mature karst development in the discharge areas. These conditions generate relevant springs (basal springs) with very high and steady discharges, only slightly influenced by seasonal variability (Petitta 2009). In this area of the Central Apennines, several stable isotope studies on the hydrogeological dynamics have been carried out: Desiderio et al. 2005 used stable isotopes to assess the origin of mineralized groundwater as well as to define isotopic mean recharge altitudes for basal springs; Barbieri et al. 2005 and Tarragoni 2006 validated the hydrogeological conceptual models of two main Apennines hydrostructures; Tazioli et al. 2007 focused their attention on both the elevation of the recharge area and recharge mechanisms; Spadoni et al. 2009 validated the digital recharge model through study of isotopes; Petitta et al. 2010 assessed the role of a rockslide dam body on the hydrodynamics of a system; Minissale and Vaselli 2011 treated karst springs as natural pluviometers defining rain isotope compositions; and, finally Nanni et al. 2013 highlighted the influence of snow melt on the isotopic signature of springs.

However, some uncertainties still arise when the stable isotopes were applied to particularly complex hydrogeological situations such as fractured aquifer ridges of the Central Apennines (Nanni et al. 2013). For this reason, to verify the accordance or the discrepancy between mean isotopic recharge areas and mean topographic recharge areas, 17 basal springs (Fig. 1) were selected and sampled in this work. All sampled springs are fed by regional aquifers, showing a deep regional flowpath in a carbonate fractured/karst environment.

Fig. 1
figure 1

Location map of selected springs and related regional aquifers selected for the study. Colours identified different hydrogeological basins. For acronyms refer to Table 1. Blue circles are springs

Full size image

The springs selected for this research show common hydrogeological and geological features as:

  1. 1.

    A very limited dependence of the stable isotopes 18O and 2H on seasonal variations, linked with a very steady discharge rate; this smooth isotopic response is reinforced by the groundwater flow being based on the fracture network and consequent relevant residence time of infiltrating waters in the aquifer. Consequently, the isotopic values of the basal springs correspond to a long-term average of groundwater flow, with limited changes over the years.

  2. 2.

    An enhanced waterflow mixing, and even a representative average isotopic value, thanks to the extent of the regional recharge areas. In fact, the larger the recharge area and the slower the flowpath from infiltration to the discharge zones, the more statistically uniform the mixing of infiltrated waters can be, including also the groundwater flow at the regional scale.

  3. 3.

    All the springs are part of well-known hydrogeological systems studied in recent times (Celico 1983; Conese et al. 2001; Petitta and Tallini 2002; Nanni and Rusi 2003; Boni and Ruisi 2006; Tarragoni 2006; Adinolfi Falcone et al. 2008; Agrillo et al. 2008; Petitta 2009; Spadoni et al. 2009; Boni et al. 2010; Allocca et al. 2014; Banzato 2014; Saroli et al. 2019). In fact, these groundwater systems and their hydrogeological boundaries are well defined by stratigraphic and tectonic boundaries with aquitards such as the Quaternary clastic continental deposits or the terrigenous deposits.

  4. 4.

    All the springs are representative of different regional aquifers, as summarized in Table 1.

Table 1 Data on the selected springs (spring numbers and aquifer abbreviations are shown in Fig. 1)
Full size table

Mean discharge values of springs (Table 1) refer to several previous studies, including some values that have not been published. The extent of each recharge area (Table 1) has been evaluated from existing studies, and the data are frequently derived from unpublished technical reports. The 17 selected springs, from north to south, are characterized by significant differences in the aforementioned parameters. Indeed, elevation varies from a minimum of 41 m above sea level (m a.s.l.; spring S17) to a maximum of 951 m a.s.l. (S2); mean discharge ranges between 0.125 m3/s (S17) and 18 m3/s (S6); the extension and mean elevation of recharge areas range respectively between 31 km2 (S11) and 1,016 km2 (S5 and S6), and between 650 m a.s.l. (S17) and 1,975 m a.s.l. (S11).

Nevertheless, the considered springs can be classified on the basis of their correspondence to an entire hydrogeological basin or to a portion of the aquifer system, as follows:

  • Basal springs fed by the entire hydrogeological system, where recharge area calculation has high reliability (as for springs S2, S3, S4, S15 and S16).

  • Basal springs representing the main exit of a wide regional aquifer, which is feeding also different springs; in this case, the uncertainty in drawing the recharge area becomes larger, but the representativeness of the calculated recharge areas remains very high (as for springs S1, S7, S9, S10, S12, S14).

  • Springs fed by a portion of the regional aquifer, feeding more than one spring by the same recharge area (e.g. S5, S6, S8, S11, S13 and S17), which implies a larger uncertainty in identifying the effective recharge area, based on subbasin internal groundwater divides, with possible limited reliability of recharge area boundaries.

δ18O and δD measurements

To obtain detailed information about stable isotope values, two sampling surveys were carried out in 2016 (June and October). The sampling was conducted in two different seasons to verify the influence of seasonal variation. Nevertheless, isotopes data collected in previous studies, about the same springs, were considered. This additional analysis has the aim to verify the relationship of the selected springs with the isotopic meteoric water lines and 18O and 2H groundwater lines, defined in the literature for each selected spring. All the samples were collected directly from each spring outlet and stored in cool boxes at 4 °C until reaching the isotopes laboratory. Isotopic composition is expressed in part per mil accordingly to the usual δ notation, with respect to the V-SMOW (Vienna Standard Mean Ocean Water). Analysis was performed in the Isotopes Laboratory of the University of Parma (Italy) using a mass spectrometer equilibrated on 25 °C as described by Epstein and Mayeda (1953) for oxygen; hydrogen measurements were carried out by reducing the water sample to hydrogen by passage over hot zinc at about 520 °C according to the procedure described by Coleman et al. (1982). The analytical precision is about ±0.1‰ for δ18O and ±1‰ for δD.

Raster layer and GIS analysis

Application of remote sensing data was made in order to evaluate characteristics of the studied recharge areas. All available information about the selected regional aquifers from previous studies were converted and geo-referenced into digital maps in geographical information system (GIS) environments, using Quantum-GIS 2.18 and SAGA GIS. Specifically, for the present study different products of digital remote sensing data, digital elevation models (DEMs), multispectral data and hyperspectral data were considered. Data from the DEMs, Advanced Spaceborne Thermal Emission and Reflection (ASTER) radiometer (Hayakawa et al. 2008) products with 30-m resolution, Shuttle Radar Topography Mission (SRTM) (Rabus et al. 2003) with 90-m resolution, and the Triangular Irregular Network of Italy (TINITALY /01, Tarquini et al. 2007) with 10-m resolution, were analysed—as an example, it was possible to assess the mean topographic elevation and the snow cover elevation, and create a hypsographic curve for each studied area. Multispectral and hyperspectral data were freely acquired from US Geological Survey (USGS) platforms for Landsat-8, Landsat-7, Sentinel-2, and the Moderate Resolution Imaging Spectroradiometer (MODIS) system to analyse the extension of snow cover. Images cover a timeframe from winter 2015 to spring 2016; to limit the atmospheric effect, only data images having cloud cover lower than 40% of the entire picture were considered, and of these, about 50 images were selected. Thus, at least four data images for each month were collected to improve data observations over time. Only for the high-medium resolution multispectral images obtained from Landsat-7-8 and Sentinel-2, a conversion to radiance value and then into surface reflectance with atmospheric correction was made by the following equation:

$$L_\lambda=M_\mathrm L\times Q_\mathrm{cal}+A_\mathrm L$$
(1)

where Lλ is the spectral radiance at the sensor’s aperture, ML is the band-specific multiplicative rescaling factor from the Landsat metadata; AL is the band-specific additive rescaling factor from the Landsat metadata, and Qcal is the quantized and calibrated standard product pixel values. Sentinel-2 images (level-1C) are already provided in the Top of Atmosphere Reflectance.

Images in radiance value can be converted to Top of Atmosphere Reflectance (ρp combined surface and atmospheric reflectance) by this formula:

$$\rho_\mathrm p=\left(\mathrm\pi\times L_\lambda\times d^2\right)/\left(\mathrm{ESUN}_\lambda\times\cos\theta_\mathrm s\right)$$
(2)

where Lλ is the spectral radiance at the sensor’s aperture, π is pi, d is the Earth-Sun distance in astronomical units (provided with the Landsat 8 metadata file), ESUNλ is the mean solar exo-atmospheric irradiances, and finally, 𝜃s is the solar zenith angle in degrees. All the values in the equation are provided by the metadata of each multispectral image.

To consider atmospheric attenuation for the calculation of a real surface reflectance, a dark object subtraction image-based technique (DOS1) was used (Moran et al. 1992):

$$\rho=\left[\mathrm\pi\times\left(L_\lambda-L_\mathrm p\right)\times d^2\right]/\left[T_\mathrm v\times\left(\left(\mathrm{ESUN}_\lambda\times\cos\theta_\mathrm s\times T_\mathrm z\right)+E_\mathrm{down}\right)\right]$$
(3)

Lp is the path radiance and Tv is the atmospheric transmittance in the viewing direction, Tz is the atmospheric transmittance in the illumination direction, and Edown is the downwelling diffuse irradiance.

In this case, using DOS1 correction, a few assumptions are made (Moran et al. 1992): Tv = 1 and Tz = 1, Edown is 0 and thus the resulting land surface reflectance is given by:

$$\rho=\left[\mathrm\pi\times\left(L_\lambda-L_\mathrm p\right)\times d^2\right]/\left(\mathrm{ESUN}_\lambda\times\cos\theta_\mathrm s\right)$$
(4)

It is worth pointing out that the accuracy of image-based techniques such as DOS1 correction is generally lower than physically based corrections, but they are very useful when no atmospheric measurements are available as they can improve the estimation of land surface reflectance.

After the aforementioned preprocessing workflow, the Normalized Difference Snow Index (NDSI) (Hall et al. 1995) was obtained from real reflectance values through raster calculator implemented in Q-GIS (Fig. 2). NDSI, with a fixed threshold of 0.4, allows one to discriminate pixels covered with snow or not and helps to discern snow/ice from cumulus clouds (Negi et al. 2009). Furthermore, the NDSI method can also help to delineate and map snow areal extension in the shadow areas in mountain zones (Kulkarni et al. 2002), which means that, in this case, the threshold value was fixed at 0.1. The NDSI formula is given in the following equation:

$$\mathrm{NDSI}=\frac{\mathrm{GREEN}-\mathrm{SWIR}}{\mathrm{GREEN}+\mathrm{SWIR}}$$
(5)

where GREEN is the green band and SWIR is shortwave infrared-bands. To optimize the raster images classification, the index S3 (Saito and Yamazaki 1999) has been used for Landsat-8 and Sentinel-2 only, adopting Eq. (6)

$$\mathrm S3=\frac{\mathrm{NIR}\left(\mathrm{RED}-\mathrm{SWIR}\right)}{\left(\mathrm{NIR}+\mathrm{RED}\right)\left(\mathrm{NIR}+\mathrm{SWIR}\right)}$$
(6)

where NIR and RED are, respectively, the surface reflectance of the near-infrared red and red band.

Fig. 2
figure 2

Example of image processing for S9 recharge area. a Landsat-8 OLI color composite with 6-5-2 bands after the pan-sharpening process, and b NDSI index surface obtained from reflectance values

Full size image

This process leads to identifying snow under vegetation, in areas where this might be difficult due to the presence of snow mixed with vegetation. For the studied aquifers, in most of the land surface above 800 m a.s.l. (where the snowpack can persist during winter and spring), naked-fissured rock, and rock with sparse thin soils are outcropping. Thus, the same workflow was carried out for each considered month, from December 2015 to April 2016, for all recharge areas. In this way, a GIS database with snow-cover shape files (over time) was created and analyzed.

Hypsographic curves

To define the mean topographic elevation of the recharge area (MTERA) and other features for each of the selected springs, hypsometric curves were delineated. The hypsographic profiles were developed using Q-GIS software (Nobedo 2.16.2 version) and TINITALY DEM with 10-m resolution cells. After defining the recharge area of each spring, these recharge areas were selected on the DEM, and using the QGIS zonal statistics plugin, MTERA values were calculated. Then every digital elevation model of the recharge area was transformed from a raster to a shape file using the polygonize function in QGIS, thusly each cell of the obtained digital elevation model was transformed into an independent polygon with respect to height. An “area” attribute was created in the attributes table of this layer, and data about elevation and area were exported in Excel format. All the records were organized from the highest to the lowest elevation and a new record “cumulative area” were calculated. The final hypsographic curves were plotted in a scatter-plot diagram comparing cumulative areas and elevation.

Results

Table 2 provides information about the MTERA obtained from the hypsographic curves, isotopic data referred to from previous studies, and the results of the 2016 sampling survey. The isotopic results (δ18O and δD) for each sampled spring were compared with mean isotopic values from previous studies (Barbieri et al. 2005; Desiderio et al. 2005; Tarragoni 2006; Agrillo et al. 2008; Petitta et al. 2009; Spadoni et al. 2009; Nanni et al. 2013; Banzato 2014; Tallini et al. 2014; Franchini et al. 2021) and unpublished data.

Table 2 Data about selected springs: MTERA (mean topographic elevation of recharge area), Summer Survey – June 2016, Autumn Survey – October 2016, Mean isotopic value from previous studies and related source of data
Full size table

For samples from the June 2016 monitoring survey (Table 2 and Fig. 3a) δ18O isotopic values range between –11.5 and –7.5‰, while δD isotopic values range between –78.2 and –46.7‰. Looking at October 2016 data (Table 2; Fig. 3b), δ18O and δD values vary respectively, between –11.3 and –7.6‰ and –47.8 and 76.2‰. Comparison between the two sampling surveys shows a very limited variability of isotopic signal with a maximum shift of ±0.2 for δ18O and ±1.1 for δD, and absolute variability values of 0.05 and 0.7 for δ18O and δD, respectively, highlighting the lack of influence of seasonal variability on the selected springs. This finding is also supported if data are compared with those obtained from previous studies (Table 2). In fact, the very low calculated standard deviation values show that data are clustered around the mean. Due to a limited available dataset, for some springs standard deviation values have been not determined, but even in this case the 2016 results and the mean isotopic values show a negligible mismatch. In addition, in both sampling surveys performed in 2016, S17 and S11 represent, respectively, the most enriched and the most depleted springs of the dataset.

Fig. 3
figure 3

δ18O–δD plots of sampling survey results from the a June 2016 and b October 2016 surveys and c previous studies isotope data. Sampling S3 was missed in the autumn survey due to the 2016 central Italy earthquake. For spring id refer to Fig. 1 and Table 1

Full size image

The recorded values are comparable with mean isotopic values obtained from previous studies (Fig. 3c), with only one outlier value characterized by a relatively more enriched signal. Data follow both the Central Italy Meteoric Water Line (C-IMWL) and the Italy Meteoric Water Line (I-MWL) proposed by Longinelli and Selmo (2003), having equations: 

$$\mathrm{C} \text{-} \mathrm{IMWL}\ \updelta \mathrm{D}=7.04{7\updelta}^{18}\mathrm{O}+5.608$$
(7)
$$\mathrm{I} \text{-}\mathrm{MWL}\ \updelta \mathrm{D}=7.61{\updelta}^{18}\mathrm{O}+9.21$$
(8)

The comparability of the datasets with local meteoric water lines is also proved by the obtained good index of agreement values (Willmott et al. 1985; Hauduc et al. 2015), which are 0.94 and 0.97, respectively for C-IMWL and I-MWL. Index of agreement values close to 1 indicate the best representativeness between the data and the model, while values close to 0 mean that the model is not representative of the dataset. An index close to 0 is instead obtained comparing the dataset with the Mediterranean Meteoric Water Line (M-MWL, Gat and Carmi 1970), highlighting the lack of influence of evapotranspiration on isotopic signal of the springs selected for this study.

In addition, the isotope data of all the analyzed data are generally depleted with respect to rain, but in line with snow isotope signals measured in previous studies. In fact, Barbieri et al. 2005 measured in two stations a mean rainfall value of –8.5‰ for δ18O and –55.7‰ for δD, while in the same area, Adinolfi Falcone et al. 2008 measured an isotope signal of snow and ice of –9.8‰ for δ18O and –65.0‰ for δD. In another recharge area, Spadoni et al. 2009 and Nanni et al. 2013 assessed, respectively, mean value of rain as –6.6‰ for δ18O and –45.3‰ for δD, and –4.1‰ for δ18O (δD values were not available).

Discussion

Isotope results of the summer sampling survey (June 2016) for each basal spring were compared to the literature on meteoric water lines referenced in previous studies. In addition, to obtain CIRE values for each spring, δ18O values were plotted on elevation–δ18O correlation lines adopted in literature for each spring or regional aquifer. Table 3 summarizes the literature correlation lines with the corresponding equations used for CIRE estimation. Specifically, CIRE values, obtained in this study, are between 806 and 2,265 m a.s.l..

Table 3 18O–Elevation correlation lines used for each spring (h= elevation in m a.s.l.). For spring id refer to Fig. 1 and Table 1
Full size table

Figure 4 shows four examples of calculation of isotopes-derived recharge area elevation (CIRE) for basal springs using correlation lines from the literature. Specifically, S2 has a medium discharge of 1.22 m3/s and a relatively limited extension of recharge area (about 41 km2). A CIRE value of about 1,780 m a.s.l. is estimated with Sibillini Mounts correlation line (Nanni et al. 2013), showing a surplus of more than 100 m with respect to the mean topographic elevation of 1667 m a.s.l..

Fig. 4
figure 4

Examples of CIRE calculation through literature correlation lines (black lines), for springs a S2, b S6, c S11 and d S16 (black dots). The red-dashed line identifies the mean elevation of the recharge area (MTERA). For spring id refer to Fig. 1 and Table 1

Full size image

S6 is fed by Giano-Nuria-Velino Mts regional aquifer, with an extension of 1,016 km2 (Boni et al. 1986; Petitta 2009). Peschiera (S6) is one of the most important springs of central Italy, with a mean discharge of 18 m3/s. The mean topographic elevation is about 1,225 m a.s.l., while the CIRE value, obtained by using the correlation line proposed by Spadoni et al. (2009), is 1,979 m a.s.l. showing a discrepancy between the two values of 754 m. The recharge area of S11 has been evaluated at 31 km2 with a mean discharge of 1 m3/s. For this spring the correlation line proposed by Nanni and Rusi (2003) for Majella aquifer has been considered to calculate the CIRE value of 2,265 m a.s.l.. Even in this case, the CIRE is higher than the mean topographic elevation of about 290 m. S16, instead, is fed by Greco-Genzana-Pizzalto-Rotella Mountains regional aquifer, with a medium discharge of 6.6 m3/s and a recharge area of about 269 km2 (Celico 1983). A CIRE value has been calculated by use of the Minissale and Vaselli (2011) correlation line, obtaining the value of 1,835 m a.s.l. for the June isotope data. The latter is 333 m higher than the mean topographic elevation of 1,502 m a.s.l..

The common feature of the aforementioned springs and generally of the entire spring dataset, is the positive discrepancy between isotopic recharge area elevations (CIRE) and the mean topographic elevations (MTERA). The comparison between these values clearly highlights this characteristic, which ranges from +106 m for S15 and +865 m for S5 with respect to MTERA (Fig. 5). This clear finding suggests that the contribution to aquifer recharge comes with different weights with respect to topographic distribution. In fact, the highest elevation areas are supposed to have more influence on the aquifer recharge.

Fig. 5
figure 5

Relationship between the mean topographic elevation of the recharge area (MTERA) and the computed isotope recharge elevation (CIRE). For spring id refer to Fig. 1 and Table 1

Full size image

As explained in the preceding section, isotopic values are representative of an average in space and time of the groundwater flow (Reddy et al. 2005), due to significant (but not quantified) residence time and the extension of the recharge area of each spring. Consequently, the application of 2016 isotope data (expressing groundwater flow of longer periods) would be misleading if compared with climate data (snow coverage and persistence) from the same year. Nevertheless, taking into account the found correlation between CIRE and MTERA and the stability of the isotopic values with time (testified by cited previous studies in each considered aquifer), the analysis conducted on an annual basis (2016) is considered representative of the relationships between recharge and discharge for each aquifer.

To shed more light on this finding, snow cover elevation and persistence during the monitoring period have been analyzed, using remote sensing information related to 2016 year. Specifically, analyses of multispectral and hyperspectral data allowed one to estimate the snow cover extension in the study area. Two parameters have been inferred and are summarized in Table 4: (1) the mean snow cover elevation (MSCE); (2) percentage of snow cover with respect to the total recharge area (evaluated for the month having the largest snow coverage). MSCE values have been assigned to each spring, based on snow cover extension and persistence on recharge areas, at different elevation ranges. These values have been estimated for each monitoring month, and a weighted average has been calculated to assign mean snow cover elevation (MSCE) values to each spring.

Table 4 Mean snow cover elevation (MSCE) and percentage of snow cover (SC) on recharge areas (*SC < 10%)
Full size table

Figure 6 shows the percentage of snow cover during the monitoring period on recharge areas considered for the selected springs. All investigated areas are generally characterized by the presence of snow cover with different distribution over time, according to the seasonal climate variability of Central Apennines. The maximum percentage of cover was reached in January 2016, followed by a general decrease from March 2016. The sole exceptions to these observations are S2, S11 and S16, which reach maximum values of snow cover, respectively, during February 2016, March 2016 and December 2015.

Fig. 6
figure 6

Histogram of percentage of snow cover of recharge areas during monitoring time. For spring id refer to Fig. 1 and Table 1

Full size image

Correlation between CIRE values and MSCE is well represented in Fig. 7. The correlation line (black-dashed line) does not include S3, S8, S14, S15 and S17 springs because snow coverage in the related recharge areas is lower than 10% (Table 4).

Fig. 7
figure 7

Correlation between calculated isotopic recharge elevation (CIRE) and mean snow cover elevation (MSCE); black dots represent recharge areas with snow coverage >10%, red dots represent recharge areas with snow coverage <10%. For spring id refer to Fig. 1 and Table 1

Full size image

For springs having at least 10% of snow cover in recharge areas, mean isotopic recharge elevation results become highly correlated with snow cover elevation showing a correlation coefficient (r) of 0.90, a coefficient of determination (R2) of 0.82 and a p-value strongly lower than 0.05, with an average discrepancy of about 100 m. It is worthy of note that the correlation line is parallel to the 1:1 line between CIRE and MSCE, confirming the strict dependence of isotope elevation values from the snow contribution to the aquifer recharge. In this framework, the ~100 m of discrepancy between the two lines would be evaluated as the uncertainty of the calculated values. The CIRE/MSCE correlation is significantly higher than the previous one between CIRE and MTERA (see Fig. 5) which showed an r-value of 0.46, R2 of 0.48 and a p-value slightly higher than 0.05, revealing the strong relationship between aquifer recharge and presence of snow cover. In terms of hydrogeological processes, this means that aquifer recharge and consequently spring discharge are clearly affected by snow contribution to the infiltration factor.

To further evaluate the different influences of snow cover on the isotopic signal of the selected basal springs, MTERA, CIRE and MSCE value positions were plotted on hypsographic curves of the recharge areas and Fig. 8 displays the findings for four of them, as representative of different relationships between the three values. MTERA logically falls in the middle of the hypsometric curves, while both CIRE and MSCE assume different positions, but in all cases, they fall in the first half of the hypsometric curves. A distinction in four classes can be inferred by this comparison, as explained in detail in the following.

Fig. 8
figure 8

Hypsographic curves of four basal springs compared with MTERA, CIRE and MSCE values. Springs a S16, b S11, c S2 and d S14. For spring id refer to Fig. 1 and Table 1

Full size image

The S16 recharge area shows a strong persistence of snowpack with time, reaching the maximum extension in December 2015 (45% of the recharge area, Fig. 6). Remote sensing analysis highlights that 28% of the recharge area lies at altitudes higher than the MSCE and this is reflected by the isotope signal of the spring. Indeed, the CIRE value is higher than the MSCE, pointing out that the melting of the long-lasting snow cover present at high altitudes clearly contributes to aquifer feeding (Fig. 8). This response is included in class A.

A similar behaviour is also recorded for S11, where the snow melting strongly affects the isotopic signal of collected water. In fact, as shown in Fig. 8, CIRE and MSCE values are almost the same with a negligible discrepancy of 24 m, suggesting that solid precipitation (snow) plays a fundamental role in aquifer recharge processes. The recharge area of this spring is snow-covered up to 5 months for more than 30%, with the maximum extension reached in March 2016 (Fig. 6). This persistence and extension over time results in a continuous aquifer feeding due to the snow-melting process, which first involves the thinner snowpack at lower altitude and then, progressively with the seasons, the thick layer located at high altitude. This response is identified as class B.

The S2 spring recharge area is representative of class C and is characterized by a considerable snow cover extension and persistence during the monitoring time (excluding April 2016) covering more than 45% of the investigated area (Fig. 6). The discrepancy among the different calculated elevations is limited, with MSCE values of 100 and 200 m respectively higher than the CIRE and MTERA values (Fig. 8). This finding confirms the influence of snow melting on the isotope signal and, at the same time, the relatively low mismatch between the snow and isotope elevations, which suggests that the infiltration process is influenced by the snow-covered portion of the recharge area, but aquifer recharge clearly also depends on rainfall coming from lower elevation.

The last class, D, is represented by S14, which shows both one of the lowest snow cover extensions (Fig. 8) and one of the lowest discrepancies between CIRE and MTERA (only 140 m). In contrast, the relatively high discrepancy measured between MSCE and CIRE (>400 m) indicates that snow cover has a negligible influence on aquifer recharge and consequently on isotopic signal. In this case, it is supposed that the aquifer is fed by classic rainfall contribution, with limited influence of the elevation gradient, not depending on snow cover peculiarity. Some springs belonging to this group are in fact characterized by relatively more enriched isotope signals, including the outlier value (see Fig. 3c), testifying to the major contribution of rainfall to the aquifer recharge.

Based on all collected information, the springs and related aquifers have been attributed to four classes, in terms of decreasing aquifer recharge dependency on snow cover:

  • Class A, strong dependency. S5, S7, S9 and S16 show a sharp discrepancy (more than 130 m) between mean snow cover elevation and CIRE, highlighting the dramatic influence of the persistence of snowpack in a large portion of recharge areas over time, and its consequent largest effect on aquifer recharge;

  • Class B, high dependency. S1, S4, S6, S11 S12 and S13, having nearly the same mean snow cover elevation and mean isotopic elevation values (discrepancy lower than 90 m), due both to the extension and persistence of snow cover during the monitoring time; in these cases, snow contribution to recharge is relevant and prominent, considering the distance with respect with the MTERA value;

  • Class C, significant dependency. S2 and S10 show an evident correlation between mean snow cover elevation and computed mean isotopic elevation, due to the persistence of snowpack over time and its impact on recharge rate;

  • Class D, negligible dependency. S3, S8, S14, S15 and S17, characterized by limited extension of snowpack (<10%) on the recharge area.

Basically, the strong correlation detected in most examined springs and aquifers highlights that snow cover represents a priority factor in the contribution to recharge of carbonate aquifers of the Central Apennines. From the hydrogeological point of view, the inferred reliance of aquifer recharge on snow cover and snowmelt infiltration processes can be characterized by classifying this relationship for most of the studied springs. This inferred dependence does not automatically affect the amount of aquifer recharge by considering a simple change of the nature of precipitation, from snow to rain, due to temperature rise. It is necessary to discuss whether an equivalent amount of precipitation as rain, relative to snow, can produce a similar infiltration rate or if this change would induce a reduction in aquifer recharge. Recent studies, mainly developed at higher latitudes, have verified how the shift from snow to rain is responsible for lowering the infiltration rate and consequent reduction in the groundwater resource. Cummings and Eibert (2018) noted that in North American mountains a lack of water storage in snowpack contributed to weak recharge of the local aquifer. Nygrel et al. (2020) affirm that large quantities of winter precipitation, stored as snowpack in the cold months, and released as meltwater in spring (Jasechko et al. 2017), clearly affects not only river discharge (Dierauer et al. 2018) but also groundwater recharge (Kløve et al. 2017). This evidence is supported also by stable isotope studies revealing that snowmelt is the main source of annual groundwater recharge (Jasechko et al. 2017). In more detail, Nygrel et al. (2020) describe a decrease in the water table over time and consequently in groundwater storage, attributed to shorter snowmelt periods and longer duration of high evapotranspiration rates, testifying to how higher temperatures (shifting from snow to rain) have driven the change, irrespective of precipitation changes. Finally, Wu et al. (2020) clearly stated that climate-driven impacts on groundwater-storage changes do not necessarily reflect the long-term trend in precipitation; instead, the trend may result from enhancement of evapotranspiration and reduction in snowmelt.

Therefore, this paper supports the hypothesis that a reduction of snow cover on carbonate fractured aquifers of central Italy results in a lowering of aquifer recharge and, of course, reduced spring discharge and groundwater resources. This inferred consequence requires careful evaluation of climate variations at a local and a global scale (Beniston et al. 2003; Kundzewicz et al. 2008; Chen et al. 2018; Opie et al. 2020; Jódar et al. 2021; Wu et al. 2020).

Conclusion

Stable isotope data collected in this study to evaluate the mean elevation of recharge areas of large aquifers, based on computed isotope recharge elevation (CIRE), clearly confirmed, at least for mountainous carbonate fractured aquifers, that recharge distribution is correlated with precipitation variation with altitude. The CIRE is found to be higher than the mean elevation of the recharge areas (MTERA) for all examined springs and related aquifers, a finding that was not obvious before analysis of the real data. For each recharge/discharge area there is strong correspondence between a high infiltration rate and the higher-elevation part of the aquifer. There can be a difference of greater than 300 m between the MTERA and CIRE. The implication is that in these cases, the aquifer recharge strongly depends (up to 75%) on the characteristics of higher portions of the hydrogeological basin.

This simple result has consequences for management of water resources, considering the potential effects of climate change, which is modulating precipitation and temperature not only throughout the year but also with elevation. Groundwater protection is also affected by the heavy contribution of recharge from high-elevation areas, and thus needs to take into consideration the concept of vulnerability correlated with travel time. Nevertheless, where fast flow contributes to only a minor percentage of the discharge, as happens in the study area, and the majority of the infiltration comes from far recharge areas, and the impact on groundwater vulnerability has to be weighted with respect to the recharge distribution with elevation.

Unfortunately, analytical evaluation of different altitude recharge zones is difficult on the basis of correlation between CIRE and MTERA, because of an inferred limited correlation (R2 = 0.5) of the two parameters. Consequently, looking at the general precipitation data, the relationship recharge/altitude is confirmed but its extent and weight remain unknown at this stage.

In any case, another and more novel finding has been inferred by this study, regarding the influence of snowpack in recharge mechanisms, as demonstrated by the correlation found between CIRE and mean snow cover elevation (MSCE) for the analyzed aquifers. The very high inferred correlation (R2 = 0.82) clearly indicates that the snowpack contributes to the infiltration rate more than the rainfall distribution with elevation. Of course, recharge happens along the whole permeable hydrogeological basin, but its distribution is strongly influenced by the snowpack in mountainous aquifers. In addition, considering not only the correlation but the similarity of CIRE and MSCE, expressed by a slope coefficient of the correlation equal to 1 (1.01), it can be affirmed that snowmelt represents the primary source of recharge on those areas, with different degrees depending on the aquifer characteristics. A preliminary distinction in four classes of snow cover distribution on the aquifer has been proposed. Based on the aforementioned recent studies, was the hypothesis that change in the nature of precipitation, shifting from snow to rain prevalence, can influence the infiltration rate due to both the karst nature of the aquifer and the persistence of the snow on the soil.

In such a conceptual model, future groundwater resources in these environments are subjected not only to a significant risk of depletion by the reduction of the extent of snow cover in Italy (at least for the Alps and Apennines), but also, in general, in the Mediterranean areas characterized by similar climate conditions and, at larger scale, the whole planet. In addition to the expected increase of evapotranspiration due to temperature rising, and of the supposed (but not unequivocally demonstrated until now) decrease of precipitation, the main pressure on groundwater availability for mountainous aquifers (where elevation is higher than the snow line) would be the reduction of snow cover, its extent and its seasonal persistence over time. The results inferred from the central Italy aquifers suggest that an impact of snow cover reduction in mountainous areas on infiltration rate cannot be excluded, even in a climate scenario with steady precipitation. Evidencing a possible dependency of spring discharge on snow cover contribution in Mediterranean climate conditions, one can conclude that it is necessary to launch a red alert warning system with respect to the future conditions of groundwater renewable resources.

To sum up, the obtained results have raised a potential risk for groundwater resource depletion in mountainous aquifers. In the Mediterranean region, these mountainous areas are the most important sources of water for drinking, environmental sustainability, and agricultural purposes. Additional efforts in research and in management are desirable, if not mandatory, for ensuring a sustainable approach to groundwater resource protection and use.