An Attempt to Characterize the Recharge of Alluvial Fans Facing the Northern Italian Apennines: Indications from Water Stable Isotopes

Nowadays, climate changes and increased water demand for human and agricultural purposes pose important questions for the groundwater management of alluvial aquifers facing the northern Italian Apennines. The large groundwater withdrawals, coupled with an overall worsening of the water quality, requires a detailed knowledge of the recharge mechanisms of these aquifers that can be useful for further adaptation measures. Concerning the recharge area of the alluvial aquifers (i.e., apices made up of gravelly materials), the present study investigates a dataset made up of 282 water samples for which stable isotopes oxygen-18 (18O) and deuterium (2H) are available. The latter involves precipitations (three rain gauges), surface water (five rivers) and groundwater (twenty wells) from five selected alluvial fans. The study confirms that the different isotopic signatures characterizing rain and river water from this area can be exploited for preliminary characterization of their significance on groundwater recharge. These results lay the foundations for the further use of a suite of environmental tracers (in which a primary role is that of water stable isotopes) at the event-scale (i.e., that of rainfall and/or flood) for eventually estimating the effective quota of recharge linked to precipitation and surface water.


Introduction
Basin aquifers worldwide are located in the vicinity of mountain chains and are usually made up of porous materials. Being located in lowland areas, i.e., where urbanization phenomena and population are concentrated, they are often intensively exploited for water supply [1]. It is well recognized that a remarkable portion of their recharge comes from the neighbouring mountain chains. In fact, due to the orographic effect, mountain areas are usually characterized by greater amounts of precipitation that are subsequently released to the lowlands representing, for the latter, an additional water quota with respect to the only zenithal precipitation for their recharge. All the processes acting to recharge the basin aquifer can be split into two sub-processes, whose intensities may vary in time and space even considering the same basin aquifer [2]. In detail, the first sub-process is represented by the sub-surface flow from the adjacent mountains (hereafter called Mountain-Block Recharge, MBR), Water 2020, 12, 1561 2 of 21 and includes all the groundwater entering the basin aquifer from the buried sector of the mountain block. This usually consists of water that travels in long and deep flow paths. The second sub-process, named Mountain-Front Recharge (MFR), is related to the streambed dispersion of mountain-sourced perennial and ephemeral rivers after these streams enter the lowlands. This process is particularly intense where the streams enter the lowland areas, i.e., where they flow over their alluvial fans. The latter are prominent sedimentary landforms made up of unconsolidated and coarse deposits released by the river as a result of a sudden change in flow velocities driven by slopes that pass from a high to a low gradient. It should be pointed out that grain sizes of the porous material composing the riverbeds along with their alluvial fans progressively decrease moving downslope; for this reason, the first part of an alluvial fan (i.e., apical part) shows coarser materials which in turn lead to higher permeability [3,4]. This means that infiltration of surface water is promoted in the apical parts of the alluvial fans allowing remarkable quotas of river flows to feed the underlying basin aquifers (up to being, in the case of arid and semi-arid regions, the dominant source of recharge to lowland aquifers; [5]). Characterization and estimation of MBR and MFR is a complex task that can be addressed with different methods (see, for instance, the recent review by [6], and the former one made by [7]). Among the others, tracing techniques based on many environmental tracers have been widely used to provide both qualitative and quantitative estimates of the two sub-processes. As an example, and without claiming to be exhaustive, environmental tracers such as noble gases [8], radioactive isotopes (e.g., [9]) and stable isotopes (e.g., [10]) have been successful in determining which one between MFR or MBR is the dominant recharge mechanism. Multi-tracer approaches exploiting several environmental tracers have proven to be effective in determining the quotas of MFR and MBR in Switzerland (noble gases, stable isotopes, and radioactive isotopes; [11]), Spain (radioactive isotopes and noble gases; [12]), and the USA (radioactive isotopes and noble gases; [13]). In all cases, from an operative point of view, it would be useful to collect water samples from all the components participating at the hydrological cycle (i.e. precipitation, surface water, groundwater) to quantify the content of tracers characterizing each component. The more the tracer content of each component is different from the other, the more the methodology turns out to be effective. Thus, it is evident that, although a suite of environmental tracers would be preferable for the amount of information that can be inferred, their use remains expensive and their interpretation difficult. Among the several available environmental tracers, oxygen ( 18 O and 16 O) and hydrogen isotopes ( 2 H and 1 H) have proven to be a powerful tool for depicting interaction between surface water and groundwater along rivers [14]. In the literature, they have allowed the determination of the alpine versus land surface recharge sources for groundwater [15,16] and to characterize and quantify the several sources of recharge in aquifers hosted in alluvial fans [17][18][19].
In Italy, the alluvial fans of rivers flowing out from the northern Apennines and Alps have always been considered as the main recharge zones for aquifers hosted in the Po Plain [20]. Concerning the southern sector of the Po Plain (i.e., the sector facing the northern Apennines), the MBR component cannot provide recharge quotas as the mountain chain bedrock is made up of impermeable clayey materials. Here, the MFR recharge is focused in the apical part of the alluvial fans [21] but, to our knowledge, it has never been quantified except with the exclusive use of modelling approaches for small size test areas [22][23][24]. Starting from the Seventies, significant land-use changes in the mountainous sector of catchments, dam constructions, river diversion developments, torrent control works, and extensive bed material mining have caused important sediment budget changes which in turn have led to significant alterations of the channel morphologies in the apical part of the alluvial fans. Here, reduction and deepening (even of tens of meters) of the streambeds have been observed for almost all the rivers, with the consequent change of river morphology from braided to single channel [25,26]. In some cases, streambeds currently lie over the impermeable formations composing the foothills of the northern Apennines [27,28]. This may have caused a deficit in the recharge of the aforementioned aquifers as a consequence of the reduced water losses through the streambeds, which has already been speculated by other authors [28]. The slight reduction of precipitations detected in the upper parts of the catchments [29], as well as land-use change (i.e., catchments re-naturalization with a subsequent increase in evapotranspiration processes), have led to a reduction in mean annual river flow rates [26].
At the same time, decreases in groundwater levels from the apical parts of the majority of the alluvial fans in the northern Italian Apennines have already been reported by [30].
Understanding and quantifying the components (i.e., precipitation and/or river streambed dispersion) of recharge for the aquifers hosted in the alluvial fans from the northern Italian Apennines is a challenge. In fact, and due to their almost impermeable catchments, here rivers respond quickly to rainfall (a few days/hours during the wettest periods), not allowing for a clear differentiation of the recharging rates by precipitation and riverbed dispersion.
In the light of this framework, in this paper we examine a dataset consisting of grabbed water samples collected in rivers, wells, and rain gauges from the northern Italian Apennines. The three-year dataset is made up of monthly (rivers and precipitations) and four-monthly (wells) values of stable isotopes of oxygen ( 18 O and 16 O) and hydrogen ( 2 H and 1 H) coupled with the monthly discharge data of rivers and the monthly amounts of precipitations.
Here, a previous study [31] had already highlighted a marked discrepancy between stable isotopic values in precipitations taking place in the alluvial fans and those in surface water from rivers. This leads these isotopes to be potentially useful tracers in an attempt to understand the current role of the dispersion from the riverbed in the recharging of the aquifers.

