The Role of Snowmelt on the Spatio-Temporal Variability of Spring Recharge in a Dolomitic Mountain Group, Italian Alps

Abstract

Springs play a key role in the hydrology of mountain catchments and their water supply has a considerable impact on regional livelihood, biodiversity, tourism, and power generation. However, there is still limited knowledge of how rain and snow contribute to the recharge of Alpine springs. This study presents a four-year investigation of stable isotopes in precipitation and spring water at the scale of a 240 km² wide dolomitic massif (Dolomites, Italian Alps) with the aim of determining the proportions of snowmelt and rain in spring water and to provide insights on the variability of these contributions in space and time. Four precipitation sampling devices were installed along a strong elevation gradient (from 725 to 2660 m a.s.l.) and nine major springs were monitored seasonally. The monitoring period comprised three extreme weather conditions, i.e., an exceptional snowpack melting period following the highest snowfall in 30 years, an intense precipitation event (386.4 mm of rain in 48 h), and one of the driest periods ever observed in the region. Isotope-based mixing analysis revealed that rain and snowmelt contributions to spring water were noticeably variable, with two main recharge time windows: a late spring–summer snowmelt recharge period with an average snowmelt fraction in spring water up to 94 ± 9%, and a late autumn–early winter period with a rain fraction in spring water up to 68 ± 17%. Overall, during the monitoring period, snowmelt produced high-flow conditions and sustained baseflow more than rain. We argue that the seasonal variability of the snowmelt and rain fractions during the monitoring period reflects the relatively rapid and climate-dependent storage processes occurring in the aquifer. Our results also showed that snowmelt fractions in spring water vary in space around the mountain group as a function of the elevation of their recharge areas. High-altitude recharge areas, above 2500 m a.s.l., are characterized by a predominance of the snowmelt fraction (72% ± 29%) over the rain contribution. Recharge altitudes of approximately 2400 m a.s.l. also show a snow predominance (65 ± 31%), while springs recharged below 2000 m a.s.l. are recharged mostly from rain (snowmelt fraction of 46 ± 26%). Results from this study may be used to develop more accurate water management strategies in mountain catchments and to cope with future climate-change predictions that indicate a decline in the snow volume and duration in Alpine regions.

1. Introduction

High-elevation mountain catchments are important sources of freshwater stored as snow and ice [1,2]. A major part of their role is played by springs, which represent an irreplaceable source of drinking water even far beyond the foot of the mountains and provide water for many uses, including hydropower generation and artificial snowmaking. Catchments in mountain areas are very sensitive environmental systems [3] and are therefore considered sentinels of climate change [4]. Various studies predicted noticeable reductions of snow cover duration and snowmelt volumes in mountain areas as a consequence of climate change [5,6,7,8]. These alterations will significantly affect groundwater recharge and water yield of springs [9,10]. The major sources to groundwater recharge in high-elevation catchments are rainfall and snowmelt but their relative contribution has been rarely quantified. At high altitudes, where most of the snowpack accumulates, precipitation stations are absent or affected by large measurement errors [11,12]. In addition, extreme weather and remoteness of the high-elevation areas prevent regular accessibility [13]. In this framework, understanding recharge processes and quantifying the contribution of rain and snowmelt to spring water is necessary to properly manage groundwater resources and cope with future changes in the spring regimes. Another major challenge of mountain catchments is that springs can be fed by recharge areas located far from the emergence point and groundwater can follow complex and unpredictable flow paths. This is particularly true for carbonate aquifers as groundwater flows in fractures, conduits, and in the rock matrix causing a complex interaction of different flow behaviors.

Stable isotopes of hydrogen and oxygen are a valid method to trace groundwater flow and are ideal tools to distinguish between different water sources in mountain catchments because the natural abundances of these isotopes in the solid (snow) and in liquid (rain) precipitation differs appreciably, e.g., [14,15,16]. Furthermore, their labelling reflects diffused and long-term water infiltration processes [17] allowing to capture annual or inter-annual recharge variations at the regional scale. In recent decades, several studies used stable isotopes to identify water sources to groundwater recharge in mountain regions. The basic concept is that the isotopic composition of the input, such as rain and snowmelt, and of the output, such as groundwater samples from springs or wells, can be included in a mixing model to determine the proportion of the inputs that contribute to the output. For instance, based on this approach, Earman et al. (2006) found that snowmelt provided at least 40–70% of groundwater recharge in high-elevation sites in the Western United States [18]. Meng et al. (2015) quantified both the melting of snow and glaciers at the edge of the Tibetan Plateau in China and concluded that approximately 35% of the groundwater derived from ice–snow meltwater sources [19]. Penna et al. (2014) calculated a snowmelt contribution up to 64± 8% to spring water in a steep catchment in Northern Italy [20]. Understanding snowmelt contributions to groundwater is further complicated by isotopic fractionation occurring during snow formation, accumulation, ablation, and melting [21]. Therefore, the isotopic composition of snowmelt water infiltrating in the sub-surface may be different from the isotopic composition of the originating snowfall [22] and a direct comparison between the recharge input and the groundwater feeding the springs is not always possible. Prior to melting, sublimation processes cause the exchange of the water vapor of the superficial layers of the snow with atmospheric moisture, producing an isotopic enrichment of the remaining snow [23,24]. During the melting season, Taylor et al. (2001) studied the isotopic evolution of the melting composition and found a progressive enrichment of the snowmelt and a dependence of snowmelt δ¹⁸O on melt rates [25]. Lee et al. (2010) showed that the initial isotopic composition of a snowpack is an important component in determining the temporal changes of the isotopic signature of both the snowpack and the melting waters [26]. In addition, melting–freezing processes before and during the snowmelt period can result in significant isotopic offsets between the original snowpack composition and that of snowmelt [27]. In general, all these processes cause slope changes in the δ¹⁸O–δ²H relationship for the meltwater samples compared to the slope of the meteoric water lines.

Often, previous works focused on single catchments and/or a few springs or groundwater wells, and very little is known about how snowmelt contributions to spring water vary in space (including along elevation gradients) and time (over seasons and years) at the scale of an entire mountain group. Even less knowledge is available on the role of snowmelt in feeding springs in karst mountain groups and catchments, where the network of fractures and conduits strongly influences the infiltration processes [28,29,30] and make the hydrological response pretty complex [31].

In this work, we used the stable isotopes of oxygen and hydrogen in precipitation and spring water sampled along a steep gradient (from 725 to 2660 m a.s.l) to investigate the spatio-temporal variability of rain and snowmelt contribution to the recharge of several springs with an estimated recharge area of over 120 km² in a karst mountain group (i.e., at the scale of the entire hydrostructure [32]) in the Dolomites, Italian Alps. We aim to address the following specific questions:

• What is the variability of the precipitation isotopic signal in the study area and how is this variability transferred to spring water?
• What is the fraction of snowmelt and rain in spring water and how does it change in space and time in this high-elevation mountain group?
• What can we infer on spring recharge dynamics in the study area?

2. Study Area

2.1. Geographical Setting

