Karst Recharge Areas Identiﬁed by Combined Application of Isotopes and Hydrogeological Budget

: The identiﬁcation of recharge areas in karst aquifers allows us to perform sustainable management of these groundwater resources. Stable isotopes ( δ 18 O and δ 2 H) have been largely used to provide information about recharge elevation in many mountainous regions. In this paper, an improved version of a recent “isotope-driven model”, for the identiﬁcation of recharge areas, was applied to Capodacqua di Spigno Spring (south of the Latium region). The model upgrade consists of a preliminary check procedure to estimate the degree of inﬂuence of the rainfall’s isotopic variability on the spring water. This additional procedure gives us an indication of the reliability of the model and its applicability conditions. Moreover, the dataset of the spring was updated to analyze the degree of reliability of the isotope-driven model. The purpose of this study was to combine the previously mentioned isotope-driven model with hydrogeological tools. A quantitative study of the basin, based on the estimation of the average monthly inﬁltration volume, was performed by using the inverse hydrogeological water budget. In this way, the qualitative model for the recharge areas’ estimation was validated by a quantitative hydrogeological tool. Both models show that, for karst mountain basins, the recharge areas decrease as the average recharge elevations increase, including areas at high altitudes.


Introduction
Karst groundwater is a primary source of freshwater in many regions around the world. About 25% of the global population is partly or entirely supplied by freshwater from karst. Most of the largest springs in the world are located in karst areas [1]. It was estimated that the percentage of karst water consumers in 2016 was 9.2% of the world's population [2]. Indeed, for a large part of the Mediterranean region, karst springs play an essential role for water supply. Karst springs provide fresh and high-quality water, which has been an important resource for human development in this region since antiquity [3]. Specifically, in Italy, the karst carbonate aquifers of the central Apennines represent the largest groundwater reservoir [3,4].
Therefore, for the protection and sustainable management of groundwater resources, it is critical to study the recharge-discharge mechanism of a karst aquifer and determine the recharge areas [5]. These aquifers are very complex systems, as each one has its own distinctive characteristics [6,7]. The high spatial heterogeneity of a karst aquifer, often characterized by poor available data, requires specific investigation techniques such as environmental tracers. At the common temperature of groundwater, there is no sensitive isotope exchange between the groundwater and minerals. Therefore, stable isotopes such as 2 H and 18 O can act as conservative tracers [8]. In this framework, the isotopic tracers play a key role. As the basic conceptual model, it is possible to claim that the isotopic composition and concentration values are affected by a number of environmental variables, which

Objectives
The present study focuses on the identification of recharge areas for karst springs by the combined application of stable isotopes (δ 18 O and δ 2 H) and hydrogeological budget. The latter is performed using the empirical method for estimating potential evapotranspiration (PET) proposed by Thornthwaite (1948) [22], applied to the inverse hydrogeological budget method (IHBM) [23][24][25]. Specifically, in the present study, the previously mentioned isotope-driven model [21] has been improved. The model extension consists of a preliminary check procedure to estimate the influence degree of seasonality combined with amount effect of rainfall (input) on spring water (output). This additional procedure aims to give an overview of the favorable conditions to the applicability of the model. Moreover, in this study, the updating of the previous dataset, regarding Capodacqua di Spigno Spring (in the southern part of the Latium region), made it possible to verify the consistency and soundness of the model.

Geological and Hydrogeological Setting of the Study Area
Capodacqua di Spigno Spring, one of the main groundwater outlets in the southeast Latium region, is sited in the province of Latina, about 4 km from the Tyrrhenian Sea [21]. This spring, located at the bottom of Petrella Mountain, has been known since Roman times. In 72 BC, the emperor Vespasian built an 11 km long aqueduct that collected spring waters to supply the city of Minturno (LT). The current aqueduct of Capodacqua di Spigno Spring feeds the larger area including Spigno Saturnia, Formia and Esperia municipalities, in Latina province [26] (Figure 1). The aquifer's altitude ranges from about 100 m a.s.l. to about 1500 m a.s.l., with an average value of 880 m a.s.l; the limits of the aquifer are assumed considering the watershed topographic pattern. The altitude of the spring, instead, is of 35 m a.s.l. This spring is located in the karst coastal region of the western Aurunci Mountains, which, together with the Lepini and the Ausoni Mountains, belong to the Pre-Apennines of Latium and form the carbonatic platform of the Volsci Ridge, separated from the Apennine ridge by the Latina Valley [21]. The spring water comes out from the permeable limestone of La Civita Mountain and Castello Mountain and flows above the Upper Miocene clays, at the lowest point of the limestone-clay contact [27]. The water table map described by previous studies [28] shows that the main direction of the aquifer is generally south-south east. Differently, there are also important local flows along the two faults, which delimit the carbonate series outcropping, especially the one to the east, which seems to represent the main water conduit towards the spring. More detailed geological information on the study area can be derived from previous studies [28]. The stratigraphic succession of the soils, from the bottom upwards, is as follows:. Deposits of the lower Pliocene (p) (continental-marine facies); • Quaternary deposits (q): river and lake deposits and groundwater debris (continental facies).
The Carbonate Series (CS1), present in the area of the Capodacqua di Spigno basin, is characterized by limestones and dolomites (Palaeocene-Upper Lias) alternations. As found in the Geological Map of Italy (scale 1:100,000, sheet no. 171 "Gaeta") [29] referring to the area under study, identified by the Civita Mountain, the aforementioned Carbonate Series (CS1) is divided into two formations: a younger Carbonate Series formed in the Cretaceous period and an older one formed during the Jurassic period. Therefore, in the geological stratigraphic succession, the Carbonate Series (CS1), Jurassic carbonate Series (JCS) and Cretaceous Carbonate Series (CCS) were taken into consideration. Figure 1 depicts a sketch of the structural setting of the Volsci chain. This scheme has been modified starting from that proposed by Baldi et al. in 2005 [28]. The new interpretation of the geology and structural characteristics of the area came from the study of both the Geological Map of Italy (scale 1:100,000, sheet no. 171 "Gaeta") [29] and the geophysical survey (tomography) performed in the area [28]. In this regard, the geological section AA' is shown in Figure 2.   Figure 1) modified from the study of the geological map of Italy [29] and geological map edited by Baldi et al., 2005 [28]: the image below depicts the structural pattern by geoelectric tomography interpretation. Legend: JCS is Jurassic Carbonate Series; CCS is Cretaceous Carbonate Series; KC is Kaolinitic Clays; Br is Clutch Breach; MS is Miocene Series and q is Quaternary deposits.   Figure 1) modified from the study of the geological map of Italy [29] and geological map edited by Baldi et al., 2005 [28]: the image below depicts the structural pattern by geoelectric tomography interpretation.   Figure 1) modified from the study of the geological map of Italy [29] and geological map edited by Baldi et al., 2005 [28]: the image below depicts the structural pattern by geoelectric tomography interpretation. Legend: JCS is Jurassic Carbonate Series; CCS is Cretaceous Carbonate Series; KC is Kaolinitic Clays; Br is Clutch Breach; MS is Miocene Series and q is Quaternary deposits.