Climatic and Hydrological Setting
The study area is located at the foothills of the northern Apennines of Italy and includes the apices of five water courses, namely, from northwest to southeast, Trebbia, Taro, Enza, Secchia and Reno ( Figure 1). Apices starts from the inside of the mountain chain and develop towards NE forming the alluvial fans (see for more detail Section 2.2 Hydrogeological setting). As a result, the elevation of the apices decreases in the NE direction and ranges from about 250 m a.s.l. to approximately 50 m a.s.l.
As reported in [29], the mean annual rainfall distribution over the period 1990-2015 is close to 900 mm/y (with the slightly higher value in the Trebbia apical of 950 mm/year). It should be pointed out that the rainfall distribution during the year is characterized by two marked positive peaks in autumn and springs seasons, when cumulative amounts may account for (400 mm) and (500 mm), respectively (in the summer months precipitations are reduced as well as during the winter ones). Moreover, potential evapotranspiration is particularly active during the summer period, when it reaches cumulative values of up to 600 mm [29].
All the rivers originate from the main watershed divide with altitudes that, in the case of the Secchia river, can exceed 2000 m a.s.l. Thus, the Secchia river is characterized by the nival-pluvial discharge regime as it is influenced by the melting of snow cover accumulated during the winter months in the upper parts of its catchments (where cumulative annual snow cover can reach up to 2-3 m). The Trebbia, the Taro, the Enza and the Reno rivers show a marked pluvial discharge regime.
More in detail, and by considering the period 2003-2017 [32], all the five rivers show perennial discharges with mean annual values at their apices between 9 m 3 /s (the Enza river) and 40 m 3 /s (the Taro river). In all cases, the flow rates vary widely during the year because of the poorly permeable bedrocks widely outcropping in the mountainous part of their catchment [33]. As a result, river discharges closely follow the rainfall distribution during the year. This means that low-flows take place in the summer-early autumn periods (August, September, and October) with minimum values of discharges that can be in the order of 2 m 3 /s (the Enza river). Floods usually occur in autumn (October and November) and spring (March and April), with discharge peaks that can be higher than 700 m 3 /s (the Secchia, the Taro, and the Reno rivers). It is worth noting that, in these two wet periods, the time lag of peak discharges after the most extreme precipitation events is usually less than 24 h. Sketch map of the five alluvial fans considered in this study and related to rivers Trebbia, Taro, Enza, Secchia and Reno along with locations in which groundwater (wells with roman letters from I to XX), precipitation (rain gauges numbered from 1 to 3) and surficial (rivers with letters from a to e) water samplings have been carried out by previous study.