The Pale di San Martino, Pale di San Lucano, and Mt. Agner (hereinafter referred together as “the Pale Mts.”) form a 240 km² wide mountain group located in the Dolomites (Eastern Italian Alps). The Pale Mts. are characterized by a typical high-Alpine landscape with narrow and deep valleys, steep peaks and vertical rock faces that abruptly switch to more gentle slopes at the base of the mountains. The central part of the area is covered by a 50 km² bare karstic high-elevated plateau, known as “Pale di San Martino Plateau” with an average altitude of 2650 m a.s.l. The highest peak is Cima Vezzana, which reaches 3192 m (Figure 1). Five main valleys deeply cut into the massif: Travignolo, Focobon, and Liera Valleys, located in the northern side, Angheraz–San Lucano Valley system to the east, and Pradidali–Canali Valley system to the south. From a climatic point of view, snow usually begins to accumulate in December and lasts until the beginning of spring in the valleys and up to late summer on the highest peaks and on the plateau. Precipitation and temperature are extremely variable due to the complex topography. For the period 2006–2016, the mean annual precipitation at the Col di Prà station (located in the upper San Lucano valley at 863 m a.s.l.; Figure 1) was approximately 1650 mm and the mean temperature 8.0 °C (official meteorological data available at www.arpa.veneto.it). For that 10-year period the highest annual precipitation was recorded in 2014 (2342 mm) while the lowest in 2015 (1076 mm).

2.2. Hydrogeological Setting

Geologically, the rocks outcropping in the Pale Mts. span from Ordovician to upper Triassic [33]. Most of the area is covered by dolomitic rocks formed from the recrystallization of Middle-Triassic reef buildups (Carbonate complex in Figure 1). These deposits are over 1000 m thick and form the main body of the mountains and the main regional aquifer. Matrix porosity is low whereas secondary porosity is high due to the numerous faults and joints that crosscut the entire thickness of the carbonate massif. Below the carbonate rocks, the stratigraphic sequence is characterized by Upper Permian to Lower Triassic deposits, including shallow marine terrigenous deposits, mainly silts, clays and evaporites which act as the main regional aquiclude (terrigenous and evaporitic complex, Figure 1). All the sedimentary cover lies on top of a rhyolite porphyry (Porphyric complex), which in turn is deposited over the Paleozoic metamorphic basement (Metamorphic complex). The contact between the carbonates and the underlying terrigenous and evaporitic complex acts as a no-flow boundary for the groundwater flow, causing the emergence of the springs [32]. Due to the numerous tectonic phases which affected the area and in particular to the Alpine compression [34,35], this boundary assumes a complex geometry providing very different spring elevations around the mountain group.

Heterogeneous quaternary deposits, including rockfall deposits, alluvial, and glacial sediments fill up the valleys and frequently cover the aquifer–aquiclude contact causing a diffuse emergence of the springs from the talus. Some springs outflow directly in the streams as gaining reaches [32,36] therefore they cannot be clearly recognized as a single outlet. These springs have not been considered in this study because we tried to sample as close as possible to the emergence point in order to capture only the groundwater contribution, and also to avoid possible secondary fractionation processes occurring in the stream in case of very low flow. Hereinafter, we will refer as “springs” considering only point springs. ANG and PRA springs are located respectively to the east and to the south of the Pale di San Martino Plateau and show the highest discharge (

Table 1. Summary of the monitoring points and sampling period. Refer to Figure 1 for the location on the map. Springs sampled seasonally refer to springs included in a periodical monitoring, while springs sampled occasionally refer to springs with fewer samples. Mean discharge values were taken from Lucianetti et al. (2019) [32].

Code	Name/Location	Elevation m a.s.l.	No. of Samples	Mean Discharge (L/s)	Sampling Period
Precipitation
SC	Scalette	725	35		September 2014–December 2018
VW	Villa Welsperg	1038	36		September 2014–December 2018
SM	San Martino	1470	33		January 2015–December 2018
RS	Rosetta	2660	20		September 2014–September 2016
Snowmelt patches	Rosetta Plateau	2650	2		July 2014 and July 2016
Periodically monitored springs
TRE	Treviso	1580	8	30	July 2014–December 2016
NER	Acque nere springs	1150	7	140	December 2014–December 2016
SCA	Scalette	725	6	10	July 2014–January 2016
LUC	San Lucano	900	6	130	July 2014–January 2016
ANG	Angheraz sx	1018	9	350	July 2014–December 2016
PRA	Pradidali springs	1456	8	210	July 2014–December 2016
TRA	Travignolo	1960	8	50	July 2014–December 2016
FOC	Focobon	1720	6	100	July 2014–July 2016
GAR	Gares	1729	6	120	July 2014–July 2016
Occasionally monitored springs
FON	Fontane Fosche	1122	5	30	December 2014–January 2016
POL	Polver	903	2	10	January 2016–July 2016
LEF	Le fontane	935	1	80	January 2016
ANG-DX	Angheraz dx	980	1	150	January 2016
VEC	Vecia	1290	5	20	December 2014–January 2016
ALB	Alberghet	1560	1	80	July 2016

). Both these springs are characterized by a main outlet and by secondary emergence points that are temporarily active during high-flow periods or after heavy rainfall events. A minor part of the discharge outflowing from the carbonate massif drains in the northwestern direction, towards the GAR, FOC and TRA springs (Figure 2a,b). Other springs, such as NER, TRE, SCA and LUC are located in the southeastern side of the mountain group (Figure 1). All the monitored springs are permanent and only some of them are used all year around for drinking water supply.

3. Materials and Methods

3.1. Monitoring Network and Sampling

Available data, such as spring cadasters and preliminary hydrogeological surveys, were used to design the stable isotopes monitoring network. Table 1 reports the location name and the type (i.e., precipitation or spring) of each monitored point. Given the numerous springs listed in the official cadasters (more than 500), we selected the ones to monitor based on two criteria: (i) to characterize the majority of the groundwater and thus focus on the largest springs; and (ii) to characterize every portion of the massif, thus opting for a well spatially-distributed monitoring network. Overall, a total of nine springs were selected for periodical monitoring, although additional samples were acquired occasionally from other springs to expand the database. When large springs were located close to others, only the largest one was chosen for seasonal monitoring, such as in the case of ANG and ANG-DX springs (Figure 1). This provided a reasonable average composition of the groundwater at the massif scale. Spring sampling was carried out between July 2014 and December 2016, approximately every three or four months depending on spring accessibility and meteorological conditions for a total of 79 samples. Precipitation, both in solid and liquid form depending on the season, was collected at four sampling stations (see Figure 1 for the location of the stations and Table 1 for their elevation). The monitoring network included also a high-elevation station located at 2660 m a.s.l. in the Pale di San Martino Plateau.