Data Collection
In addition to the 14 samples already analyzed, whose results have been previously published [21], the isotope geochemistry laboratory of the University of Parma (Italy) tested 10 more groundwater samples, which were taken always from Capodacqua di Spigno Spring. This means that we considered 70% more samples. The IRMS (Isotoperatio mass spectrometry) continuous low-equilibration method, with CO 2 for δ 18 O and H 2 for δ 2 H, was applied. The results of all samples collected at the spring, shown in Table 1, are expressed as per mil deviations from the internationally accepted standard VSMOW (Vienna Standard Mean Ocean Water). The "¦Ä" notation is a way to express the relative abundances of oxygen ( 16 O, 18 O) and hydrogen ( 1 H, 2 H) isotopes in the water molecules [30]. The analytical error is ±0.2‰ for δ 18 O and ±1‰ for δ 2 H. During the same time period the Department of Civil, Building and Environmental Engineering (DICEA) team took spring water samples and carried out discharge measurements. Specifically, the authors made measurements along the river downward of the spring to determinate the excess flow rate of the spring, which has to be added to the discharge, collected by the local water supply agency Acqualatina Spa, in order to obtain the total amount of the spring outflow, Qs, as represented in Table 1. Measurements were conducted by the application of traditional current-meter method [21]. The stream flow velocity was measured through a SEBA horizontal axis current-meter F1, according to EN ISO 748:2007 requirements. This activity consisted mainly of measuring the excess flow rate of the spring in a discharge channel [26]. The current meter measures velocity at a one point, so it requires determination of the mean velocity in each of the selected verticals. When possible, to measure water velocity, the "Three-Point Method" described by the "National Environmental Monitoring Standards (NEMS)" was used [31].
The morphometric parameters were defined by using the digital elevation model (DEM) discretized into finite square elements (FSE), with the size of 100 × 100 m, using the open source software Q-GIS. In this way the elevation classes distribution of Capodacqua di Spigno basin was extracted from the resampled DEM.
The value of the monthly supply volume of the aquifer was calculated through the inverse hydrogeological water budget method application [23][24][25]. For the average monthly groundwater recharge estimation (i.e., the effective infiltration), data referring to four meteorological stations near the study area, for a time series of 40 years (1959-1999) for rainfall and 6 years (2012-2018) for temperature, were used. The thermo-pluviometric stations are Gaeta, Itri, Esperia and Santi Cosma e Damiano, which cover an area of about 220 km 2 ( Figure 4) [32].

Average Recharge Areas Elevation (I-Elevation)
In a given hydrogeological basin, groundwater sources may be fingerprinted using isotopes, such as 2 H and 18 O, which depend on infiltration altitude of rainfalls. Higher infiltration altitude will have precipitations, which are depleted in 2 H and 18 O. The altitude effect is used in hydrogeological studies because it allows researchers to determine at which altitude recharge areas are located. However, some studies on this topic have emphasized that the value assumed by the gradient of isotopic content varies according to meteo-climatic characteristics, too [33]. In Italy, the average values of vertical isotopic gradient are given by Grad O = −0.24‰/100 m of elevation for δ 18 O [12] and Grad H = −1.8‰/100 m of elevation for δ 2 H [12]. Therefore, in this study, the average elevation values of the recharge area (I-elevation) were estimated starting from δ 18 O and δ 2 H values detected for Sabaudia Station (17 m a.s.l.) [11]. These values are computed by solving the following formulas: Considering that ∆q = difference between the elevation value of the mean recharge area and the altitude of the Sabaudia station, it results that:

Average Recharge Areas Elevation (I-Elevation)
In a given hydrogeological basin, groundwater sources may be fingerprinted using isotopes, such as 2 H and 18 O, which depend on infiltration altitude of rainfalls. Higher infiltration altitude will have precipitations, which are depleted in 2 H and 18 O. The altitude effect is used in hydrogeological studies because it allows researchers to determine at which altitude recharge areas are located. However, some studies on this topic have emphasized that the value assumed by the gradient of isotopic content varies according to meteo-climatic characteristics, too [33]. In Italy, the average values of vertical isotopic gradient are given by Grad 18 O = −0.24‰/100 m of elevation for δ 18 O [12] and Grad 2 H = −1.8‰/100 m of elevation for δ 2 H [12]. Therefore, in this study, the average elevation values of the recharge area (I-elevation) were estimated starting from δ 18 O and δ 2 H values detected for Sabaudia Station (17 m a.s.l.) [11]. These values are computed by solving the following formulas: (1) Considering that ∆q = difference between the elevation value of the mean recharge area and the altitude of the Sabaudia station, it results that: where: q Sab is the Sabaudia station elevation (17 m a.s.l.) q2 H and q18 O are the mean elevation values of the recharge area related to 2 H and 18 O; Now it is possible to define q as the average elevation of the recharge areas between the mean elevation, respectively, related to 2 H and 18 O.

Isotope-Driven Model
The isotope-driven model proposed by Iacurto et al. [21] is based on the assumption that the isotopic composition of spring water is mainly due to the elevation effect, linked to the topographic variation of the aquifer basin surface. This assumption suggests that the variation in the isotopic concentration of rainfall, given by meteo-climatic factors, is indifferent or has little influence on the spring isotope composition. In order to demonstrate the existence of this assumption, a preliminary check procedure is introduced. The latter must be used to define the influence of meteo-climatic factors (seasonal and amount effects) on the isotopic composition variability of the spring water.

Preliminary Check of the Rainfall Influence
In the area where the model is applied, a preliminary support study, which aims to define the degree of influence of the rainfall's isotopic variability on the spring's isotopic concentration, has to be performed. A check procedure has been set up with this aim, and it assumes that the isotopic concentration of the springs changes in relation to the input rainfall ones because of temperature changes (seasonality effect) and magnitude of rainfalls (amount effect). These meteo-climatic factors act in a discordant way on isotopic concentration of rainfall [34,35]. Thus, the isotopic enrichment is due to temperatures increasing as well as to rainfall intensity and duration decreasing. On the contrary, isotopic depletion occurs when temperatures decrease and rainfall intensity and duration become bigger. In addition, based on these natural trends, intense meteoric events in summer (in the Northern Hemisphere) produce rainfalls depleted in isotopic concentration in comparison to others more commonly seen in the same season [36].
With the aim of better taking in account all these factors, any range of distribution of considered parameters has been normalized, referring to the maximum and minimum registered value. As a consequence of this, we have introduced partial relative indexes:

•
Temperature index: where ∆T 20 is the variation from 20 • C, which is usually defined as the mean condition (amount effect threshold) [36]; • Monthly rainfall index: where R m i is the monthly cumulative rain of the i-month (periodic rainfall amount), defined in terms of rain height.
• Forty eight hour rainfall index: where R 48h i is the 48 h cumulative rainfall (rainfall event amount), defined in terms of rain height.
As temperature and rainfall magnitude act discordantly on isotopic concentration of rainfall, and considering as equal the influence of monthly rainfall and short rainstorm events, a combined rainfall isotopic index, IR, can be defined in the following expression: IR can assume values in the range [−1, 1] in relation to the isotopic depletion and enrichment effects influencing rainfall water, e.g., heavy rain in winter and drought in summer, respectively, define the aforesaid extreme conditions.
In the meantime, an IS index related to the isotopic concentration of the spring is introduced to compare it with IR. Thus, IS defines the isotopic enrichment and depletion in the range [−1, 1], such as IR, and for an i-th value of isotopic concentration of the spring IS has the following expression: For any further comparative analysis, it should be considered that the isotopic concentration of groundwater at the spring is the result of the infiltration of different rainfalls which have previously involved the aquifer basin. Starting from this assumption, an analysis has been performed to highlight the relationship between the isotopic concentration trend, represented by (IS), and the combined trend of the temperature and amount effects (IR) on the isotopic concentration of multimonth rainfalls. This last trend was therefore calculated by using the average values computed for different time intervals (1, 2, 3, 4 and 5 months), of ∆T 20 i , R m i , R 48h i . The time intervals are considered as lags from the month of spring water sampling.

The Model
By using the isotope-driven model [21], it was possible to identify the main recharge areas according to the variation of isotopic rates detected in groundwater samples taken at the spring. The model is based on the ascertained linear relationship between the spring discharges and the relative average recharge area elevations, I-elevation, detected by oxygen-18 and deuterium isotopes [21]. Specifically, the linear regression ( Figure 5) suggests a main trend that associates I-elevation with an increase of the flow rate. In this way, other secondary effects, probably due to seasonality and impulsive/cumulate rainfall events influencing the isotopic concentration, are not considered in the model:

1.
With the increase in the rainfall intensity without an increase in the recharge area involved, the flow rate changes but the average recharge area elevation remains invariant; 2.
The same rainfall intensity involves different recharge areas; the flow rate remains invariant while the average recharge area elevation changes. With the aim of linking the single flow rate value to only one recharge area value, linked in turn to one value of average altitude, we introduced the equivalent geometric area (Aeq) [21]. This equivalent area is therefore calculated in terms of the weighted average of the elevation classes, as described by the following formula.
where Ai is the area related to ith-elevation hi.
In the model, it is considered that the maximum discharge of the time series is associated with the total topographic elevation distribution of the area, as wholly involved in the recharge. It is possible to assume that: The topographic average elevation is therefore defined as: where htop and heq_max are, respectively, the topographic average elevation and the average recharge area elevation associated with the maximum discharge. In addition, considering the inverse linear relationship between the average recharge area elevations and the spring discharge values, it is possible to assume that the minimal value of spring discharge is associated with a smaller equivalent recharge area (Aeq_min). This has a higher elevation (average recharge area elevation) detected by isotopic ratio values. Aeq_j is the generic equivalent recharge area associated with the generic jth-discharge, which corresponds to an average elevation hq_j, obtained by isotopic concentration values.
In this way, the analytical model identifies the recharge area associated with jth discharge as a minor quantity of the total aquifer area. The recharge area is redistributed across the real aquifer through the Weibull cumulative probability function [37], starting with covering the highest altitudes ( Figure 6). With the aim of linking the single flow rate value to only one recharge area value, linked in turn to one value of average altitude, we introduced the equivalent geometric area (A eq ) [21]. This equivalent area is therefore calculated in terms of the weighted average of the elevation classes, as described by the following formula.
where Ai is the area related to i th -elevation h i . In the model, it is considered that the maximum discharge of the time series is associated with the total topographic elevation distribution of the area, as wholly involved in the recharge. It is possible to assume that: The topographic average elevation is therefore defined as: where h top and h eq_max are, respectively, the topographic average elevation and the average recharge area elevation associated with the maximum discharge. In addition, considering the inverse linear relationship between the average recharge area elevations and the spring discharge values, it is possible to assume that the minimal value of spring discharge is associated with a smaller equivalent recharge area (A eq_min ). This has a higher elevation (average recharge area elevation) detected by isotopic ratio values. A eq_j is the generic equivalent recharge area associated with the generic j th -discharge, which corresponds to an average elevation h q_j , obtained by isotopic concentration values.
In this way, the analytical model identifies the recharge area associated with j th discharge as a minor quantity of the total aquifer area. The recharge area is redistributed across the real aquifer through the Weibull cumulative probability function [37], starting with covering the highest altitudes ( Figure 6).  (3)). The non-uniform redistribution was made using the Weibull cumulative probability function (Equation (4)), with the aim to detect the recharge area involved in the real aquifers.
The model uses the two parameters for the calibration of Weibull function (k and λ), these are found by iterative process; up to satisfy the following condition: where: A_(eq_weibull) is the weighted average area computed using the Weibull function; hi is the central value of the class elevation associated with the topographic ith-class area Ai; k and λ, respectively, are the shape and scale modelling parameters of the Weibull function. In this specific model, their behavior shows that k has a regular non-linear decrease with discharge and a regular non-linear decrease with elevation; while λ has a regular non-linear increase with elevation and a regular non-linear decrease with discharge.
As previously mentioned, the new version of the isotope-driven model includes integrations based on meteo-climatic data analysis, as shown in Figure 7.  (3)). The non-uniform redistribution was made using the Weibull cumulative probability function (Equation (4)), with the aim to detect the recharge area involved in the real aquifers.
The model uses the two parameters for the calibration of Weibull function (k and λ), these are found by iterative process; up to satisfy the following condition: where: A _(eq_weibull) is the weighted average area computed using the Weibull function; h i is the central value of the class elevation associated with the topographic i th -class area A i ; k and λ, respectively, are the shape and scale modelling parameters of the Weibull function. In this specific model, their behavior shows that k has a regular non-linear decrease with discharge and a regular non-linear decrease with elevation; while λ has a regular non-linear increase with elevation and a regular non-linear decrease with discharge.
As previously mentioned, the new version of the isotope-driven model includes integrations based on meteo-climatic data analysis, as shown in Figure 7.   Figure 8 shows the oxygen-18 and deuterium ¦Ä values for Capodacqua di Spigno Spring together with the Global Meteoric Water Line (GMWL), the Mediterranean Meteoric Water Line (MMWL) and two local meteoric water lines for central Italy, i.e., the two Central Italy Meteoric Water Lines (CIMWL) [11,12].