Hydrogeological Setting
Alluvial fans from the northern Apennines of Italy are oriented S-N and have developed over a basal aquiclude made up of marine clays deposited during the Pliocene (see Figure 2). The alluvial fans are composed of sedimentary particles (grain size ranging from clays to gravels) that were made available by the continuous dismembering of the mountain chain substrata. Since the long-term variations in the discharge of rivers that occurred during the Pleistocene have implied a change in their solid transport, several lenses of gravelly deposits have propagated northwards (distal parts can be located up to 20-25 km from the foothills) in the most humid periods. On the contrary, in the driest periods, the coarsest deposits have stopped in the vicinity of the foothills of the northern Apennines and were followed, toward the north, by finer deposits made up of silty-clayey materials ( [34]; Figure 2).
Water 2020, 12, x FOR PEER REVIEW 4 of 23 from I to XX), precipitation (rain gauges numbered from 1 to 3) and surficial (rivers with letters from a to e) water samplings have been carried out by previous study.
As reported in [29], the mean annual rainfall distribution over the period 1990-2015 is close to 900 mm/y (with the slightly higher value in the Trebbia apical of 950 mm/year). It should be pointed out that the rainfall distribution during the year is characterized by two marked positive peaks in autumn and springs seasons, when cumulative amounts may account for (400 mm) and (500 mm), respectively (in the summer months precipitations are reduced as well as during the winter ones). Moreover, potential evapotranspiration is particularly active during the summer period, when it reaches cumulative values of up to 600 mm [29].
All the rivers originate from the main watershed divide with altitudes that, in the case of the Secchia river, can exceed 2000 m a.s.l. Thus, the Secchia river is characterized by the nival-pluvial discharge regime as it is influenced by the melting of snow cover accumulated during the winter months in the upper parts of its catchments (where cumulative annual snow cover can reach up to 2-3 m). The Trebbia, the Taro, the Enza and the Reno rivers show a marked pluvial discharge regime.
More in detail, and by considering the period 2003-2017 [32], all the five rivers show perennial discharges with mean annual values at their apices between 9 m 3 /s (the Enza river) and 40 m 3 /s (the Taro river). In all cases, the flow rates vary widely during the year because of the poorly permeable bedrocks widely outcropping in the mountainous part of their catchment [33]. As a result, river discharges closely follow the rainfall distribution during the year. This means that low-flows take place in the summer-early autumn periods (August, September, and October) with minimum values of discharges that can be in the order of 2 m 3 /s (the Enza river). Floods usually occur in autumn (October and November) and spring (March and April), with discharge peaks that can be higher than 700 m 3 /s (the Secchia, the Taro, and the Reno rivers). It is worth noting that, in these two wet periods, the time lag of peak discharges after the most extreme precipitation events is usually less than 24 h.

Hydrogeological Setting
Alluvial fans from the northern Apennines of Italy are oriented S-N and have developed over a basal aquiclude made up of marine clays deposited during the Pliocene (see Figure 2). The alluvial fans are composed of sedimentary particles (grain size ranging from clays to gravels) that were made available by the continuous dismembering of the mountain chain substrata. Since the long-term variations in the discharge of rivers that occurred during the Pleistocene have implied a change in their solid transport, several lenses of gravelly deposits have propagated northwards (distal parts can be located up to 20-25 km from the foothills) in the most humid periods. On the contrary, in the driest periods, the coarsest deposits have stopped in the vicinity of the foothills of the northern Apennines and were followed, toward the north, by finer deposits made up of silty-clayey materials ([34]; Figure  2).

Figure 2.
Simplified hydrogeological cross-section of a representative alluvial fan from the northern Apennines of Italy (modified after [34]).

of 21
As a result, each alluvial fan is characterized, in the vicinity of the mountain foothills, by an almost undifferentiated gravelly pack (apical part) that, by moving northward, evolves into an overlay of poorly permeable (aquitards and aquicludes composed of clayey and silty materials) and permeable materials (aquifers made up of gravels and sands). Although the stratigraphy of each alluvial fan may be different due to specific tectonic conditions (such as the presence/absence of buried thrusts), their thickness usually increases towards the north from some tens of meters in the apical part up to more than 700 m in the lowlands [35][36][37].
Due to the abovementioned hydrogeological features, the apical part of each alluvial fan acts as a monolayer aquifer hosting unique groundwater, whose flow-paths slowly move towards the north. Its recharge is supplied both by rainfall percolation from the topsoil and the streambed dispersions (i.e., MFR); as a result, groundwater levels change slightly during the years (up to more than 5 m; see Figures S1 and S2 in Supplementary Materials for comparison of precipitations, surface water and groundwater from the Taro and the Enza rivers, respectively). It is worth noting that here Pliocene clays composing the foothills of the northern Apennines do not allow the development of flow-paths; thus MBR is not promoted. Downstream of the apical part, aquifers become multi-compartmentalized as they are superimposed to aquitards and aquicludes; groundwater levels change little during the year.
It should be pointed out that the dynamics of water consumption (due to both water supply for human purposes and agriculture) have caused decreases in the piezometric levels over the whole alluvial fans, which are more pronounced in the deep and confined aquifers from the distal parts. In the apical parts of the fans, a number of hydrochemical features also indicate direct and fast recharges. Firstly, wells from the apical parts present a high content of nitrates driven by intensive agriculture practices ( [38,39]). Detailed isotopic analysis involving both nitrogen and boron have demonstrated that high nitrate concentrations are correlated here with pig farming since manure spreading represents one of the primary nitrate sources in groundwater from agriculture, the other being synthetic fertilizers [40]. Moreover, nitrate contents vary during the year and are strongly dependent on rainfall patterns ( [21]; see also Figure S3 in Supplementary Materials for nitrates contents in three wells from the Enza river). Secondly, as in the case of nitrates, total dissolved solids also change during the year (Ca (Mg)-HCO 3 facies with Total Dissolved Solids TDS ranging from 200 to 1500 mg/l) with decrement concentrated in the recharge period (from November to April; [21]). Thirdly, Tritium values in groundwater are in the order of those noticed in rains, indicating short transit times of percolating water within the unsaturated zones (lower than one year; [41]). Fourthly, the redox potential (Eh) in groundwater from the apical parts is always high (up to 350 mV), indicating the very recent exposure of these water to the atmospheric oxygen [41].