Snowmelt sampling for isotopes analyses can be performed using different methods, such as passive capillary samplers [37,38], lysimeters [39], or snow pits [40]. These sampling techniques require maintenance and can be costly, becoming unsuitable for remote areas, such as high-elevation mountain catchments. Snowmelt samples need to be preferred over snow samples due to the modification of the isotopic composition of snowmelt with respect to that of the original snowpack [25,26]. However, previous studies have shown that unheated precipitation samplers characterized by “extended funnels” represent a valid alternative to these methods and are sufficient to capture the general isotope signal over the course of a snowfall and snowmelt period [18,41]. According to Earman et al. (2006) samples from the extended funnel collectors were “[…] more representative of infiltration water isotope composition than a fresh snow sample because the snow in the extended collector is subject to many of the snow metamorphism effects that impact snow on the ground” [18]. Despite this method is not accurate when considering melting isotopic variations at the daily time resolution, it can be considered sufficient to capture the general isotope signal over the course of the longer seasonal melting periods [18]. Considering the remoteness and the temporal resolution of this study, in this paper we used extended funnel gauges samplers and we considered them representative of the isotope signal of snowmelt. A total of 124 bulk precipitation samples were collected approximately at monthly intervals between September 2014 and December 2018. During each precipitation sampling occasion, the amount of water in the precipitation collector was also measured. These bulk amounts were used in the study for all the calculations involving precipitation, as they showed to be more accurate than the official weather stations, especially during the winter months. The high-elevation station was monitored only for two years due to the difficult accessibility during the winter months. Sampling of rain and snow was carried out using a homemade device, built following the IAEA isotope hydrology laboratory design [42]. The device consists of a high-density polyethylene bucket with an inner tube connected to a 20-cm funnel and an outer tube used for pressure equilibration. This device was proved to prevent evaporation without the need to use paraffin. In winter periods, a cylindrical tube, also known as “extended funnel gauge” (Figure 2c) was placed on top of the funnel to also collect snowfall. Snowfall was sampled without heating and was transferred in a sealed bucket to guarantee melting without evaporation and was left to melt at roughly 20 °C. Only two snowmelt samples collected directly from the water dripping from residual winter snowpack are available.

A widespread availability of rain and snowmelt samples collected at different times and locations is advisable in complex and rugged terrain to characterize the infiltrative input into the sub-surface. Unfortunately, our sampling collection was limited by the accessibility of the area, which proved to be challenging and required the use of helicopters on a few occasions. Consequently, we had to select a limited number of sampling stations. Moreover, sampling at higher temporal frequency (i.e., at least once a week or even more frequently) might have allowed us to explore some short-time responses in tracer concentration and distinguish between quick-flow dynamics in the karst conduits and slower groundwater circulation in the dolomitic fractures. However, the purpose of the research was to represent the dominant flow processes at the seasonal and regional scale without claiming to distinguish the different flow components.

3.2. Laboratory Analysis

Water samples were analyzed for their ¹⁸O and ²H composition at the isotope geochemistry laboratory at the University of Parma, Italy. The oxygen isotopic composition was measured by means of the water–CO₂ equilibration technique [43] using the Finnigan MAT HDO automatic equilibration device. The same device was used for the measurement of the deuterium/hydrogen ratios, equilibrating directly ultrapure hydrogen with water by means of a catalyzer. A Finnigan Delta Plus mass spectrometer was used to carry out the isotopic measurements that are reported here in terms of δ units (‰ deviation of the isotope ratio from the V-SMOW isotopic standard). The laboratory standards were periodically calibrated, as recommended by the IAEA, versus the international water standards V-SMOW, GISP, and SLAP. Instrumental precision is ± 1.0 ‰ for δ²H and ± 0.05 ‰ for δ¹⁸O.

3.3. Climate and Spring Discharge Data

Temperature data used in this work are publicly available and were acquired by different Weather Agencies: Trento Province stations (data available at www.meteotrentino.it) were used to represent the temperature at our SM and VW stations, Veneto Region Environmental Protection Agency data (ARPAV; www.arpa.veneto.it) were used for the temperature at SC and the station installed by Meteotriveneto Association (www.meteotriveneto.it) was considered representative for the temperature at RS station. Precipitation amounts used in this study generally refer to the monthly bulk volumes measured in our precipitation sampling devices with the exception of long-term time series or daily precipitation data which were taken from the Col di Prà station (ARPAV), which is the station that has the longest record. Daily snow depth data were provided by ARPAV at the Malga Losch station. Average discharge data were taken from Lucianetti et al. (2019) and were integrated with additional unpublished discharge data of three of the largest monitored springs, to compare them with the isotope values deriving from the monitoring network. Discharge measurements were performed measuring the spring sectional area and the flow velocity by an OTT Nautilus C2000 current meter equipped with a Senza Z300 velocity indicator. The discharge measurements are affected by a systematic error of about 10–11% [44].

3.4. Data Analysis

3.4.1. Dual Isotope Analyses

A linear-regression analysis was applied to δ²H and δ¹⁸O values of precipitation samples collected between 2014 and 2018 to determine the Local Meteoric Water Line (LMWL). The deuterium excess (hereafter d-excess) of each sample was calculated following Dansgaard (1964) [45]:

d-excess = δ²H − 8 δ¹⁸O

(1)

Graphically, the d-excess represents the slope of the Global Meteoric Water Line (GMWL). Changes in the slope of the LMWL compared to the GMWL indicate fractionation processes during the evolution of the air masses, from their origin to the study site. For example, a differential isotopic fractionation of ¹⁸O with respect to ²H can be associated with evaporation during rainfall or sublimation during snowfall, resulting in a lower slope for the LMWL. d-excess was used to study these effects in the study area. The δ²H and δ¹⁸O of spring water were plotted against the LMWL (Figure 3). The isotopic signal transmitted to groundwater is that of recharge and does not necessarily reflect the precipitation range [46]. Therefore, differences or analogies between the isotopic composition in precipitation and spring water can give information on recharge sources and processes. Stable isotopes measured in precipitation at different elevations were used to estimate the elevation of the recharge area feeding the springs. For the application of this method it is assumed that spring water represents an integrated sample of all the recharge along the slope to a certain distance upgradient to the sampling point [47]. To evaluate the average recharge altitude of each spring we used the relationship between δ²H in precipitation vs. elevation and the monitoring period in which all the four stations were functioning (September 2014 to December 2016). This relationship was obtained by performing a linear fitting for each monthly campaign and then averaging all the fitted lines to obtain a mean altitude gradient line. Standard error was used to identify the level of uncertainty in the fitting. Recharge elevation were estimated both mathematically and graphically by matching the average spring isotopic composition for the same time period with the precipitation gradient line. Given the numerous assumptions, recharge altitudes computed with this method should be considered to be a rough estimate, but they still provide useful ranges to compare the different springs throughout the massif.

3.4.2. Seasonal Variability of Spring Water Isotopic Composition

The seasonal variation of the isotopic composition in precipitation represents an input signal that infiltrates in the sub-surface, recharges the aquifer, and contributes to the catchment hydrological response. The different amplitude between the seasonal isotopic variation of precipitation and groundwater and the shift between the two signals have been used in the literature to estimate the mean transit time of the groundwater [17,48]. The method is based on fitting the isotopic signal to a sine function of appropriate amplitude and frequency and calculating the phase shift and difference in amplitude between the input and output. Due to the sampling gap in our spring dataset this fitting was not possible. Hence, to quantify the isotope damping and to compare it among the sampled springs, we used a simple statistical coefficient. Coefficients have been successfully adopted as proxy for catchment travel times by previous authors, when the dataset was limited or presented gaps (e.g., [49,50,51]). Analogously, we used the Damping Ratio (DR, Soulsby, 2015 [50]), defined as the ratio of the coefficient of variation (CV) of the damped signal to the CV of precipitation. In previous studies the output signal was referred to stream water, while in our case it consists of the variation of spring water, as follows:

DR = CV springs/CV precipitation

(2)

DRs were calculated for the periodically monitored springs using the volume-weighted precipitation measured at the highest stations (RS). Since the interest of this approach was to compare the DRs between different springs and not their absolute value, we report and discuss this calculation only for this station.

3.4.3. Two-Component Mixing Model

To estimate the contribution of snowmelt and rain to the recharge of spring water, a two-component mixing model was performed [52]. The choice of a simple two-component mixing model was motivated by data availability and by the intention to perform preliminary tracer-based analysis of the spring water system in the study area through a relatively simple approach. This approach somehow assumes that there is a single underground reservoir that is recharged by precipitation and that supplies water to the springs. This is of course a strong simplifying assumption, especially in karst hydrological settings, where quick-flow components in the karst conduits can be present. However, considering the seasonal time scale of this work, this approach can be deemed reasonable. In addition, the two-component mixing model assumes that spring water results exclusively from the mixture of the two sources (endmembers). This condition is not met when exchanges occur with adjacent aquifers or spring water receives inflows from external water bodies. However, the dolomitic aquifer in our study area lies on an impermeable aquiclude throughout all its extension, creating an isolated hydrostructure without exchanges with neighboring aquifers, and it can, as such, be considered a sort of individual simple reservoir.

The computation of the snowmelt and rain contribution to spring water is based on the tracer mass balance, as follows:

1 = sm + r

(3)

δ_spring = sm × δ_sm + r × δ_r

(4)

where sm and r represent the spring water fraction due to snowmelt and rain, respectively and δ refers to the isotopic composition of each component. In this case, both δ¹⁸O and δ²H were chosen for the mixing model, but given the similar results obtained from the use of the two tracers only the results for δ²H are reported. Equations (3) and (4) can be solved to obtain the snowmelt fraction, as follows:

sm = (δ_spring-δ_r )/(δ_sm-δ_r )

(5)

Mixing models in snow-dominated catchments require particular care in the definition of the endmember compositions [53]. For this reason, first of all we used a temperature-based method to decide whether to attribute the bulk precipitation sample to the snowmelt endmember or to the rain endmember. We partitioned the measured bulk precipitation amount measured in our sampling devices (mm) in rain and snow based on a melt factor (F) which in turn depends on the average monthly temperature measured at the station [54]:

rain = F × P

(6)

snow = (1 – F) × P

(7)

A simple approach to calculate F is:

F = 0   T ≤ 0 °C;

(8)

F = 0.167 × T   0 °C < T< 6 °C

(9)

F = 1    T ≥ 6 °C

(10)

where P is the bulk monthly precipitation amount, and T is the average monthly temperature recorded at the given station. In this way, when average monthly temperature was above 6 °C the entire precipitation was considered to be rain, when average monthly temperature was below 0 °C all precipitation was considered to be snow [54]. When temperature was between 0 and 6 °C rain and snow were partitioned based on the temperature value. Intuitively, temperatures closer to zero resulted in larger snow proportions and vice versa. Using the obtained monthly rain and snow values, when the rain proportion was greater than the snow proportion we attributed the sample to the rain endmember, when the rain proportion was smaller than the snow proportion we attributed the sample to the snowmelt endmember.

The next step of the mixing modelling analysis involved the choice of using time-invariant or time-variant isotopic values for the endmembers [55]. Using a time-invariant isotopic values for the endmembers (i.e., the average value of precipitation isotopic composition collected during the monitoring period) would inadequately represent seasonal variations in the precipitation input, implying that springs could be recharged not only by precipitation fell before the sample was collected but also by precipitation fell after the spring sample was collected, which is of course physically not meaningful. Therefore, we defined time-variant endmembers based only on precipitation data antecedent to the time when the sample was collected. Only for the first spring sampling campaign (July 2014), for which antecedent rain samples were not available, we used average rain values for the entire period. Since spring sampling occurred between July 2014 and December 2016, the precipitation data used in the mixing model refer to the same period to have a coherent dataset. We used two different mixing scenarios based on the gauging stations included in the mixing model. The first one considered all the precipitation stations (scenario 1), while the second one (scenario 2) considered only the two stations at highest elevations (SM at 1470 m and RS 2660 m a.s.l.), which are assumed to well represent the recharge altitude of the springs in this hydrogeological setting (the stations are located at elevations above most of the springs). Uncertainty in the separation of the two components at 70% confidence was estimated using the method of Genereux (1998) [56]. The more the components in the mixing model are isotopically similar and the larger their variability, the greater the uncertainty. In this study, the two components are clearly distinguishable (see Figure 3), but the variability (standard deviation, SD) of each endmember is high due to the intrinsic variability of precipitation in such a rugged topography.

4. Results

4.1. Dual Isotope Analysis

Precipitation samples show a large variability during the four years of sampling, ranging between −131.2 ‰ and −19.8 ‰ for δ²H and between −17.28 ‰ and −3.47 ‰ for δ¹⁸O (Figure 3). Snowmelt samples collected from the extended funnel gauges are depleted in heavy isotopes relative to rain and clearly distinguishable on the δ¹⁸O−δ²H plot. Snowmelt collected from the snowmelt patches overlaps with snowmelt samples from the extended funnel gauge, indicating a comparable composition.