Oxygen-18 and Deuterium Isotopes Analyses
coastal area about 100 km away from Capodacqua di Spigno. In fact, at this spring dry months result in a greater heavy-isotope enrichment compared to the other warm months and vice versa. Therefore, the seasonal effect influences the variation in the isotopic content in the rainfall water, but is not evident in groundwater [10]. This is due to the fact that groundwater is the storage water coming from several rainfall events, which occurred over the course of several months. Additionally, with regard to the so-called amount effect, the isotopic values of the spring waters seem to have a discordant trend with what usually occurs in rainwater. In fact, in the rainy months, the rainwater is characterized by the depletion of isotopes [34][35][36], contrary to what was found in Capodacqua di Spigno Spring. The wettest months, in the case study, are in fact the winter ones, as well as the month of May 2019, which are characterized by spring waters being more enriched in heavy isotopes than during the driest months ( Figure 8). However, the main trend common to all values depends on the variation of the infiltration altitude. Actually, as expected, the altitude-topographic effect has a relevant role in the mountain areas.  The graph reported in Figure 8 highlights what was found in the first phase of monitoring and isotopic analysis of the waters of Capodacqua di Spigno Spring [21], i.e., the variation in the isotope values is connected to the variation in the spring flow rates. When Capodacqua di Spigno spring delivers maximum flow rates, δ 18 O and δD change and become enriched in 18 O and 2 H.
The altitude and (to a lesser extent) the latitude are the main geographical factors affecting the isotopic signature of precipitation in Italy [12]. Considering that in this case we are analyzing the isotopic content changes within the same study area, such a change is mainly due to the altitude-topographic effect. As for the effect of seasonality on groundwater, the enrichment in both heavy isotopes of meteoric waters is irrelevant in the warm months. In fact, for Capodacqua di Spigno Spring, in the warm months ( Figure 8) the water show lesser δ-values than in the others. In this regard Figure 8 shows that spring waters have a general δ-values trend characterized by a poor or discordant relationship with the seasonal effect characterizing the rainfall water, encountered by P. Bono et al. [38] in a coastal area about 100 km away from Capodacqua di Spigno. In fact, at this spring dry months result in a greater heavy-isotope enrichment compared to the other warm months and vice versa. Therefore, the seasonal effect influences the variation in the isotopic content in the rainfall water, but is not evident in groundwater [10]. This is due to the fact that groundwater is the storage water coming from several rainfall events, which occurred over the course of several months. Additionally, with regard to the so-called amount effect, the isotopic values of the spring waters seem to have a discordant trend with what usually occurs in rainwater. In fact, in the rainy months, the rainwater is characterized by the depletion of isotopes [34][35][36], contrary to what was found in Capodacqua di Spigno Spring. The wettest months, in the case study, are in fact the winter ones, as well as the month of May 2019, which are characterized by spring waters being more enriched in heavy isotopes than during the driest months ( Figure 8). However, the main trend common to all values depends on the variation of the infiltration altitude. Actually, as expected, the altitude-topographic effect has a relevant role in the mountain areas.
The spring water samples collected on 29 May 2019, 18 December 2018 and on 28 January 2020, regarding the events of maximum flow rate, show similarly high content of heavy isotopes. This occurred despite the different rainfall regime and season that characterized these three events. Therefore, this effect can be mainly attributed to the "altitude effect", i.e., to the infiltration of the precipitation at a lower elevation than the water collected in all the other samples. If the altitude effect is considered to be predominant, the samples that appear to be more depleted in heavy isotopes are those associated with a higher average recharge elevation. These samples were collected during spring, in the months when the spring flow rate was lower. In this regard, it is possible to affirm that the isotopic variation of the spring depends on the flow rate rather than on rainfall seasonality.

Rainfall Influence
About the case study, Figure 9 shows that the spring isotopic variability curve, represented by IS, is lowly correlated or opposed to IR that represents the variability of the meteo-climatic factors influencing the rainfall isotopic ratio. This effect is observed both considering the month of spring water sampling (1 month) and up to four previous months. In observation of the scope of the introduced check procedure, this condition is in favor of the topographic effect (altitude effect) as the main factor characterizing the water isotopic concentration of the spring. Besides, this condition meets the behavior already observed at the extremes of the series. The evidence is in the fact that the higher isotopic concentrations (isotopic enrichment) in the spring water are related to the greater rainfall magnitudes that occurred in the cold periods prior to sampling; e.g., 27 November 2019 and 18 December 2018. Different and significant is the minimum related to the 30 October 2019, in proximity to the maximum, which occurred at about the same temperature but in a period with very low rainfall. In addition, it seems that the rainstorm events (cumulate at 48 h) have a remarkable inverse correlation with the isotopic enrichment peak measured at the spring (e.g., 29 May 2019).

Recharge Areas
The spring discharge values, measured on the same days as the spring water sampling, were compared with the average elevation values of the recharge area, I-elevations, obtained using the oxygen-18 and deuterium isotopes ( Figure 10).
In December 2019, extraordinary rainfall events occurred, which on 23 December 2019 caused the river to overflow, flooding of the surrounding area ( Figure 11). It was therefore not possible to access the usual groundwater sampling point (drainage room). However, water was sampled from a different point (water fountain), where water extracted from the Capodacqua di Spigno Spring water comes directly out. This sample (red dot in Figure 10), as expected, is different from the main trend and can be considered an outlier due to the different sampling point. Indeed, in post processing analysis we realized that we sampled with the same modality performed in the drainage room. This state was not permitted to obtain an adequate purging of the pipe, considering that the water fountain is rarely open and it is quite far from the spring. With a good probability, the water stored had an isotopic concentration different from the water in the spring. Trend of the combined rainfall isotopic index (IR) and trend of the index related to the isotopic concentration of the spring (IS)-1 month = month of spring water sampling, 5 months = from the same month up to 4 months prior to sampling.