Isotopic Data in Precipitation, Rivers, and Groundwater
Stable oxygen and hydrogen isotopic data from precipitation (3 rain-gauges), surface water (5 rivers) and groundwater (20 water wells) are derived from [31] and [41]. The final dataset consists of 282 isotopic values, of which 61 from precipitation, 144 from surface water, and 77 from groundwater. The three rain-gauges (namely, Parma, Lodesana and Langhirano) are located in the proximity of alluvial fan of the Taro River. Following [42], such rain gauges can be considered as representative of the isotopic features of precipitation in the whole study area (see Section 3.2 "Reference Isotopic Compositions in Precipitations from the foothills of the northern Italian Apennines"). Being positioned at higher elevation, the rain gauge from Langhirano (220 m a.s.l,) can be ascribed to the apical area of the alluvial fans while Lodesana (150 m a.s.l.) and, more specifically, Parma (65 m a.s.l.) have to be considered as characterizing the more distal part of the alluvial fans.
The five river gauges in which surface water samples for isotopic analyses were collected are located in the apical parts of the corresponding alluvial fans (Figure 1) as well as groundwater from 20 water wells. The latter are selected as characterized by both low-depth groundwater levels close to the river channels, and are distributed as follows: 5 from Trebbia, 4 from Taro, 3 from Enza, 5 from Secchia, 3 from Reno alluvial fans.
The dataset from rain gauges and rivers is represented from monthly data, and water wells were characterized by grabbed four-monthly samples. With reference to rain gauges, isotopic datasets last over the period from January 2003 to December 2006. Isotopic data from surface water and groundwater cover the period from January 2004 to December 2007.
As reported by [31] and [41], the isotopic analyses were carried out by using Isotope Ratio Mass Spectrometry (IRMS) and are reported as a deviation of the sample from the standard (Vienna Standard Mean Oceanic Water: V-SMOW). They are presented in the δ-notation as per mil (% ), where δ =  [43]), which is calculated as: Corresponding uncertainty on d E are assessed as ±0.7% . All the data are reported in the form of Supplementary Material.