The LMWL build for the study site shows a strong linear correlation between δ¹⁸O and δ²H (R² = 0.99, p < 0.001; n = 124), with a slope of 7.9 and deuterium intercept of 9.8. This equation is very similar to the Global Meteoric Water Line (GMWL, Craig, 1961 [57]) and to other water lines found for Northern Italy sites (Longinelli and Selmo, 2003, y = 7.7x + 9.4 [58]; Penna et al. 2014, y = 8.1x + 10.3 [16]). Penna et al. (2017) also defined a LMWL in the Dolomites, but the intercept found in their study is slightly higher (13.5) [20]. All the δ¹⁸O and δ²H pair values measured in springs (n = 79 samples) plot in a narrow range relative to the LMWL (Figure 3). Overall, the isotopic signature of the springs appears less scattered than the isotopic composition of precipitation and shows a much smaller variation, with a range in δ²H between −71.0‰ and −99.0‰. Spring samples fall closer (i.e., are more similar in their isotopic composition) to snowmelt samples than to rain samples. The isotopic signature of precipitation shows a distinct variation with elevation. The δ²H in precipitation vs. elevation relationship is reported as an inset in Figure 3. During the computation of the gradient line three samples at the highest station showed a more enriched composition compared to the samples at lower altitude stations. This produced an increase of the isotopic composition with increasing elevation (reversed gradient) and overall a decrease in the slope of the mean gradient line for the study area. By analyzing these three samples we observe that they are all snow samples with lower d-excess compared to the average of the other precipitation samples (respectively +7.5‰, + 11.0‰ and +8.4‰ vs. an average of the other samples of +13.1‰) suggesting the occurrence of secondary fractionation processes in the snowpack, such as snow sublimation which is known to reduce d-excess values [23]. Including these samples in the computation of the elevation isotopic gradient produced overestimated and unrealistic recharge elevations, i.e., above the maximum elevation of the study area and thus we decided to exclude them from the computation. Removing these enriches samples returns a mean gradient line that matches very well the topographic elevation of the massif (inset in Figure 3). The computed lapse rate of −1.0 ‰ in δ²H per 100 m corresponds to the lapse rate reported by Flaim et al. 2013 for the Trentino–South Tyrol Region of the Eastern Alps (−1.04 δ²H %/100 m) [59] and is similar to the lapse rates reported by previous authors for other high-elevation sites (-1.2 δ²H ‰/100 m found by Jeelani et al. 2003 in Kashmir Himalayas [60]; −0.8 δ²H ‰/100 m found by Paternoster et al. 2008 for the Vulture volcano in Southern Italy [61]) and lower than the lapse rate of −1.6 δ²H ‰/100 m found by Penna et al., 2014, in a mountain valley in South Tyrol [16]. Using this gradient, we estimated the recharge elevations (

Table 2. Mean elevation of the recharge areas of the periodically monitored springs calculated using the δ²H vs. elevation in precipitation (inset in Figure 3) and the mean δ²H in the springs. Values in parenthesis indicate the minimum and maximum elevation, respectively, considering (i.e., adding and subtracting) the computed standard error.

Spring	Mean Recharge Elevation (m a.s.l.)
TRA	2710 (2143–3391)
FOC	2601 (2041–3273)
GAR	2535 (1980–3202)
ANG	2429 (1881–3087)
PRA	2372 (1827–3025)
SCA	1937 (1420–2553)
NER	1856 (1345–2466)
LUC	1853 (1342–2463)
TRE	1782 (1276–2386)

).

Figure 4a shows a box plot of the δ²H ‰ composition in precipitation at each station. A clear elevation effect is visible also in this representation, with the mean of the highest station showing the most depleted value (δ²H ‰ = −77.5 ‰) and the lowest station exhibiting the most enriched mean composition (δ²H ‰ = −60.9 ‰). All the stations are characterized by a high variability in isotopic composition due to heterogeneous nature of precipitation, which changes from rain to snow repeatedly throughout the year. Figure 4b shows the variation of d-excess in precipitation calculated for each station. d-excess values and their variability increase with elevation, up to a mean value of +12.5‰ calculated for the highest station (RS, 2660 m a.s.l.). At SM station (1470 m a.s.l.) the mean is +10.5‰, while it decreases to +9.5‰ for the VW station (1038 m a.s.l.). However, this relationship with elevation does not hold for SC station, which shows values similar to VW but with a greater variability. The average value of d-excess in all the spring samples is +12.2‰ with the highest absolute value found in the spring outflowing at the highest elevation (TRA spring, 1960 m.a.s.l., +14.9). On average, the d-excess values of spring waters are comprised in the range of d-excess values calculated for the two high-elevation precipitation stations.

4.2. Analysis of Seasonal and Spatial Variations of δ²H