Recharge Areas
The spring discharge values, measured on the same days as the spring water sampling, were compared with the average elevation values of the recharge area, I-elevations, obtained using the oxygen-18 and deuterium isotopes ( Figure 10). In December 2019, extraordinary rainfall events occurred, which on 23 December 2019 caused the river to overflow, flooding of the surrounding area ( Figure 11). It was therefore not possible to access the usual groundwater sampling point (drainage room). However, water was sampled from a different point (water fountain), where water extracted from the Capodacqua di Spigno Spring water comes directly out. This sample (red dot in Figure 10), as expected, is different from the main trend and can be considered an outlier due to the different sampling point. Indeed, in post processing analysis we realized that we sampled with the same modality performed in the drainage room. This state was not permitted to obtain an adequate purging of the pipe, considering that the water fountain is rarely open and it is quite far from the spring. With a good probability, the water stored had an isotopic concentration different from the water in the spring. In reference to the previous values at the base of the isotope-driven model, the new isotopic vs. discharge values confirm the linear relationship between the spring discharges and average recharge area elevations detected by isotopic concentration [21]. This condition admits that the recharge area increases with the maximum flow rates of the spring; in this way their average elevation moves down, and vice versa. In addition, the assumption at the base of the model is that the recharge area associated with the maximum flow rate involves the whole area with higher elevation than the spring point. In this regard, as written in Equation (13), the topographic average elevation is set equal to the average recharge area elevation associated with the maximum discharge. The average recharge area elevation associated with the maximum discharge in May 2019 is similar to that of December 2018. Differently, the maximum discharge occurs in November 2019 gives a lesser I-elevation than this latter. This condition, which was probably effected by extraordinary 48 h rainfall event (Figure 9), was caused by a great runoff on the basin surface that shifted the infiltration at lower elevations. On the contrary, rainfall events more distributed over time are not characterized by great runoff on the topographic surface and therefore do not lead to a lowering of the average recharge area elevation. How- In reference to the previous values at the base of the isotope-driven model, the new isotopic vs. discharge values confirm the linear relationship between the spring discharges and average recharge area elevations detected by isotopic concentration [21]. This condition admits that the recharge area increases with the maximum flow rates of the spring; in this way their average elevation moves down, and vice versa. In addition, the assumption at the base of the model is that the recharge area associated with the maximum flow rate involves the whole area with higher elevation than the spring point. In this regard, as written in Equation (13), the topographic average elevation is set equal to the average recharge area elevation associated with the maximum discharge. The average recharge area elevation associated with the maximum discharge in May 2019 is similar to that of December 2018. Differently, the maximum discharge occurs in November 2019 gives a lesser I-elevation than this latter. This condition, which was probably effected by extraordinary 48 h rainfall event ( Figure 9), was caused by a great runoff on the basin surface that shifted the infiltration at lower elevations. On the contrary, rainfall events more distributed over time are not characterized by great runoff on the topographic surface and therefore do not lead to a lowering of the average recharge area elevation. However, this lesser I-elevation does not produce substantial change in the slope of the I-elevation vs. discharge regression line ( Figure 10). This condition permits us to use the isotopic driven model proposed in the previous study [21].
Starting from the elevation classes distribution of the topographic surface of Capodacqua di Spigno Spring basin, it was possible to apply the isotope-driven model to simulate the spring recharge areas. Figure 12 shows the recharge area distribution related to the maximum discharge (total area of the basin) and the recharge areas' distributions related to the two minimum discharges of 2018 and 2019, simulated by a model with a Weibull function.

Estimation of the Aquifer Monthly Supply Volume
With the aim of having a hydrological validation of the isotope-driven model, a comparison was made between the estimated recharge areas, with the aforesaid model, and the monthly variability of the infiltration volume. This comparison was possible in that the monthly average recharge (monthly infiltration amount) of the Capodacqua di Spigno basin was distributed across the areas identified by different elevation classes.
The monthly average infiltration in the Capodacqua di Spigno basin has been evaluated, using the empirical method for estimating potential evapotranspiration (PET) proposed by Thornthwaite (1948) [22] applied to the inverse hydrogeological budget method (IHBM) [23][24][25]. The equation only requires mean monthly air temperature and mean daily daylight hours for each month, which can be calculated based on the latitude at the site of interest. With the aim of estimating the average monthly recharge of Capodacqua di Spigno basin, it was necessary to calculate the real monthly evaporation (RET), starting from potential monthly evapotranspiration (PET). For each cell, we compared the rainfall, through the rainfall-altitude relationships, with the PET. The real monthly evaporation was made equal to potential monthly evaporation if the rainfall value was bigger of the PET value, while the real monthly evaporation was imposed as equal to rainfall value if the rainfall value was less than the PET value.
First, the annual heat index (i), for each of the four rainfall stations (Figure 4) was calculated, so in the calculation of the potential evapotranspiration (PET) the average of the four annual heat index values was used (Imean). In this regard, a summary table of the calculated values, necessary to estimate the PET of Thornthwaite, for each cell of the grid, is reported below ( Table 2). Table 2. Monthly Thornthwaite heat index (i), annual heat index (I), annual mean heat Index (Imean) and α coefficient for Capodacqua di Spigno basin.