Reference Isotopic Compositions in Precipitations from the Foothills of the Northern Italian Apennines
Many authors have dealt with the isotopic composition of precipitations from the foothills of the northern Apennines of Italy [42,[44][45][46]. All the authors have confirmed the temperature-linked variability of δ 18 O (and in turn δ 2 H), with lower values in winter than summer months and a final isotopic signal from each rain gauge that commonly approximates a sinusoid. δ 18 O-δ 2 H pairs from precipitations are usually aligned along a regression line everywhere. If all the world-scale rain gauges are considered, δ 18 O and δ 2 H pairs lie along the regression line called the "Global Meteoric Water Line": GMWL [47], which is defined as: Although this relation is valid everywhere, more accurate regression lines can be obtained by selecting δ 18 O and δ 2 H values from restricted areas. This is due to the specific isotopic fractionation processes (i.e., vapor pressure and temperature conditions) controlling the precipitation over each area. If rain gauges from the Mediterranean area alone are considered, a regression line comes out with a slightly higher intercept than that of the GWML (the Mediterranean Meteoric Water Line MMWL; [48]): Table 1 summarizes the Local Meteoric Water Lines (LMWLs) obtained by the abovementioned authors for the study area. Slopes of the alignments are close to that of GMWL and MMWL, being included between 7.56 and 8.20 (the latter at a higher altitude, i.e., 495 m). It is worth noting that, for the foothills of the northern Italian Apennines, the regression slope slightly increases with altitude. Moreover, [45] have evidenced that such slopes become greater also if anomalous summer δ 18 O-δ 2 H pairs are excluded.
Except for Bologna, all the intercepts (between 8.6 and 11.6) are close to that of GMWL (10.0), indicating a dominant control on precipitations occurring in this area by air masses from the Atlantic Ocean. As already indicated by [31] lower values of intercept (4.2) detected in Bologna (we recall that this dataset is older than the others and regarded the period 1996-2000) was related to the proportion of air masses that had originated over different areas and in particular Central Europe (see, for example, the recent Meteoric Water Lines obtained by [49] and [50] for the Balkan zone: δ 2 H(% ) = 7.4 δ 18 O + 2.7 and δ 2 H(% ) = 7.74 δ 18 O + 5.6, respectively). These main areas of origin of the vapor masses condensating over the foothills of the northern Italian Apennines have been further confirmed employing the back trajectories calculation with different atmospheric tracers by [51].
It is worthwhile stating that the overall Italian maps of mean annual values of δ 18 O [44,52] agree in reporting similar isotopic contents for the whole study area.

Precipitation, River Discharge, and Groundwater Levels
Precipitation and river discharge data from the three rain and five river gauges in which sampled water for isotopic analyses come from [32]. They consist of mean monthly precipitation and discharge data from the considered time period (January 2004-December 2006 in case of rainfall and January 2005-December 2007 for discharge). It is worth noting that, upward of the apices, all the five catchments are characterized by a low degree of urbanization and no presence of remarkable artificial reservoirs and or diversions. Moreover, data of groundwater levels are made available by [32] and involve a reduced number of wells. They consist of four-monthly data that cover different periods. For the sake of convenience, in this paper we provide examples (as Supplementary Materials) of groundwater level changes in two wells from the apical part of Taro and Enza alluvial fans.

Methodology
The study consists of two steps involving univariate and bivariate statistical analyses, respectively. In both cases, statistical analysis are carried out by taking into account isotopic datasets of precipitations (rain gauge gauges numbered from 1 to 3: 1: Parma; 2: Lodesana; 3: Carpineti), surface water (letters from a to e: a: Trebbia; b: Taro; c: Enza; d: Secchia; e: Reno) and groundwater (isotopic values from wells belonging to each alluvial fan are aggregated and groundwaters are numbered with a capital letter from A to E; A: Trebbia, wells from I to IV; B: Taro, wells from V to IX; C: Enza, wells from X to XII; D: Secchia, wells from XIII to XVII; E: Reno, wells from XVIII to XX).
Firstly, such isotopic datasets are compared by utilizing univariate analysis involving box plots and some statistical estimates (means and quartiles, together with weighted and unweighted means; see Section 4.1). This step allows us to highlight the possible presence of discrepancies between the isotopic values of rainfall, surface water, and groundwater (all the isotopic data from wells belonging to the same alluvial fan are aggregated), which may be due to pre-infiltrative fractionation processes.
Successively, bivariate regressions among δ 18 O-δ 2 H pairs are exploited to verify which component (rainfall or river discharge) actively participates in recharging the groundwater hosted in the apical parts of alluvial fans (Section 4.2). Some regression characteristics (such as slopes of regression lines) may change if an evapotranspiration or evaporative process has occurred before infiltration and can be supposed to act more intensively in river water.
It is worth noting, in both steps, the role of the precipitations and river discharge amounts during the years taken into account. This is carried out by preliminary application of specific weighting procedures to isotopic datasets.

Univariate Analysis and Comparison between Isotopic Compositions in Precipitation, Surface Water and Groundwater
δ 18 O, δ 2 H, and d E of precipitations, surface water, and groundwater are first of all compared by utilizing box-plots. The latter has been widely used in hydrology and hydrogeology for identifying discrepancies among datasets by reporting isotopic values as vertical columns of univariate data [53,54]. In detail, each vertical column, which is centered on its median value, graphically depicts quartiles of its distribution (1st quartile 1Q and 3rd quartile 3Q, while the 2nd quartile is called median; see [55] for further statistical details). The spacings between the different parts of the box (i.e., 1Q and 3Q quartiles) indicate the degree of dispersion (spread) and skewness in the data together with outliers.
Successively, averages are computed for isotopic time series of precipitation and river discharge. In detail, the statistical analyses of isotope ratios (δ) and d E are carried out by both (4) arithmetic mean M A and (5) weighted mean M w as follows: where W i and δ i are the amount of precipitation (or river discharge) and the measured monthly isotopic compositions in the i-th month of the series, respectively. N is the number of observations.

Bivariate Analysis and Comparison among Water Lines
Water lines from precipitations (MWLs), rivers (RWLs), and groundwater (GWLs) are compared to verify possible discrepancies among the different datasets. These discrepancies are usually promoted by post-condensation and temperature-driven processes such as evaporation and evapotranspiration, which lead to a fractionation between the different isotopologues of water [56,57]. In fact, lighter water molecules ( 1 H 2 16 O) vaporize faster than heavier ones ( 2 H 2 18 O) inducing an enrichment of the latter into the liquid residual. As a water parcel evaporates, its isotopic composition usually evolve with a δ 18 O-δ 2 H regression line whose slope is lower than those of MWLs. Thus, a change in slopes from RWLs and GWLs concerning those of MWLs may be used as a groundwater recharge marker. In all cases (MWLs, RWLs, GWLs), the ordinary least squares regression (OLSR) approach is used (see [55] and [58] for more-in-depth statistical details). Here, we stress that the OLSR approach minimizes the sum of the squared vertical distances between the y data values and the corresponding y values on the fitted line (the predictions). Thus, the OLSR design assumes that there is no variation in the independent variable (x). It is noteworthy to say that this approach gives equal weighting to all the data points regardless of the amount of precipitation or discharge they represent.
To take into account the precipitation and discharge patterns, we have obtained the precipitation-weighted (WMWLs) and discharge-weighted (WRWLs) lines by applying the OLSR approach to the δ 18 O and δ 2 H time series in precipitation and rivers that are previously weighted with precipitation amount and discharge values, respectively. More in detail, data prior to weighting procedure of the former δ 18 O and δ 2 H time series in precipitation follows the procedure reported in [59], where weighted δ 18 O and δ 2 H time series (i.e., wδ 18 O(i)) is obtained by weighting the δ 18 Oi of monthly precipitation with its volume (P, in mm): where N is the number of months for which isotopes have been collected while δ 18 O G is the long-term mean input isotope ratio, which is calculated as: As in the case of the weighted δ 18 O time series from precipitation, also the weighted δ 18 O time series from river discharge can be obtained by taking into account the river discharge pattern. The wδ 18 Oout(i) time series for each river was reconstructed by using the Equations mentioned above (6) and (7), in which Q i (discharge at the i-month) is used instead of P i .
Both in the case of unweighted and weighted regression lines, the predictive performances of the OLSR models are tested utilizing Pearson's correlation coefficient coupled with the t-distribution test, in which the selected threshold p has been set at 0.05.

Univariate Analysis and Comparison between Isotopic Compositions in Precipitation, Surface Water and Groundwater
Box plots involving δ 18 O values from rain gauges (1, 2, 3), rivers (a, b, c, d, e) and wells from the apical part of the corresponding alluvial fans (groundwater: A, B, C, D, E) are shown in Figure 3 (δ 2 H not reported as being characterized by the same pattern of the δ 18 O). Table 4 summarizes median and quartiles 1Q and 3Q of each distribution. By considering the rain gauges, it can be noted that their δ 18 O distributions are broad and median values become lower with altitude, passing from −8.62% (1) to −9.99% (3). Table 2. Unweighted (MWLs) and Weighted Meteoric Water lines (WMWLs) for rain gauges considered in this study obtained from unweighted and weighted isotopic composition (in bold), respectively. The number of samples used for interpolations (ordinary least squares regression-OLSR) are also reported along with statistical parameters (R 2 and p-value).   From a general point of view, dispersions of δ 18 O values from rivers and wells are strongly reduced compared to those of precipitations ( Figure 3). This reduction is more marked in the case of groundwater. Moreover, δ 18 O distributions of rivers and wells are always included within those of precipitations, but some discrepancies can be evidenced. Except for the rain gauge (1), median values of surface water ranged between −7.56‰ (b) and −8.88‰ (d) lying above the corresponding values in precipitation. Moreover, Figure 3 highlights that δ 18 O distributions from surface water and groundwater are almost staggered (i.e., 1Q and 3Q from surface water do not include the corresponding groundwater quartiles).

Meteoric Water Lines (MWLs
Passing to the d E , all distributions are skewed while a marked increase in box marks and whiskers in surface water and groundwater can be noticed in comparison with δ 18 O (Figure 3). Letters from a to e and numbers from 1 to 3 recall surface water and precipitation, respectively. Capital letters from A to E refer to groundwater from the apical parts of alluvial fans (in form of aggregated isotopic values from wells). The whiskers represents the 10th and 90th percentiles, the box limits indicate 1st and 3rd quartiles (1Q and 3Q) and the line within the box mark is the median. For further details, readers are referred to Table 2 (precipitation) and Table 3 (surface water and groundwater). Table 4. Results of the statistical analyses (median, 1st and 3rd quartiles, carried out from δ 18 O and dE time series of precipitations (rain gauges; numbers from 1 to 3), surface water (rivers; letters from a to e) and groundwater from the apical pars of alluvial fans (in form of aggregated isotopic values from wells; capital letters from A to E). Median, 1st and 3rd quartiles are those from box plots (Figure 3). Means (MA and weighted means (Mw) are also reported. For further details, readers are referred to From a general point of view, dispersions of δ 18 O values from rivers and wells are strongly reduced compared to those of precipitations ( Figure 3). This reduction is more marked in the case of groundwater. Moreover, δ 18 O distributions of rivers and wells are always included within those of precipitations, but some discrepancies can be evidenced. Except for the rain gauge (1), median values of surface water ranged between −7.56% (b) and −8.88% (d) lying above the corresponding values in precipitation. Moreover, Figure 3 highlights that δ 18 O distributions from surface water and groundwater are almost staggered (i.e., 1Q and 3Q from surface water do not include the corresponding groundwater quartiles).
Passing to the d E , all distributions are skewed while a marked increase in box marks and whiskers in surface water and groundwater can be noticed in comparison with δ 18 O (Figure 3). Median values of d E from precipitations range between 11.8% (1) and 9.0% (2). In the event that surface water is taken into account, d E varies from 11.4% (c) to 12.3% (a) while groundwaters show slightly lower values, which range from 10.3% (A) to 11.7% (D). Table 2 also summarizes the mean (M A ) and weighted (Mw) values of δ 18 O and d E . δ 18 O-M A from precipitation decreases with altitude by passing from −9.04% in (1) to −9.67% in (3). The corresponding weighted averages are slightly lower than the aforementioned unweighted means being between −10.26% (1) to −11.12% (3). The same behaviour can be noticed in case of surface water, whose δ 18 O-M A vary from −7.70% (a) to −9.02% (d) while δ 18 O-Mw is between −7.75% (a) to −9.36% (d). Groundwater shows δ 18 O-M A in the range of −8.26% (B) and −8.55% (E). d E -Mw from precipitations are higher than the corresponding d E -M A, with values that vary from 11.2% (c) to 13.0% (a) and from 9.8% (c) to 11.0% (a), respectively. Weighted averages lead to an increase of d E also in case of surface water, with d E -M A included between 10.7% (d) and 12.3% (b) and d E -Mw between 11.1% (d) and 13.0% (c). Groundwaters are characterized by d E -M A values between 10.6% (A) and 12.1% (E).
The above-mentioned difference between M A and M W values is because precipitation and river discharge amounts, as well as their isotopic composition, are not uniformly distributed during the year (Figure 4 for δ 18 O and Figure 5 for d E ).  Concerning d E in surface water, the signal is more disturbed than those seen in the water from rain gauges. At first sight, four samples belonging to (d) and collected from January 2017 to April 2017 are anomalously low in such parameters (values up to 2.3% in April 2007). Higher values are usually found in winter seasons while lower values in the summer ones.
Both δ 18 O and d E varied during the considered period in groundwater. It must be specified that, for the sake of convenience, Figures 4 and 5 report one selected well for each apical (here, we recall that they have to be considered representative of each groundwater group A, B, C, D and E).

Bivariate Analysis and Comparison among Water Lines
In all cases (precipitations, surface water, and groundwater), δ 18 O and δ 2 H are correlated. The δ 18 O-δ 2 H relationships are summarized in Table 2 (Meteoric Water Lines MWLs from rain gauges), Table 3 (River Water Lines RWLs from rivers) and Table 5 (Ground Water Lines GWLs), in which slopes, intercepts, and coefficients of determinations (R 2 ) are reported separately. Tables 2 and 3 also include the corresponding weighted relationships for precipitations (WMWLs) and surface water (WRWLs), while supplementary materials ( Figure S4) provide plots for visual inspection. Table 5. Groundwater lines (GWLs) for the alluvial fans (apices) considered in this study. The number of water wells and samples used for interpolations (ordinary least squares regression-OLSR) are also reported along with statistical parameters (R 2 and p-value). Slopes obtained interpolating δ 18 O-δ 2 H pairs from rain gauges (MWLs) range between 7.5 (1) and 7.8 (2,3). If weighted δ 18 O-δ 2 H pairs are considered, slopes of WMWLs of rain gauge (2) and (3) increase to 8.4. Intercepts are always positive being between 6.2 (1) and 7.8 (2) in case of MWLs and 3.3 (1) and 13.1 (2,3) if WMLs are taken into account.

Groundwater Lines (GWLs)
RWLs slopes from surface water are slightly lower and are included between 4.3 (d) to 6.5 (b). As already pointed out by [31], (d) is characterized by R 2 value (equal to 0.33), indicating weaker correlations. The latter is due to 4 monthly samples from January 2017 to April 2017 lying far below the alignment of the others (see also Figure S4 in Supplementary Materials) as a result of sublimation processes acting on the snow cover in the uppermost part of its catchment [31]. In case these four anomalous samples are not taken into account, the slope of regression is still lowered (5.1), while statistical parameters improve remarkably (R 2 = 0.85 and p-value now lower than 0.0001).
In case weighted pairs of δ 18 O-δ 2 H are taken into account, slopes from WRWLs are between 1.9 (d) and 6.6 (c). Once again, (d) is characterized by the lowest R 2 (=0.15) with the consequent exceeding of the p-value threshold (that has been set equal to 0.05, as anticipated in Section 4.2). If the four samples from January 2017 to April 2017 are removed from the time series, the slope becomes equal to 3.3 (the most reduced slope among the five considered linear regressions from rivers) while p-value is now significant (lower than 0.0001 and R 2 = 0.81).
In almost all cases, intercepts are negative, with the lowest values detected in (4) for both unweighted (−22.8) and weighted (−46.0) case. As already evidenced by [31], such negativization is because these waters have undergone evaporation.
Passing to the groundwater, the δ 18 O-δ 2 H alignment show slopes falling between those of precipitations and rivers water, being included between 6.2 (D) and 8.0 (A). It should be pointed out that three groundwaters show intercepts close to that of GMWL and LMWLs from the area (A, B, E; values between 9.3 and 11.1). On the contrary, (C) and (D) are characterized by negative intercepts.

Discussion
The univariate statistical analyses allow us evidencing some discrepancies among the different datasets (precipitations, surface water, and groundwater). Firstly, box-plots make evident that the dispersion of δ 18 O in rivers and groundwater is reduced (Figure 3). This is due to the prevailing homogenization effect occurring within catchments and aquifers, which has been already highlighted in several areas worldwide [54,56,60,61].
Except for (c, C), box-plots involving δ 18 O from surface water and groundwater indicate that the majority of samples (i.e., those within the whiskers; Figure 3) are staggered with each other. This fact suggests that the contribution of surface water to groundwater recharge in these alluvial fan apices is somehow reduced.
Secondly, weighting procedures taking into account the precipitation and discharge patterns highlight that δ 18 O-M W are lower for both precipitations and surface water. Both in weighted (Mw) or unweighted (M A ) cases the means are taken into account, the isotopic signatures from surface water are always higher than those of precipitations. The different δ-signal in rivers with respect to that of precipitations is due to both the different isotopic content of air masses in the highest portions of the catchments (only them, Central Europe or Tyrrhenian; [31,62]) and temperature-driven fractionation phenomena occurring at the catchment scale and leading the heavier water molecules ( 2 H 2 18 O) to be enriched in the surface water [31]. As described in more detail below, such processes take place when drops reach the topsoil and are mainly related to partial evaporation from the soil moisture and evapotranspiration. The above-mentioned discrepancies among δ-values in precipitation and surface water are successively exploited in a bivariate analysis to verify whether a component of recharge for groundwater hosted within the apices is prevalent or not. The Local Meteoric Water lines (LMWLs) obtained by processing the pairs of δ 18 O and δ 2 H from samples collected from rain gauges were similar, with slopes ranging from 7.6 (1) to 7.8 (2,3). These values are close to that (equal to 8.0) characterizing both the World Meteoric Water Lines [47] and the Mediterranean Meteoric Water Line [48] and, more specifically, the foothills of the northern Italian Apennines (7.33: [45]; 7.73: [42]; 7.56: [46]). If we consider the weighted precipitation values, slopes forming WMWLs slightly increase up to 8.4 in the rain gauges (2,3) located at higher altitudes (i.e., the same altitude of the apical parts of the investigated alluvial fans). As already highlighted by [45] for the study area, such an increase in slope values is related to rainfall occurring between autumn and spring, which took place at lower temperatures. In all cases, slopes from rivers (RWLs) are much lower than those from precipitations (LMWLs; slopes included between 4.3 and 6.5), even taking into account the discharge weighted values (WRWLs; slope values between 1.9 and 6.6). As anticipated, it is widely recognized that the main causes inducing a lowering in slope from surface water compared to that of meteoric water have to be searched for in evaporative and evapotranspiration processes [56,57]. The latter play an evident role in temperate, continental and arid areas worldwide (see for instance [57,63,64]). Following [31], such a lowering in slopes from δ 18 O-δ 2 H relationships in rivers from the northern Italian Apennines is due to the base-flow fed into the area by a large number of low-yield springs whose water have been affected by pre-infiltrative evaporation/evapotranspiration. The latter is promoted by the spread of clay-rich material, in which low permeability values allow moisture to be partially evaporated before infiltration. It must be added that the same authors have excluded instream evaporative processes as no relationships between δ 18 O values with flow-lengths (and catchment areas) of rivers were found. Moreover, and for the northern Italian Apennines, [65] have recently identified marked diel fluctuations in stream water levels during the evapotranspiration period (from April to September) suggesting a dominant control of tree transpiration on streamflow, which in turn promotes the lowering of the slopes mentioned above in rivers (see, also the benchmark paper by [66]). Concerning groundwater, aggregated wells from some alluvial fans (A, B, E) show δ 18 O-δ 2 H relationships with slopes within the range of those from rain gauge water (from 7.8 to 8.0). The similarity between slope values leads us to confirm that almost all the recharge in these apices is linked to rainfall in the area. On the contrary, slopes from δ 18 O-δ 2 H pairs collected in groundwater from alluvial fans (C, D) are slightly lower (6.4 and 6.2, respectively), suggesting an important quota of recharge also from streambed dispersion. This is in agreement with intercept values that, with reference to the groundwater (A, B, E), are positive and similar to those of precipitation. On the contrary, isotopic values from aggregated wells in (C) and (D) are negative (−2.8 and −4.5, respectively), indicating that these waters have undergone evaporation as the case of river waters (c) and (d). If δ 18 O-d E pairs from groundwater are reported in Figure 6 together with the corresponding mean isotopic values from rain gauge (rain gauge 3; we recall that this is the most representative rain gauge for the apical parts of the alluvial fans) and rivers, further confirmation can be obtained. It can be evidenced that, in the case of (A, B, E), the majority of the dots are related to precipitations (i.e., they fall within the area identified by the standard deviations (1σ) of precipitations and not that of river water). On the contrary, δ 18 O-d E pairs from apices of alluvial fans (C, D) are more dispersed with several dots that may be related either to precipitations and river waters.
Although such a scattered isotopic dataset allowed us to preliminarily discriminate against the recharge mechanisms, their quantification is still a challenge that requires further effort. The latter should be based on a more comprehensive approach, which involves a set of different environmental tracers and in-continuous monitoring of river discharge, precipitations, and groundwater levels. In particular, we believe that in-continuous monitoring of all the components is necessary and should involve both quantities (i.e., precipitation amount, river discharges, and groundwater levels) and some selected physico-chemical parameters of water (such as TDS and temperature). Such approach allows us to analyse each recharge event and, in particular, to differentiate those related to the recharge by stream dispersion alone (river flood generated by rainfall in the upper part of the catchments) from those due to zenithal recharge by precipitation alone (occurring only on the apical part of the alluvial fan and with river fed by base flow). Furthermore, the implementation of an in-continuous monitoring network would be preparatory for the application of a suite of environmental tracers (stable isotopes of water and tritium) to be obtained at the event scale (that of flood and or rainfall with, at least, hourly sampling activities). This allows end-member mixing analysis to be carried out at the event scale and aquifer recharge quotas from each component to be better constrained.

Conclusions
This study has qualitatively demonstrated that three of five alluvial fans from the northern Italian Apennines (those related to the Trebbia, Taro, and Reno rivers) are likely to be mainly recharged by precipitations infiltrating in their apical areas. Wells located in the apical part of the alluvial fans from the other two rivers (Enza and Secchia) have shown a more complex behaviour suggesting either precipitations and river dispersions are currently participating in the aquifer recharge. These results are obtained by comparing the isotopic content (δ 18 O, δ 2 H, d E ) of precipitations, river water, and groundwater, with univariate (box-plots) and bivariate (slopes from δ 18 O-δ 2 H regressions and δ 18 O-d E plots) statistic approaches. The usefulness of these isotopic tools is here somehow enhanced by the isotopic modification (mainly related to sublimation and evaporative/evapotranspiration processes) occurring in water from rivers when traveling along their catchments before the apical part of the alluvial fans. That allows exploitation of the already available time series of scattered water isotopes (monthly samples of precipitations and river waters, four-monthly samples in the case of groundwater) and a first characterization of the mechanisms of recharge of the aquifers hosted within these alluvial fans. We think that the results achieved by this study are important for two main reasons. Firstly, they represent the basis for the subsequent quantification of the recharge of aquifers facing the northern Apennines of Italy. This information is needed by stakeholders in charge of groundwater management for further adaptation plans, such as managed recharge and the development of agricultural best practices capable of inducing an increase in groundwater levels and water quality. Secondly, they demonstrate that the fractionation processes affecting water isotopes from rivers can be exploited to verify whether a component of recharge for groundwater is prevalent or not. Since water of a large number of rivers from temperate, continental and arid areas worldwide have proven to be affected by fractionation processes, we think the use of the presented approach (based on scattered water samples for which stable isotopes of water are available) can be useful for other hydrogeological studies aimed at determining the interaction between surface water and groundwater.
Author Contributions: Conceptualization, F.C. and G.M.; methodology, F.C.; formal analysis, F.C.; data curation, F.C. and A.D.; writing-original draft preparation, F.C.; writing-review and editing, F.C., G.M. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.