The time series of δ²H in precipitation and in the monitored springs are represented in Figure 5 to highlight seasonal variations during the monitoring period. The monitoring period comprised very different weather conditions, including an important rainfall event in November 2014, (the 48 h cumulative rainfall was 386.4 mm) and a dry period during November and December 2015 (64 days without rainfall, Figure 5a). The first spring sampling campaign was performed during an exceptional melting period which reflected the exceptionally thick snow cover of the previous 2013–2014 winter season (see snow depth in Figure 5a). At the Malga Losch station (Figure 1) the maximum monthly snow depth value measured in March 2014 was 368 cm, the highest value ever recorded since the beginning of the monitoring (year 1989). The station of Arabba, close to the study area, also experience the highest value of snow depth since 1980 (https://www.arpa.veneto.it/temi-ambientali/climatologia/dati/inverno-2013-14#clima).

The δ²H composition of the precipitation samples shows a large seasonal variability (Figure 5a), with the highest range at the SM station (−102.3 ‰ between the minimum and maximum value). As expected, more depleted isotopic composition in precipitation was recorded during winter months, such as March 2015, while summer composition shows enriched values. The isotopic composition of precipitation at the highest station (RS, 2660 m a.s.l) is generally more depleted than at the other sampling points, consistently with the well-known isotope altitude effect in precipitation [62]. SC station, which is the lowest station among all, exhibits an enriched signature, with the highest recorded value of −19.8‰. December 2015 was characterized by the complete absence of precipitation in all the stations and during the previous month, November 2015, only two stations recorded a very small amount of precipitation (22 mm in the SC station and 75 mm in the RS station). Also in December 2016 cumulative monthly precipitation was null. The isotopic signature of the springs is significantly damped compared to that of precipitation, with the greatest seasonal range of −19.8‰ recorded at the TRA spring. Despite the clear decrease in the amplitude of the spring signal compared to precipitation, a distinctive seasonal variation is still recognizable (Figure 5b). Results of the damping ratio quantification (DRs) are reported in

Table 3. Damping ratios (coefficient of variation of springs to precipitation) calculated using the δ²H precipitation time series of the highest station (RS, 2660 m a.s.l.) and the periodically monitored springs.

Spring	Damping Ratio (DR)
TRA	0.21
FOC	0.20
GAR	0.18
LUC	0.16
SCA	0.16
PRA	0.12
TRE	0.11
ANG	0.10
NER	0.09

. The lowest damping, corresponding to the largest DRs, are found in the TRA, GAR, and FOC springs with a maximum value of 0.21 for TRA. These springs are in the northern side of the mountain group and are characterized by the highest elevations. The smallest damping occurs in the NER spring and in the ANG springs which are located to the south and to the east of the massif, respectively.

Results from discharge measurements (Figure 5c) show that for the largest springs of the Pale Mts. (ANG and PRA) the highest discharge values were measured during the months of July (July 2014 and July 2016). Conversely, the lowest values were recorded in January 2016, during the monitoring survey that followed two dry months. Thus, these discharge values likely represent baseflow conditions for these two springs. Also, May 2015 shows relatively low values as the beginning of the melting period occurred later and some springs were still in low-flow conditions. For the TRA spring, discharge data exhibit the maximum in July as well. In this case, very low discharge values are already recorded in December as the spring is located at very high elevation (1960 m a.s.l. compared to 1018 and 1456 m of ANG and PRA, respectively) and snowfall may occur earlier compared to other parts of the Pale Mts., preventing recharge. Discharge values measured in December and November change appreciably from year to year, depending on the intensity of autumn rainfall events, but are generally lower than summer high flows.

4.3. Mixing Model Results

A summary of the time-variant isotopic composition of the endmembers used in the two-component mixing model is reported in

Table 4. Summary of the time-variant isotopic composition of the rain and snowmelt endmembers of the mixing model used to calculate the snowmelt and rain fraction for each spring sampling campaign for the two different scenarios. For the first spring sampling campaign (July 2014), for which antecedent rain samples were not available, we used average rain values for the entire period.

Date	Scenario 1				Scenario 2
	Rain δ2H (‰)		Snowmelt δ2H (‰)		Rain δ2H (‰)		Snowmelt δ2H (‰)
	Weigh. avg.	SD	Weigh. avg.	SD	Weigh. avg.	SD	Weigh. avg.	SD
14 July	−54.9	16.4	−91.2	na	−55.1	15.5	−91.2	na
14 December	−61.5	10.5	−117.1	14.1	−63.8	8.1	−105.3	10.0
15 May	−60.6	11.6	−102.4	16.8	−64.5	6.7	−102.5	18.5
15 July	−52.4	13.7	−102.4	16.8	−52.6	12.5	−102.5	18.5
15 November	−55.5	13.9	−98.0	17.2	−57.2	12.6	−94.8	18.8
16 January	−55.5	13.9	−92.5	21.3	−57.2	12.6	−91.9	20.4
16 July	−56.8	13.4	−89.3	19.6	−57.5	12.3	−88.2	18.4
16 December	−54.9	16.4	−90.8	19.5	−55.1	15.5	−90.4	18.3

. Standard deviation for scenario 1 (n = 4 precipitation stations) is slightly higher than for scenario 2 (n = 2 precipitation stations) due to the increased number of samples considered. Weighted averages are very similar indicating that increasing the number of stations from two to four does not significantly change the results. This is because high-elevation stations receive more precipitation amount and thus have a greater impact on the weighted average isotopic composition compared to valley-floor stations. This suggests the importance of a monitoring network which includes high-elevation stations to characterize the precipitation isotopic composition of a mountain massif, while lower altitudes stations are less relevant.

Panel (a) and (b) in Figure 6 show the snowmelt fraction from the mixing model for scenario 1 and 2, respectively. Results were grouped according to the spring sampling campaigns (each box represents a sampling campaign), to show the seasonal variation of the snowmelt fraction in the monitored springs. Uncertainty related to the mixing models is very high due to the marked variability in space and time of the endmember composition, as already reported in previous studies in snow-dominated catchments [53]. Uncertainty ranges between ±2% and ±61% (average: ±29%) and between ±2% and ±63% (average: ±31%) for scenario 1 and 2, respectively. It should be taken into account that in the method used for uncertainty calculation the relatively high standard deviation of the snowmelt and rain samples (Table 4), as often observed in high-elevation hydrological systems (e.g., [63,64]), is multiplied by the corresponding Student t-value [56]. This method produces large uncertainties, yet can capture all the variability of the precipitation conditions in the study area. As expected from the similarity between the endmembers used in the two scenarios, the calculated snowmelt fractions are highly comparable, with a maximum average difference of 0.05 (5%) for the December 2014 and May 2015 sampling campaigns. The July 2014 campaign shows the highest snowmelt contribution, with an average fraction of snowmelt of 0.94 (94%) and on average only 0.06 (6%) rain contribution. The variability of the snowmelt fraction, expressed by the box and whiskers, is also very high. This campaign was performed during an intense melting period in which the thick snow cover accumulated in the previous winter season was still melting (see Figure 5a for weather variables). During the following campaign of December 2014, the snowmelt fraction in spring water dropped significantly to approximately 30%. In this campaign spring water is characterized mainly by rain in both scenarios and the variability of the snowmelt fraction is relatively low. From a climatic point of view, the previous month of November 2014 was characterized by intense rainfall with a cumulative monthly rainfall of 594.6 mm recorded at the Col di Prà station. Moving to the following campaigns, seasonal differences are less evident. During May 2015, some springs were already subject to snow melting while others were still in baseflow conditions, indicating that this campaign shows an intermediate behavior. Even so, an increase of the snowmelt fraction can be seen compared to the previous December 2014 campaign, but the snowmelt fraction is still below the average due to the strong influence of the rain fell during the previous months. For the same reason, coupled with a scarce snow cover during the winter season, a below average snowmelt fraction is found also in July 2015. The November 2015 campaign is characterized by a decrease in the snowmelt fraction in favor of a greater rain contribution. During January 2016 all the springs were in baseflow conditions due to the prolonged dry period in the previous months. The snowmelt fraction calculated for this baseflow campaign is on average approximately 60%. July 2016 and December 2016 are characterized by a larger snow contribution and a large rain contribution, respectively. Panel (c) and (d) of Figure 6 show the same mixing model results as the panel (a) and (b) but in this case all individual springs are color coded. This representation highlights both the seasonal variations and the spatial variations of the snowmelt fraction in the different springs. For the TRA, FOC, and GAR springs, emerging at high elevations at the base of the northwester slopes of the Pale Mts., snowmelt fractions are very high, on average 75%, 71% and 68% for scenario 1, respectively. In the case of July 2014 for these three springs the mixing model was unable to represent the snowmelt fraction, giving a snow contribution above 100% in a few springs. The physically-not-meaningful calculated snowmelt fractions to spring water over 100% found in these three samples for the first sampling campaign (Figure 6c,d) indicate an inappropriate selection of the snowmelt endmembers in the mixing model. However, it must be noted that those samples were taken from the three highest elevation springs in the study area (TRA, FOC, and GAR, Figure 1) following a 30-year record snowfall period. This observation suggests that the snowmelt endmember signal, successfully used for the rest of the spring samples during the first campaign and for all samples during all other sampling campaigns, was not representative for those exceptional conditions as well. These high proportions are likely related to heterogeneous snow cover conditions of the exceptional previous winter which were not able to be thoroughly represented with just one sampling point. On the contrary, low snowmelt fractions are found in the spring emerging at the base of relatively lower peripheral ridges, in the southeastern side of the mountain group, such as TRE, NER, LUC, and SCA springs. An intermediate snowmelt fraction is found in the largest springs of the Pale Mts., ANG and PRA. These springs are characterized by the highest discharge among all the monitored springs and are located to the east and to the south of the Pale di San Martino Plateau.

5. Discussion

5.1. Variability in the Isotopic Composition of Precipitation in the Study Site

The isotopic composition in precipitation in the study area is close to the average global composition, as demonstrated by the similar slope between the LMWL and the GMWL. This trend is consistent with a prevailing continental North European Atlantic origin of the air masses originating precipitation [65]. The dual isotope plot of all samples shows that samples taken from the extended funnel gauge and from dripping snow patches have similar δ¹⁸O and δ²H values. Although the very limited number of the latter, this observation suggests that snow sampled in the extended funnel gauges represents well the bulk meltwater input, supporting the finding of Earman et al. (2006) who showed that the isotope signature from extended funnel collectors gave a robust indication of the infiltration signal [18]. Both rain and snowmelt samples plot with a comparable slope along the LMWL (Figure 3). This behavior also reflects similar precipitation origin during winter and summer. Despite numerous processes such as evaporation and sublimation, freezing and unfreezing are known to alter the isotopic composition of water from the recharge area to the springs [25,26,27] in our case we believe that these processes are of negligible importance in the study area. Indeed, evaporation and sublimation, freezing and refreezing processes would impose a shift of the dual isotopic composition of the springs compared to the meteoric water line. Our data do not show any shift, as the springs plot very close to the LMWL. Absence of evaporation is typical of Alpine settings where temperature is below 0 °C for several months of the year. It should be also considered that the carbonate aquifer under investigation is mainly formed by permeable bare karst without any soil cover, increasing water availability for recharge rather than water lost for evapotranspiration [66]. In addition, other isotope-based studies in regions close to our study area have highlighted that sublimation is not a relevant process [67] and that evapotranspiration is little for high-altitude lakes above >1500 m a.s.l., an altitude comparable to our carbonate outcrops [59]. Furthermore, in the δ¹⁸O−δ²H plot spring samples fall closer to snowmelt samples indicating that snowmelt plays a greater role in groundwater recharge compared to rain.

The isotopic signature of precipitation shows a distinct variation with elevation, an effect that is known in the literature as “altitude effect” [45,65] and is typical of mountain areas where the orographical uplift cools the air masses precipitating preferentially the heavier isotopes. The d-excess computed for precipitation shows an altitude effect as well, with increasing d-excess with elevation. Froehlich et al. (2008) suggested that the continuous increase of the d-excess from spring to late summer is associated with moisture recycling during air mass movements from the Atlantic area [68]. It is interesting to note that according to the results of these authors this effect was supposed to be limited to the northern side of the Alpine ridge, while our study site is located to the Southern side. The mean value of d-excess of +12.5‰ and +10.5‰ calculated for the RS and SM stations, respectively, is in good accordance with the average d-excess composition of all the spring samples (+12.2‰). Cervi et al. (2017) in a mountainous site which is not far from our study area, found large differences between the d-excess of groundwater and that of precipitation (in the range of +0.8‰ to +1.4‰ compared to precipitation) [67]. On the contrary, in our case, the average groundwater d-excess of +12.2‰ is comprised between the average value of precipitation recorded at the highest stations, respectively +12.5‰ for RS and +10.5‰ for SM station. These different results are likely due to the fact that we sampled precipitation along a steep elevation gradient whereas samples from the higher elevations (i.e., more representative of the recharge elevation) could not be taken in the work by Cervi et al. (2017) [67].

5.2. Seasonal Variability of the Snowmelt Contribution to Spring Water

The noticeable damping of the spring isotopic signal compared to precipitation (Table 3) indicates that newly infiltrated precipitation undergoes a mixing with pre-stored groundwater. Considering the hydrogeological setting, this storage likely occurs in the interconnected fractures of the dolomitic aquifer. Despite the presence of this well-mixed flow component, seasonal variations in the spring isotopic signature can still be recognized during the multi-annual monitoring period. According to Treck and Zojer (2009), seasonal variations of stable isotopes in groundwater are not appreciable for mean residence times exceeding five years because the isotopic signature is totally damped [69]. On the contrary a preserved seasonal amplitude of the spring water indicates a residence time shorter than one year [17,69]. In our case, the variations seem to point towards the latter option, although the significant damping of the signal might be related to residence times longer than one year. This should be confirmed by a more extensive sampling of the spring water within the study area, mainly in terms of more frequent sampling in order to better capture short-scale variability in the isotopic composition and therefore better infer hydrological and hydrogeological dynamics. The presence of clearly recognizable seasonal variations in spring water indicates that the recharge inputs have a relatively rapid influence on the groundwater storage and ultimately on spring discharge. This seasonal character of the storage highlights a climate control on the storage processes and has been recognized in other snow-dominated catchments [66,70,71]. In the case of karst, fractured aquifers a complex partitioning of the storage is due to the overlap of different flow components with different water age distributions [17]. Although water in the open fractures and conduits is rapidly drained and stored, groundwater in the poorly drained fissures moves more slowly allowing for the occurrence of mixing. The more evident seasonal cycles identified in the spring waters of the Pale Mts. are related to the input of late spring–early summer snowmelt and to the late autumn–early winter rainfall events (Figure 7). According to our results, the highest fractions of snowmelt is recognizable in the summer spring waters, contemporary to the snowmelt season (average of July campaigns 76% see Figure 7a). The fraction of snowmelt is greater after winters characterized by particularly abundant snowfalls and consequent thick snowpack. A clear example in our dataset is represented by the July 2014 sampling campaign, during which an exceptionally thick snowpack from the previous winter was still melting. In that period the snowmelt fraction in spring water was extremely high, reaching on average 94 ± 9%, increasing the hydraulic head in the saturated zone of the dolomitic massif, as observed for other carbonate aquifers [72]. The main contribution of rain to spring water recharge is concentrated in late autumn and early winter, when the snowpack is absent favoring rainfall infiltration. Indeed, the average contribution of rain to spring recharge in the autumn–early winter period is 53± 28% (Figure 7b) and spring discharge is mainly depending on the intensity of rainfall events. In case of extreme rainfall events, such as the November 2014 event, the rain contribution to spring recharge can reach almost 70% (68 ± 17%). The potential for intense precipitation to recharge the groundwater system was also observed in a recent extreme precipitation event, which hit north-eastern Italy at the end of October 2018 and that was able to sustain the baseflow during the following winter season (“Vaia” storm, see Lucianetti et al., 2019 [73]). In the winter months, recharge is prevented due to the thick snow cover (Figure 7c). These months are characterized by a drop in discharge in accordance with the relatively fast release of groundwater storage which characterizes the studied aquifer. The snowmelt fraction during the baseflow period under investigation is greater than the rain fraction (Figure 7c). This suggests that overall, during the monitoring period, snowmelt contributed more than rain to sustain baseflow. It should be considered that the last two months of 2015 in the studied region were characterized by the absence of rainfall that lasted until early 2016, resulting in one of the driest periods ever observed (meteotrentino.it). Therefore, the January 2016 campaign was able to capture a good example of baseflow conditions.

The conceptual model showed in Figure 7 offers a simplified view on the recharge and discharge processes occurring in the study area. This view is related to the selection of an inherently simple two-component mixing model assuming spring water to be only a mixture of rain and snowmelt, therefore not able to accurately represent the complex, multi-component system where rain water and snowmelt mix with sub-surface water of different age previously stored in the karst structure and potentially mobilizing under specific conditions. Multi-tracer experimental studies coupled to numerical modelling would be required to achieve a more comprehensive understanding of the mechanical behavior of this mountain group.

5.3. Spatial Variability of the Snowmelt Contribution to Spring Water

The distribution of the snowmelt fractions in spring water changes noticeably in space within the mountain group. This variability is not surprising since snow coverage and rainfall patterns are highly heterogeneous over a rugged topography. Despite this variability, the snowmelt fraction in spring water seems to follow three main spatial clusters that can be clearly recognized when comparing the map of the area with the computed snowmelt fractions.

The three spring clusters are represented with different colors in Figure 8. Cluster 1 (yellow) springs characterized by a relatively low snowmelt fraction (on average 46 ± 26%); cluster 2 (red): springs with an average snowmelt fraction of 65 ± 31%; cluster 3 (cyan): springs characterized by a high snowmelt fraction (on average 72% ± 29%). The three identified spring groups reveal that snowmelt fractions are not casually distributed but follow a spatial pattern. In particular, group 3 is located in the northwestern part of the massif, which is the most elevated of the mountain group, while the springs belonging to group 1 are in the southeastern portion of the massif, at the base of relatively low-elevated peripheral ridges. To better analyze the elevation–snowmelt fraction relationship, we explored the pattern of the computed mean recharge areas of the springs (Figure 9). The average recharge elevation of the springs can be graphically estimated connecting the spring isotopic composition with a straight line to the precipitation gradient line (see the example depicted for the PRA spring in Figure 9).

By looking at this figure the same three clusters recognized from the snowmelt fraction distribution can be identified. This indicates that the spatial variability of the snowmelt fraction is related to the elevation of the springs recharge area. Group 1 (yellow in Figure 8 and Figure 9), comprising LUC, NER, SCA, and TRE springs, has an average recharge altitude below 2000 m a.s.l. Group 3 (cyan), including TRA, GAR, and FOC springs, is characterized by mean recharge areas above 2500 m a.s.l. PRA and ANG springs (group 2, in red) have an intermediate recharge elevation of approximately 2400 m a.s.l. From this observation it can be concluded that springs with high-altitude recharge areas, above 2500 m a.s.l., are characterized by a predominance of the snowmelt fraction over the rain contribution (72% ± 29%). A greater snowmelt contribution to spring recharge is observed also for springs fed by a recharge altitude of approximately 2400 m a.s.l (65 ± 31%). On the contrary, springs recharged below 2000 m a.s.l. are recharged mostly from rain (46 ± 26%). This pattern suggests that in the same mountain group large differences in snowmelt contributions should be expected based on the springs recharge altitude. A clear relationship between the snowmelt fraction and the mean spring discharge cannot be recognized (Figure 8b). However, it is interesting to note that the two springs with the highest discharge of the Pale Mts. show a very similar snowmelt fraction and a comparable recharge elevation, despite their different geographical location. In addition, their recharge area is expected to be wide due to the elevated recharge rate. These indications suggest that the spring recharge area could be the nearby Pale di San Martino Plateau, which is a wide and nearly flat plateau covered by a thick snow cover for most of the year. Goldscheider and Holtzl (1999) highlighted the important hydrogeological role of the “plateau-like karst massifs” and differentiated them from the “folded Alpine karst” [74]. These authors indicate that plateau-like karst massifs develop deep karst systems and are often surrounded by steep slopes and major springs. In our case, the plateau portion of the Pale Mts. corresponds well to the behavior described for plateau-like karst massifs, but this interpretation needs further investigations. In any case, thanks to the three spatial recharge clusters identified in this study, the main hydrostructure of the Pale Mts. [32] can be divided into three different hydrogeological structures which can serve as a starting point for a more precise delineation of spring recharge areas. It is interesting to note that the springs belonging to group 3 (cyan in Figure 8 and Figure 9) are characterized by a high variability of the snowmelt fraction (see box plots in Figure 8). These springs also show the highest DRs, indicating low damping between the precipitation and the spring signal. This similarity between high snowmelt fractions and DRs suggests that the high variability of the isotopic signal recorded in the TRA, GAR, and FOC springs could be related to the greater snowmelt input rather than to shorter residence times. The higher standard deviation found in the snowmelt samples compared to the rain endmember (Table 3) seems to support the hypothesis of larger damping for snow-dominated springs and smaller damping for rain-dominated springs. However, it cannot be excluded that the aquifer geometry, such as the thickness of the saturated zone or the size of the groundwater storage, might influence the spring isotopic damping and thus the calculated DRs, favoring the mixing of the isotopic seasonal variability of the recharge.

6. Conclusions

This study determined the relative contribution of snowmelt and rain in the spring water of an entire mountain group providing new understanding on how this variability was distributed in space and time at the regional and seasonal scale. The use of stable isotopes proved to be a valuable tool for this purpose even in this very complex hydrogeological and topographical setting. Results of the present study can be summarized as follows:

• Snowmelt and rain are both substantial sources of recharge to the springs in this Alpine Dolomitic mountain group. This contribution is seasonally variable and it is related to two main recharge periods: the first in late spring to summer during which recharge is mainly provided by snowmelt and the second one in late autumn to early winter dominated by rain recharge. Accordingly, during the monitoring period snowmelt fractions were high during mid-summer (on average 76 ± 28%) and low during the late autumn (on average 47 ± 28%). Extraordinary meteorological events were recorded in the monitoring period and resulted in an average snowmelt fraction in groundwater of 94 ± 9% after a 30-year record snowfall winter and in a rain contribution of approximately 70% (68 ± 17%) after an intense precipitation event (386.4 mm of rain in 48 h). A greater snowmelt fraction (59 ± 40%) compared to rain was calculated during baseflow conditions following an extraordinary dry period (64 days without precipitation), indicating the predominant role of snowmelt during periods of water scarcity.
• The snowmelt contribution in spring water is spatially variable throughout the mountain group according to the elevation of the spring recharge areas. High-altitude recharge areas, above 2500 m a.s.l., are characterized by a predominance of the snowmelt fraction over the rain contribution (72% ± 29%). A greater snowmelt contribution to spring recharge is observed also for recharge altitudes of approximately 2400 m a.s.l (65 ± 31%). Springs recharged below 2000 m a.s.l. are recharged mostly from rain (snowmelt fraction of 46 ± 26%).
• We argue that the relatively rapid transmission of the recharge pulses to the spring waters reflects a strong seasonal and climate-dependent character of the storage processes. Nevertheless, the attenuation of the isotopic variations compared to precipitation does not exclude the presence of a well-mixed groundwater component with longer residence times.
• Rain-dominated springs have a more stable isotopic composition compared to snow-dominated springs and this element may be considered when using seasonal variability statistical coefficients as residence time proxies.

Climate-change predictions generally agree on a shift in the seasonality of the water cycle [75,76]. This study suggests that both rain and snowmelt are substantial sources of recharge in this Alpine Dolomitic area and that their contribution changes according to specific time windows. This seasonal partitioning of the recharge may be altered in the future affecting both water availability and groundwater-dependent ecosystems. Springs emerging from the same mountain may be more sensitive or buffer climate-change effects based on their location and in particular according to their recharge altitude. Also, during the monitoring period snowmelt had a predominant contribution both during high-flow conditions and during baseflow. This finding has implications for water availability, as during water scarcity periods the main contributor to river and stream discharge is via baseflow [77]. Also, baseflow intrusion in rivers sustains aquatic ecosystem processes and biodiversity [78]. Despite most of our results suffer from the limitation of the intrinsically simple two-component mixing model used in this work, the findings of this study offer preliminary information that can serve as a solid starting point for future hydrological and hydrogeological investigations and can be used to plan water management strategies under the current climate-change conditions.