Estimation of the Aquifer Monthly Supply Volume
With the aim of having a hydrological validation of the isotope-driven model, a comparison was made between the estimated recharge areas, with the aforesaid model, and the monthly variability of the infiltration volume. This comparison was possible in that the monthly average recharge (monthly infiltration amount) of the Capodacqua di Spigno basin was distributed across the areas identified by different elevation classes.
The monthly average infiltration in the Capodacqua di Spigno basin has been evaluated, using the empirical method for estimating potential evapotranspiration (PET) proposed by Thornthwaite (1948) [22] applied to the inverse hydrogeological budget method (IHBM) [23][24][25]. The equation only requires mean monthly air temperature and mean daily daylight hours for each month, which can be calculated based on the latitude at the site of interest. With the aim of estimating the average monthly recharge of Capodacqua di Spigno basin, it was necessary to calculate the real monthly evaporation (RET), starting from potential monthly evapotranspiration (PET). For each cell, we compared the rainfall, through the rainfall-altitude relationships, with the PET. The real monthly evaporation was made equal to potential monthly evaporation if the rainfall value was bigger of the PET value, while the real monthly evaporation was imposed as equal to rainfall value if the rainfall value was less than the PET value.
First, the annual heat index (i), for each of the four rainfall stations (Figure 4) was calculated, so in the calculation of the potential evapotranspiration (PET) the average of the four annual heat index values was used (I mean ). In this regard, a summary table of the calculated values, necessary to estimate the PET of Thornthwaite, for each cell of the grid, is reported below ( Table 2).   (Figure 13). For the rainfall-altitude and temperaturealtitude relationships, linear interpolation has usually been used [24]. However, in the case of mountain basins, a logarithmic or exponential rainfall-altitude relationship leads us to obtain much more plausible values for the higher altitudes [32,39].
In this way, it was possible to estimate, for each finite square element (FSE), the value of the monthly potential evaporation (PET) and the real monthly evaporation (RET).
By multiplying the single FSE infiltration value by the value of its area it was possible to obtain the average monthly infiltration volume (monthly average recharge). Figure 14 shows the estimation of the aquifer's monthly average recharge (M.A.R.) for each month of the year. Figure 15 shows the distribution of the average recharge volume (M.A.R.) in the different months of the year. In Figure 15 it is evident that November (line blue) is the month with the highest infiltrated volume. In this month the recharge volume increases with the altitude, according to the rainfall-altitude relationship, based on the area available at each altitude. On the other hand, in May (line red), which, except for months with no infiltration, is the month with lowest infiltrated volume, only part of the basin is involved in the infiltration process. At low altitudes, in fact, the infiltration volume is zero. This is in accordance with the findings of the model, i.e., with the estimate of the recharge areas' distribution associated with the minimum and maximum flow rates ( Figure 12).
Regarding Capodacqua di Spigno Spring, in a previous study the cross-correlation coefficient (C PiQ ) between the monthly spring discharge series and monthly rainfall series was analyzed [26]. C PiQ numerically indicates how much each single previous monthly rainfall value influences the total spring flow rate in relation to the specific month considered. For Capodacqua di Spigno Spring, the highest values of the correlation coefficients are from the same month up to four previous months (from i = 0 to i = 3) with respect to the spring flow rate occurrence. Considering that the minimum flow rate of the Capodacqua di Spigno Spring usually occurs between the end of August and mid-September, it can be assumed that the minimum discharge is mainly due to the infiltrations that occurred in the months of May, June, July and August. Therefore, it is possible to state that the minimum flow rate of the spring comes mainly from the infiltration that takes place in May, since the infiltrations of the months of June, July and August are equal to zero ( Figure 14). ues (A.M.R.M. and A.M.T.) it was possible to obtain the rainfall-altitude and temperature-altitude relationships for each month (Figure 13). For the rainfall-altitude and temperature-altitude relationships, linear interpolation has usually been used [24]. However, in the case of mountain basins, a logarithmic or exponential rainfall-altitude relationship leads us to obtain much more plausible values for the higher altitudes [32,39].
In this way, it was possible to estimate, for each finite square element (FSE), the value of the monthly potential evaporation (PET) and the real monthly evaporation (RET). By multiplying the single FSE infiltration value by the value of its area it was possible to obtain the average monthly infiltration volume (monthly average recharge). Figure 14 Figure 15 it is evident that November (line blue) is the month with the highest infiltrated volume. In this month the recharge volume increases with the altitude, according to the rainfall-altitude relationship, based on the area available at each altitude. On the other hand, in May (line red), which, except for months with no infiltration, is the month with lowest infiltrated volume, only part of the basin is involved in the infiltration process. At low altitudes, in fact, the infiltration volume is zero. This is in accordance with the findings of the model, i.e., with the estimate of the recharge areas' distribution associated with the minimum and maximum flow rates ( Figure 12). with the highest infiltrated volume. In this month the recharge volume increases with the altitude, according to the rainfall-altitude relationship, based on the area available at each altitude. On the other hand, in May (line red), which, except for months with no infiltration, is the month with lowest infiltrated volume, only part of the basin is involved in the infiltration process. At low altitudes, in fact, the infiltration volume is zero. This is in accordance with the findings of the model, i.e., with the estimate of the recharge areas' distribution associated with the minimum and maximum flow rates ( Figure 12). Regarding Capodacqua di Spigno Spring, in a previous study the cross-correlation coefficient (CPiQ) between the monthly spring discharge series and monthly rainfall series was analyzed [26]. CPiQ numerically indicates how much each single previous monthly We therefore confirm what was simulated by the isotope driven model. It means that the recharge area associated with the minimum discharge is smaller and includes areas at higher altitudes. As a matter of fact, the graph in Figure 16 shows how the infiltration volume increases according to recharge area, since all the altitudes of the basin are affected by the infiltration.
Water 2021, 13,1965 22 of 26 rainfall value influences the total spring flow rate in relation to the specific month considered. For Capodacqua di Spigno Spring, the highest values of the correlation coefficients are from the same month up to four previous months (from i = 0 to i = 3) with respect to the spring flow rate occurrence. Considering that the minimum flow rate of the Capodacqua di Spigno Spring usually occurs between the end of August and mid-September, it can be assumed that the minimum discharge is mainly due to the infiltrations that occurred in the months of May, June, July and August. Therefore, it is possible to state that the minimum flow rate of the spring comes mainly from the infiltration that takes place in May, since the infiltrations of the months of June, July and August are equal to zero ( Figure 14).
We therefore confirm what was simulated by the isotope driven model. It means that the recharge area associated with the minimum discharge is smaller and includes areas at higher altitudes. As a matter of fact, the graph in Figure 16 shows how the infiltration volume increases according to recharge area, since all the altitudes of the basin are affected by the infiltration.

Isotope-Driven Model vs. Hydrogeological Budget Applications
Using an open source GIS plugin (GEarthView) it was possible to represent the results described above on Google Earth with the images presented in Figure 17. These figures show the simulated area involved by the infiltrations of May, identified by the inverse hydrogeological balance, and the recharge area related to Qmin carried out by the isotope driven model, respectively. It is quite evident that the two areas are very similar:

Isotope-Driven Model vs. Hydrogeological Budget Applications
Using an open source GIS plugin (GEarthView) it was possible to represent the results described above on Google Earth with the images presented in Figure 17. These figures show the simulated area involved by the infiltrations of May, identified by the inverse hydrogeological balance, and the recharge area related to Q min carried out by the isotope driven model, respectively. It is quite evident that the two areas are very similar: they are smaller than the total surface of the Capodacqua di Spigno basin and involve only the highest altitudes.

Discussion and Conclusions
The isotope-driven model, set up with the aim of reaching a reliable means of id tifying aquifer recharge areas, turned out to be a powerful tool, as this target is very portant in several studies addressing groundwater protection. In this paper the isot set up model presented in a previous paper [21] has been improved with a prelimin check procedure. This procedure is an important supplement to the model. This aim compare the trend of the possible natural factors that influence the isotopic concentra of rainfall with the trend of the spring's isotopic concentration. In relation to these tren a discordance or absence of correlation permits us to assume the elevation effect as dominant in the spring waters, and so the "isotope-driven model" can be run.
Several of the isotope-based techniques involve interpreting differences between isotopic compositions of rainfall and groundwater samples [10,40,41]. However, in paper, the preliminary check procedure analyses how the single affecting factor (seaso and amount effect) that leads the variability of the rainfall's isotopic composition can fluence the isotopic composition of the spring waters. In karst mountain aquifers, it is w known that the altitude effect is predominant compared to the other effects on sprin isotopic composition, and allows us to define the average recharge altitude [12,17, However, this does not provide an indication of the recharge area. In this regard, the p

Discussion and Conclusions
The isotope-driven model, set up with the aim of reaching a reliable means of identifying aquifer recharge areas, turned out to be a powerful tool, as this target is very important in several studies addressing groundwater protection. In this paper the isotope set up model presented in a previous paper [21] has been improved with a preliminary check procedure. This procedure is an important supplement to the model. This aims to compare the trend of the possible natural factors that influence the isotopic concentration of rainfall with the trend of the spring's isotopic concentration. In relation to these trends, a discordance or absence of correlation permits us to assume the elevation effect as predominant in the spring waters, and so the "isotope-driven model" can be run.
Several of the isotope-based techniques involve interpreting differences between the isotopic compositions of rainfall and groundwater samples [10,40,41]. However, in this paper, the preliminary check procedure analyses how the single affecting factor (seasonal and amount effect) that leads the variability of the rainfall's isotopic composition can influence the isotopic composition of the spring waters. In karst mountain aquifers, it is well known that the altitude effect is predominant compared to the other effects on spring's isotopic composition, and allows us to define the average recharge altitude [12,17,42]. However, this does not provide an indication of the recharge area. In this regard, the proposed model estimates the main recharge areas that feed the different spring discharges during the hydrological year.
The model has also been strengthened by the updating of new 2 H and 18 O data, with respect to the previous one [21], so it was possible to verify its soundness. Indeed, the 10 more spring water samples have made it possible to apply the model on a more extended isotopic dataset. This allowed us to analyze several maximum and minimum annual flow rates conditions.
On the other hand, as the application of tracers to groundwater characterization is generally preferable in a multitracing approach framework, in this case the set up model has been validated by a quantitative hydrogeological tool. It means that it has applied the Inverse Hydrogeological Budget method, modified for a monthly scale application, with the aim of finding the distribution of the mean infiltration volume in each area. Specifically, the monthly infiltration amount was distributed on the areas identified by different elevation classes. In this way, the minimum recharge areas identified with the isotope-driven model were compared with the monthly infiltration areas of those months related to the minimum spring flow rate. The research highlights that the distribution of the recharge areas on elevation classes given by hydrological water balance is very similar to the one resulting from the application of the isotope-driven model. Specifically, the application of the hydrogeological balance has confirmed that the total area of the basin is not always involved in the recharge process throughout the hydrological year. In dry periods, there is no infiltration in low altitude areas. In fact, considering May-August, i.e., the 4 months preceding the minimum flow rate of the spring, the corresponding recharge area is not the whole basin. This leads, in mountain karst basins such as the Capodacqua di Spigno one, to an increase in the average recharge altitude, which affects the isotopic values found in the spring water. The proposed model, which uses groundwater's isotopic signature, can be a preliminary tool for estimating the infiltration area and, therefore, improving the protection management of the aquifer, in reference to possible contamination sources affecting the feeding area.