Inferring Hydrological Information at the Regional Scale by Means of δ¹⁸O–δ²H Relationships: Insights from the Northern Italian Apennines

Abstract

We compared five regression approaches, namely, ordinary least squares, major axis, reduced major axis, robust, and Prais–Winsten to estimate δ¹⁸O–δ²H relationships in four water types (precipitation, surface water, groundwater collected in wells from lowlands, and groundwater from low-yield springs) from the northern Italian Apennines. Differences in terms of slopes and intercepts of the different regressions were quantified and investigated by means of univariate, bivariate, and multivariate statistical analyses. We found that magnitudes of such differences were significant for water types surface water and groundwater (both in the case of wells and springs), and were related to robustness of regressions (i.e., standard deviations of the estimates and sensitiveness to outliers). With reference to surface water, we found the young water fraction was significant in inducing changes of slopes and intercepts, leading us to suppose a certain role of kinetic fractionation processes as well (i.e., modification of former water isotopes from both snow cover in the upper part of the catchments and precipitation linked to pre-infiltrative evaporation and evapotranspiration processes). As final remarks, due to the usefulness of δ¹⁸O–δ²H relationships in hydrological and hydrogeological studies, we provide some recommendations that should be followed when assessing the abovementioned water types from the northern Italian Apennines.

1. Introduction

Oxygen (¹⁸O and ¹⁶O) and hydrogen isotopes (²H and ¹H) of water are commonly used in surface and subsurface hydrology [1,2,3]. They are considered environmental tracers in the form of δ¹⁸O and δ²H, where $\mathsf{\delta }(‰)=\left[\left(\frac{R_{sample}}{R_{standard}}-1\right)\times 1000\right]$ and R is the corresponding isotopic ratio (¹⁸O/¹⁶O or ²H/¹H) in the water sample or in the standard (usually V-SMOW, i.e., the Vienna Standard Mean Oceanic Water). If a multitude of water samples is collected from the same source (a rain gauge, a river, a spring, or a well), the corresponding δ¹⁸O–δ²H pairs in a Cartesian graph will be aligned along a regression line in the form of y = mx + q, where y is δ²H, x is δ¹⁸O, m is the slope, and q the intercept. This fact was first noted by [4] when reporting several δ¹⁸O and δ²H values from precipitation waters worldwide, which allowed definition of the so-called “global meteorological water line” (GMWL) equal to δ²H(‰) = 8.0 · δ¹⁸O + 10.0. The authors of [5] have substantially confirmed this slope and intercept by processing more recent data. Although this relationship is valid everywhere, more accurate regression lines (with different slope and intercept) can be obtained by selecting δ¹⁸O and δ²H values from restricted areas. This is due to the specific isotopic fractionation processes (i.e., vapour pressure and temperature conditions) controlling precipitation over each area. Moreover, since further post-condensation and temperature-driven processes such as evaporation and evapotranspiration could act prior to infiltration and/or during runoff, δ¹⁸O–δ²H regressions from rivers (river water lines, RWLs) and groundwater (groundwater lines, GWLs) may also differ from that of the precipitation occurring in their recharge areas (meteoric water lines, MWLs). In fact, evaporation and evapotranspiration lead to a fractionation between the different isotopologues of water, with lighter water molecules (¹H₂¹⁶O) vaporising faster than heavier ones (²H₂¹⁸O) and inducing an enrichment of the latter into the liquid residual. In this case, as a water parcel evaporates, its isotopic composition usually evolves with a δ¹⁸O–δ²H regression line whose slope is lower than those of MWLs.

It is evident that such changes in slopes from RWLs and GWLs concerning those of MWLs can be used to infer information on the hydrological processes occurring at the slope and the catchment scales. As an example, and without claiming to be exhaustive, the following studies highlighted that a change in slope from RWLs and GWLs concerning those of MWL can be used to:

• Display the role of the riparian zone in feeding base flow in low-relief and forested catchments [6];
• Highlight pre-infiltrative evaporation/evapotranspiration that acted by modifying the waters before their infiltration towards the aquifer and, subsequently, to the base flow of rivers [7];
• Demonstrate that groundwater may be fossil (related to other climate recharge condition) or actually recharged by losses from the streambed or even a mixing among these two components [8];
• Reveal instream or lake evaporation in a nonarid environment [9,10];
• Elucidate which component (precipitation or losses from the streambeds) is exclusive or prevalent in recharging alluvial aquifers [11].

Starting with the pioneering works by [4,12] regarding meteoric water and up to this day, δ¹⁸O–δ²H regression lines are usually carried out by means of the ordinary least squares (OLS) method, i.e., an approach that minimises 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 OLS design assumes that there is no variation in the independent variable (x) and is considered as the simplest method among the several available linear regression models. By focusing on the aforementioned isotopes of water, we should also take care of the variations associated with variable (x) as the same or different isotopic fractionation processes, which may have developed even at different rates and may have affected both δ-values.

For this reason, [13] proposed a more complex linear regression approach for obtaining MWLs lines that consider associated errors with both dependent (y) and independent variables (x), such as the reduced major axis (RMA) and the major axis (MA). In the end, it is found that MA approaches usually led to the smallest discrepancies between the estimated and predicted values (a measure of goodness of fit usually described with the well-known coefficient of determination R², i.e., smallest discrepancies are identified by higher values of R²) and larger slopes in MWLs calculation than those obtained with RMA and OLS [14]. The recent attempt made by [15] on several water types (river water, groundwater, soil water) confirmed that RWLs and GWLs were also characterised by larger slopes in the case of MA regression while the lowest values were noted when using OLS. The same authors found that in all cases (MWLs, RWLs, GWLs), the higher the R² between δ¹⁸O and δ²H values, the smaller were the differences in the slopes obtained by MA, RMA, and OLS. This was due to the different sensitiveness of the linear regression approaches to outliers and large measurement errors rather than temperature-driven post-condensation processes.

This study aims to verify whether such discrepancies in slopes and intercepts from different regression methods are present (thus significant) or not in four water types (precipitation, surface water, groundwater collected in wells from the lowlands, groundwater from low-yield springs) from the northern Italian Apennines. For this reason, we exploited datasets already published in the literature, e.g., [7,11,16,17,18], and we carried out visual inspection (heat maps) and statistical comparison of results from the three aforementioned approaches (OLS, RMA, MA) already tested in [15]. With reference to OLS approach, we further verified whether preliminary weighting of the isotopic data to the corresponding values of discharge or precipitation may have induced changes to our results or not. In addition, we tested two methods, namely, Prais–Winsten (PW) and robust (R), to investigate possible influence on the final δ¹⁸O–δ²H alignments of nonstationary processes (here, we recall that the OLS, RMA, and MA approaches are based on the assumption that data are not serially correlated, thus a δ¹⁸O–δ²H pair from a determined time period is not correlated with the earlier one, while properties as mean, variance, and autocorrelation are constant over time, i.e., stationary, while recent papers in the literature highlighted the possible nonstationary behaviour of such series of isotopic data [19]) or even outliers (i.e., anomalous isotopic values) within the series of isotopes.

Furthermore, we provided a possible explanation for the geographic and climatic factors (i.e., catchment descriptors) influencing the several regressions and finally we made some considerations concerning their applicability in the context of mountainous areas such as the northern Italian Apennines.

2. Overview of the Climatic, Geomorphological, and Hydrogeological Features of the Study Area

The study area is located in the northern Italian Apennines and belongs to the Emilia–Romagna Region (Figure 1). It includes nine catchments between the Trebbia River and the Savio River, with river gauges (in which the samples were collected) located close to the foothills of the mountain chain. The area has an overall extension of 6261 km². Maximum altitudes are in the southern sectors, where the main watershed divide lies (with main peaks showing elevations higher than 2000 m a.s.l., such as Mt. Cimone with its 2165 m a.s.l.) is the southern border of the Emilia–Romagna Region. Elevation decreases towards the NE direction to approximately 40 m a.s.l. at the Savio River gauge.

All the nine rivers originate from the main watershed divide and flow towards the NE. Six rivers (namely, Trebbia, Nure, Taro, Enza, Secchia, and Panaro) are tributaries of the Po River while the other three (Reno, Lamone, and Savio) enter the Adriatic Sea. Catchment areas are between 193 km² (Lamone) and 1300 km² (Secchia), while flow lengths range from 28 km (Enza) to 85.2 km (Secchia). Mean annual discharges during the period 2006–2016 are included between 8.4 m³ s⁻¹ (Savio) and 30.4 m³ s⁻¹ (Secchia).

From a hydrogeological point of view, a report [20] grouped bedrock outcrops in the aforementioned catchments into six main classes (or hydrogeological complexes, namely: clay, marl, flysch, foreland flysch, ophiolite, and limestone). Those composed of poorly permeable or impermeable materials (clay, marl, and flysch hydrogeological complexes; see Figure 1) are the most represented in terms of areal coverage, leading to a runoff response of rivers that closely follows the rainfall distribution during the year (pluvial discharge regime). Rivers originating from the most elevated parts of the main watershed divide (Secchia, Panaro) are characterised by a nival–pluvial discharge regime as they are influenced by the melting of snow cover accumulated during the winter months in the upper parts of their catchments [21].

In the vicinity of the foothills (therefore close to the corresponding river gauges), several wells drilled in the alluvial fans of the Trebbia, Taro, Enza, and Secchia rivers continuously pump groundwater for both agricultural and drinking purposes. As previously highlighted by [11], by exploiting water stable isotopes, wells pumping water from confined aquifers in Trebbia and Taro alluvial fans are also likely to be recharged by zenithal precipitation infiltrating through gravels and sands that outcrop at the foothills of the northern Apennines (i.e., apical part of the alluvial fans). On the contrary, an important quota of recharge also occurs from streambed dispersion (focused on the apical part of their alluvial fans, see [22]) seems to affect wells located in the alluvial fans of the Enza and Secchia rivers. Two groups of low-yield springs from the Secchia River catchment (namely, Pietra di Bismantova and Montecagno) were also considered. These springs should be considered as representative of the common ones in the northern Italian Apennines, whose discharges are closely related to the rainfall pattern while outflows are strongly reduced (often in the order of 1 L·s^–1 or less) at the end of the summer periods (shallow groundwater flow paths; for more details see [20,23,24]).

From the climatic point of view, and as already reported in [25], the mean annual rainfall distribution during the period 1990–2015 exceeds 2200 mm/y near the main watershed divide and progressively decreases to about 900 mm/y in the foothills. The rainfall distribution during the year is characterised by a marked minimum in the summer season and two maxima during autumn (the main one) and spring. Close to the main watershed divide, the cumulative annual snow cover can reach 2–3 m at the end of the winter season [21]. Potential evapotranspiration ranges from about 500 up to 650 mm/y in the lowlands and is mainly active during the summer months [25].

3. Methodology

The methodology used in this study consists of five steps involving data inspection and statistical comparison on datasets, including rainfall (4 rain gauges, namely: Parma, Lodesana, Langhirano, Berceto), surface water (9 rivers, namely: Trebbia, Nure, Taro, Enza, Secchia, Panaro, Reno, Lamone, Savio) groundwater from wells (aggregated isotopic values from wells belonging to 4 alluvial fans, namely: Trebbia—4 wells, Taro—5 wells, Enza—3 wells, and Secchia—5 wells), and groundwater from springs (aggregated isotopic values from springs belonging to 2 areas from the Secchia catchment, namely: Pietra di Bismantova and Montecagno). Firstly, a check on the assumption of stationary behaviour of each series of stable isotopes was carried out by means of conventional statistical tests coupled with inspection of standardised residuals.

Secondly, slopes and intercepts from δ²H–δ¹⁸O alignments were obtained by using 5 different regression approaches. Moreover, we further considered 2 different regressions applied to δ-values from rain gauges and rivers that had been preliminary weighted to the corresponding quota of precipitation and discharge, respectively. Thirdly, slopes and intercepts were visually inspected by means of heat maps to identify discrepancies among the several regression methods. Fourthly, slopes and intercepts were compared through bivariate (correlation matrices) and multivariate analyses (hierarchic clustering, i.e., dendrograms) to identify linear correlations and similarities. Fifthly, and with reference to the only surface water, we made a comparison between differences in slopes and intercepts with some selected catchments to verify linear or nonlinear correlations among these variables.

For convenience (see Figure 1 for location of the sampling points and

Table 1. Main features of the sampling points from which isotopic data were derived. For the corresponding map locations, readers are referred to Figure 1.

Location	Type	Code	Number of Samples	Timing of Sampling	Length of Time	References
Parma	Precipitation	a	41	Monthly	3 January–6 December	[7,16]
Lodesana	Precipitation	b	18	Monthly	3 December–5 May	[7,16]
Langhirano	Precipitation	c	18	Monthly	3 December–5 May	[7,16]
Berceto	Precipitation	d	14	Monthly	4 September–5 October	[7,16]
Trebbia	Surface water	1	36	Monthly	5 January–7 December	[16]
Nure	Surface water	2	24	Monthly	6 January–7 December	[16]
Taro	Surface water	3	36	Monthly	5 January–7 December	[16]
Enza	Surface water	4	24	Monthly	6 January–7 December	[16]
Secchia	Surface water	5	24	Monthly	6 January–7 December	[16]
Panaro	Surface water	6	36	Monthly	5 January–7 December	[16]
Reno	Surface water	7	24	Monthly	6 January–7 December	[16]
Lamone	Surface water	8	24	Monthly	6 January–7 December	[16]
Savio	Surface water	9	36	Monthly	5 January–7 December	[16]
Trebbia	Groundwater from wells	A	66	Four-Monthly	4 January–7 December	[16]
Taro	Groundwater from wells	B	23	Four-Monthly	4 January–7 December	[16]
Enza	Groundwater from wells	C	23	Four-Monthly	4 January–7 December	[16]
Secchia	Groundwater from wells	D	33	Four-Monthly	4 January–7 December	[16]
Pietra di Bismantova	Groundwater from springs	α	32	Monthly Two-Monthly Three-Monthly	14 January–15 December	[17]
Montecagno	Groundwater from springs	β	21	Monthly Two-Monthly Three-Monthly	14 March–15 December	[18]

for further details), we report below rain gauges signified by letters (from a to d: “a”—Parma, “b”—Lodesana, “c”—Langhirano, “d”—Berceto); surface water locations as numbers (from 1 to 9: “1”—Trebbia, “2”—Nure, “3”—Taro, “4”—Enza, “5”—Secchia, “6”—Panaro, “7”—Reno, “8”—Lamone, “9”—Savio); groundwater from wells as capital letters (“A”—Trebbia, “B”—Taro, “C”—Enza, “D”—Secchia); and groundwater from springs as Greek letters (“α”—Pietra di Bismantova, “β”—Montecagno).

3.1. Isotopic Datasets

The dataset from 4 rain gauges and 9 rivers consists of monthly isotopic data while 17 water wells were characterised by grabbed four-monthly samples. All the data are derived from [7,16]. With reference to rain gauges, isotopic datasets lasted over the period from January 2003 to December 2006. Isotopic data from surface water and groundwater covered the period from January 2004 to December 2007. Further monthly, two-monthly, and three-monthly isotopic data from 2 groups of nearby low-yield springs were considered. These data were published in [17,18] and included the period January 2014 to December 2016.

The final dataset consists of 553 isotopic values, of which 91 are from precipitation, 264 from surface water, 145 from groundwater collected in wells, and 53 from groundwater collected in springs. Precipitation and river discharge that were used for further weighting procedures (see Section 3.3, “Linear Regression Types”) come from [26].

As reported in the previous works of [7,11,16], the isotopic analyses were carried out by using isotope ratio mass spectrometry (IRMS) while instrument precision (1σ) was on the order of ±0.05‰ for δ¹⁸O and ±0.7‰ for δ²H. With reference to groundwater from springs [17,18], the corresponding isotopic analyses were carried out by mixed technique involving IRMS for δ¹⁸O and cavity ring-down spectroscopy (CRDS) for δ²H. Instrument precision (1σ) was assessed as ±0.1‰ for δ¹⁸O and ±1.0‰ for δ²H.

3.2. Verifying the Stationary Behaviour of Isotopic Data Series

As anticipated in the introduction, three of the five linear regression approaches are considered in this study, namely, ordinary least squares (OLS), reduced major axis (RMA) and major axis (MA), which are based on the assumption that a series of stable isotopes are stationary [27,28]. This means that statistical properties as mean and variance remain constant along each δ¹⁸O–δ²H alignment, i.e., modelling errors are normally distributed and homoscedastic. Moreover, the closest pairs of δ¹⁸O–δ²H must not be affected by autocorrelation phenomena as the modelling errors (i.e., residuals) from the regressions must be independent [27].

The presence of outliers or heteroscedasticity (i.e., modelling errors have not the same variance over the alignment) or autocorrelation lead the assumption of stationarity to be violated, thus slopes and intercepts from the abovementioned regression may not be meaningful. In this work, we exploited conventional statistical tests for verifying multivariate normality (i.e., presence of outliers inducing non-normality; Doornik–Hansen test [28,29]) heteroscedasticity (Breusch–Pagan test [30]), and autocorrelation (Durbin–Watson [31]).

All of the abovementioned tests are based on a comparison of the corresponding statistics’ p-value results with a threshold value (level of significance α set at 0.01) to decide whether the null hypotheses have to be rejected (p < 0.01) or not (p > 0.01). The following null hypotheses were selected: multivariate normality (Doornik–Hansen), homoscedasticity (Breusch–Pagan), no autocorrelation (at a lag of 1) in the residuals (Durbin–Watson). With reference to the Durbin–Watson test, it must be specified that values are included between 0 and 4. Values close to 0 indicate almost total positive autocorrelation while results in the proximity of 4 indicate a total negative autocorrelation. Values between 1.5 and 2.5 show there is no autocorrelation in the data.

As suggested by [27], in the case of rejection of the null hypothesis of multivariate normality or homoscedasticity, we further verified the standardised residuals for identifying outliers (i.e., those values for which standardised residuals fall outside the interval from −4 to 4). Moreover and always following [27], if residuals were found to be serially correlated, we made autocorrelation functions of the standardised residuals to further confirm the lag length of autocorrelation.

3.3. Linear Regression Types

In this work, 5 types of linear regressions were tested, namely, ordinary least squares (OLS), reduced major axis (RMA), major axis (MA), robust (R), and Prais–Winsten (PW). To avoid excessive mathematical details, we provide a cursory examination of the methods and the reader is referred to more specific literature [14,15,27,32]. For convenience, we report all the corresponding equations by replacing δ¹⁸O with x and δ²H with y, while n is the number of samples and r is the Pearson correlation coefficient.

The OLS regression assumes that the x values are fixed (i.e., it is commonly used when x values have very few associated errors) and finds the line that minimises the squared errors in the y values. The slope of the linear regression (slope_OLS) is calculated as follows:

$$slop{e}_{OLS}=r\times \frac{s{d}_y}{s{d}_x}$$

(1)

where sd represents standard deviations calculated for x variables (sd_x) and y variables (sd_y). Unlike OLS, RMA and MA try to minimise both the x and the y errors [33]. In the case of RMA, the corresponding slope_RMA can be obtained with:

$$slop{e}_{RMA}=sign\left[r\right]\times \frac{s{d}_y}{s{d}_x}$$

(2)

where sign[r] is the algebraic sign of the Pearson coefficient. The slope_MA is calculated by:

$$slop{e}_{MA}=-A+\frac{\sqrt{r^2+A^2}}{r}$$

(3)

where A can be obtained as:

$$A=0.5\times \left(\frac{s{d}_x}{s{d}_y}-\frac{s{d}_y}{s{d}_x}\right)$$

(4)

As anticipated in introduction, the PW regression [34] has never been used for developing δ¹⁸O–δ²H alignments, as series of isotopic data have always been considered to the present as time-invariant (i.e., stationary). Recently, [19,35] highlighted that multiannual series of such isotopic data may be affected by nonstationary processes (such as, for example, trends in the means or presence of far-off values). In this case (nonstationary multiannual series of δ¹⁸O–δ²H pairs), the use of common regression methods such as OLS, RMA, and MA could induce residuals to be larger and characterised by stronger serial correlations. PW is commonly used for data with serially correlated residuals of the estimates [36]. As a matter of fact, this approach takes into account AR1 (i.e., autoregression of the first order) serial correlation of the errors in a linear regression model. The procedure recursively estimates the coefficients and the error autocorrelation of the model until sufficient convergence is reached. All the estimates are obtained by the abovementioned OLS.

As in the case of PW approach, the R method has also never been tested for δ¹⁸O–δ²H regressions. This method is an advanced Model I (in which x is always the independent variable) regression which is less sensitive to outliers than OLS estimates. Having less restrictive assumptions, R is recognised to provide much better regression coefficient estimates than OLS when outliers are present in the data. In particular, this approach has proven to be successful in the case of “almost” normally distributed errors but with some far-off values. This happens as outliers usually violate the assumption of normally distributed residual in OLS method. The algorithm is “least trimmed squares” reported by [37], in which the method consists of finding that subset of x–y pairs whose deletion from the entire dataset would lead to the regression having the smallest residual sum of squares. As in the case of PW approach, estimates of each subset are calculated owing to the OLS method. It must be added that, depending on the size and number of outliers, R regression conducts its own residual analysis and downweight or even these x–y pairs; this fact deserves an accurate inspection of the outliers made by the operator prior to any removal in order to decide whether these x–y pairs have to be considered or not.

For all 5 different regression approaches, the corresponding intercept is obtained with the following:

$$intercept=\frac{{{\displaystyle \sum }}_{i=1}^n{y}_i}{n}-slope\frac{{{\displaystyle \sum }}_{i=1}^n{x}_i}{N}$$

(5)

In all the abovementioned linear regressions, each observation has an equal influence of the orientation of the fitted line. As a matter of fact, it is well recognised that some isotopic data may be more important than others as related to a higher amount of water (for example, a flood in the case of a river or a high discharge event of a spring or high rainfall amount during a storm event). In this case, greater influence in the regression should be given to these isotopic data. In order to also take this effect into account, OLS were applied to isotopic datasets that had been previously weighted on the corresponding monthly amount of precipitation (rainwater) and discharge (freshwater and river water) by means of two different methods, i.e., the classical one that simply involves multiplying each y_i by the water amount (see [28] for further details; hereafter called W) and as reported in [38] (see [7] for the formulation; hereafter called B).

3.4. Comparison among Regressions

Initially, slopes and intercepts from all the regressions were compared by means of heat maps. The heat maps are matrices of fixed cell size showing the magnitude of difference among values with a selected binary colour ramp (in our case from red to green, respectively, indicating the lowest value and highest value within the isotopic dataset), in which the colour intensities provide visual cues to the reader about discrepancies between the data. The goal was to verify any presence of clusters to be further investigated by means of bivariate and hierarchical cluster analyses.

As a bivariate analysis, we carried out scatterplot matrices to determine if linear correlation between multiple variables (slopes and intercepts obtained by the different regression approaches) were present or not. Tests were carried out highlighting the level of significance, which was set as p < 0.01.

Furthermore, a hierarchical cluster analysis was performed to identify similarities among the series of slopes and intercepts from the whole datasets. Clustering was done according to the unweighted pair–group average (or centroid) method, in which each group consisted of slopes (or intercepts) from a determined regression approach. The method was based on a step-by-step procedure in which series of slopes (or intercepts) were grouped into branched clusters (dendrogram) based on their similarities to one another. As a result, the two most similar series of slopes (or intercepts) were selected and linked based on the smallest average distance among the values of all slopes (or intercepts). Progressively more dissimilar series were linked at greater distances; in the end, they all were joined to one single cluster. The cophenetic coefficient was used as a measure of similarity between each pair of clusters; more than 2 time series being analysed, the dendrogram was supported by a cophenetic distance matrix. Further details on this method can be found in [28].

3.5. A Focus on the Differences among Regression Approach from River Water: Comparison with Catchment Characteristics

As suggested by [15], we investigated whether some selected catchment characteristics (also called descriptors) could have affected differences among values of slopes and intercepts as obtained by the different approaches reported in Section 3.3. In order to make all slopes and intercepts comparable, we followed the approach proposed by [15] that consisted of prior computed differences in the slopes (as slope_OLS–slope_{RMA/MA/R/PW/W/B}) and intercepts (as intercept_OLS–intercept_{RMA/MA/R/PW/W/B}). Then, and following again the procedure reported in [15], we applied the Spearman ranking correlation matrix in which the abovementioned differences in slopes and intercepts were compared with 9 catchment characteristics. It should be added that this approach is a nonparametric measure of rank correlation that provides a statistical dependence between the rankings of two variables at a time. Unlike the Pearson coefficient, the Spearman ranking assesses how well the relationship between two variables can be described using a monotonic function, even if their relationship is not linear [27]. In particular, several authors have highlighted that many hydrological processing occurring at both slope and catchment scales are nonlinear (see for instance [38,39]) and such behaviour was in turn seen in some descriptors calculated by means of time series of stable isotopes of water [7,40,41,42].

In order to take into account the linearity among the variables, we also considered the linear correlation by providing the Pearson correlation coefficients. Here, we recall that Pearson and Spearman matrices reflect the magnitude of similarity among the parameters by means of r (the Pearson correlation coefficient described in Section 3.3) and r_s coefficients, respectively. Both correlation coefficients (r and r_s) describe the strength and direction between the two variables and return a closer value to 1 (or −1) when the two different datasets have a strong positive (or negative) relationship. The significance probability (p-value) for both the Spearman and Pearson correlation coefficients calculated in this study was set at 0.01, meaning that p-values lower than 0.01 represented statistically significant relationships. Readers are referred to [28] for further details on statistical formulations.

The 9 catchment characteristics (or descriptors, see

Table 2. The 9 catchment characteristics included in the analysis. For further details on the catchment characteristics from a single catchment, see Table S1 in Supplementary Materials.

Acronym	Variable	Units	Minimum	Mean	Maximum
A	Catchment area	km2	193	696	1303
Hmin	Altitude of stream gauge	m	43	171	421
Hmax	Maximum altitude	m	1158	1784	2165
Hmean	Mean altitude	m	526	754	944
P	Precipitation	mm	924	1090	1304
F	Flow length	km	20.9	55.5	85.2
q	Specific mean annual runoff	L s−1 km−2	2.2	15.0	36.3
q95	Specific runoff exceeded for 95% of the time	L s−1 km−2	0.0	1.0	1.7
Fyw	Young water fraction	%	9.3	13.7	22.9

and Supplementary Materials Table S1 for further details) were those already considered by [7], namely: catchment area (A); elevation (H); precipitation (P); flow length (F); specific mean annual runoff (q); specific river runoff exceeded for 95% of the observation period (q95; this is a well-known low flow index that is used worldwide for the regionalisation procedure and can be estimated even from a relatively short time series of daily runoffs [43]); and the young water fraction (F_yw proposed by [42]; this is considered the percentage proportion of catchment outflow younger than approximately 2–3 months and was estimated from the amplitudes of seasonal cycles of stable water isotopes in precipitation and stream flow that had been already calculated in [7]).

It should be noted that 4 descriptors (P, q, q95, F_yw) were obtained by processing daily precipitation (42 rain gauges homogeneously distributed over the study area for P), discharges (for both q and q95), and water isotopes time series (for F_yw) lasting over the same time period. Moreover, flow length (F) was derived from a 5 × 5 m gridded digital terrain model created by the digitalisation and linear interpolation of contour lines represented in the regional topography map at a scale of 1:5000.

4. Results

4.1. Stationary Behaviour of Isotopic Data Series

Table 3. Results from the three statistical tests aimed at verifying the compliance with the stationary assumption, namely: Doornik–Hansen (multivariate normality), Breusch–Pagan (homoscedasticity), and Durbin–Watson (autocorrelation). * Null hypotheses rejected as p < 0.01.

Location	Type	Code	Doornik–Hansen	Breusch–Pagan	Durbin–Watson
Parma	Precipitation	a	1.81	0.17	1.57
Lodesana	Precipitation	b	1.54	2.40	3.65 *
Langhirano	Precipitation	c	1.45	0.76	1.43
Berceto	Precipitation	d	3.83	0.46	0.84
Trebbia	Surface water	1	5.44	7.84 *	1.07 *
Nure	Surface water	2	1.50	0.36	1.41
Taro	Surface water	3	3.23	0.09	1.31
Enza	Surface water	4	1.94	0.90	0.90 *
Secchia	Surface water	5	10.30 *	0.19	1.21
Panaro	Surface water	6	8.98 *	0.86	0.69 *
Reno	Surface water	7	5.89	0.15	1.02 *
Lamone	Surface water	8	7.58	6.95 *	1.44
Savio	Surface water	9	41.42 *	0.00	2.30
Trebbia	Groundwater from wells	A	1.30	1.14	2.05
Taro	Groundwater from wells	B	3.67	0.18	2.28
Enza	Groundwater from wells	C	8.21	3.74	2.42
Secchia	Groundwater from wells	D	6.15	0.41	1.41
Pietra di Bismantova	Groundwater from springs	α	7.98	2.61	2.17
Montecagno	Groundwater from springs	β	7.76	0.37	1.21

summarises the results from the three statistical tests (Doornik–Hansen for multivariate normality, Breusch–Pagan for homoscedasticity, and Durbin–Watson for autocorrelation) used for assessing the compliance with the stationary assumption. Isotopic series from rivers were those mainly affected by problems of non-normal behaviour (rivers “5,6,9”) and autocorrelation at a lag of 1 (rivers “1,4,6,7”). The latter were positive (we recall here that values closer to 0 identify positive autocorrelation phenomena) and more intense in the case of rivers “4,6”. Moreover, the Breusch–Pagan test suggested residual homoscedasticity for river “8”. By considering the plots of standardised residuals (see Figure S1 in Supplementary Materials), the presence of outliers was further confirmed for rivers “5,6,9” (river “5”, 3 outliers; river “6”, 4 outliers; river “9”, 1 outlier) as well as the increase of variance of standardised residuals along estimates for river “1” (i.e., heteroscedasticity). Autocorrelation functions carried on standardised residuals (see Figure S2 in Supplementary Materials) allowed for demonstrating the presence of serial correlations, although with different lag lengths (river “1”, 2 lags; river “4”, 2 lags; river “6”, 3 lags; river “7”, 2 lags).

4.2. Slopes and Intercepts

The δ¹⁸O–δ²H relationships are summarised in in Supplementary Materials containing slopes (a; see Table S2), intercepts (b; see Table S3), standard deviation of the estimates (c; see Table S4), standard deviations of the estimates and coefficient of determinations (d: see Table S5) coefficient of determinations R². By viewing all the results reported in the form of heat maps (see Figure S3 in Supplementary Materials), the substantial invariance of slopes (from 6.9 to 13.1) and intercepts (from 7.2 to 8.4) from rain gauges located at lower altitudes (a, b, c), with high performance of the regression (R² always close to 0.99) was noticed. These results are in agreement with the GMWL (we recall that this line is characterised by slope and intercept equal to 8.0 and 10.0, respectively; see Figure S4 in Supplementary Materials) with no evidence of outliers. When the two weighting approaches were taken into account, no changes among the unweighted values of slope and intercept were found with the exception of intercepts in the rain gauge “a“ (intercepts remarkably lower in the case of weighting procedures in the order of +4.0).

With reference to the rain gauge “d“, i.e., that located near the main watershed divide, the R approach provided remarkably higher values of both slope (8.1) and intercept (11.5) than those obtained with the other regression approach (we recall that all values of intercepts from “d” were negatives). It must be highlighted that the standard deviations of the estimates are slightly higher than those obtained with the other regressions (see Table S2 in Supplementary Materials).

By considering the surface water, the RMA and MA approaches almost provided slightly higher values of slopes and intercepts (up to 13.0 and 55.2, respectively, in the case of river “5”). In the case of rivers “1,2,3,4,7,8”, values of slopes are in the range of those obtained by weighting procedures. On the contrary, intercepts showed a larger variability among the regression methods investigated. It should be highlighted that in the case of rivers “5,6,9” the values of slopes remarkably varied as well, in particular if MA and R were used. As in the case of the abovementioned rain gauge “d”, the δ¹⁸O–δ²H alignments from “5,6,9” were characterised by the lowest values of R² and the larger values of standard deviations of the estimates (see Table S4 in Supplementary Materials).

It must be noted that the discrepancies reported for these points (i.e., “5,6,9”) affected water with the presence of several outliers and/or serial correlations of residuals (see Figures S2 and S3 in Supplementary Materials and Section 4.1), which violated the stationary assumption.

Akin to the cases of rain gauges and rivers, RMA and MA approaches carried out on groundwater from wells and springs were characterised by larger values of both slopes and intercepts than in the case of OLS. With the exception of groundwater from “C” and “D”, the R and PW approaches induced larger variations in both slopes and intercepts, which were particularly marked in the case of groundwater from β (i.e., water whose δ¹⁸O–δ²H alignment was characterised by low values of R² and large standard deviations of the estimates, see Table S2 in Supplementary Materials).

In

Table 4. Correlation matrix reporting associations among the slopes from different regression approaches considered in this study (namely: OLS, RMA, MA, R, W, B). Progressively darker green colour is associated with a higher correlation coefficient. * Significant as p < 0.01.

	OLS	RMA	MA	R	PW	W
RMA	0.86 *
MA	0.35	0.77 *
R	0.50	0.29	−0.06
PW	0.99 *	0.87 *	0.38	0.53
W	0.87 *	0.29	−0.54	0.39	0.79 *
B	0.84 *	0.31	−0.47	0.27	0.74 *	0.97 *

, the matrix reporting correlation coefficients between pairs of slopes indicated that the largest degree of association was found (p < 0.01) between OLS–PW and W–B. High values of correlation (with slightly lower degree of association but still p < 0.01) also characterised the following relations: OLS–W, OLS–RMA, OLS–B, RMA–PW, RMA–MA, and PW–W. A significant degree of association (p < 0.01) was also found for PW–B. It should be highlighted that in several cases regarding R and MA, the degree of associations was very low and, in some cases, even negative (i.e., an increase in the value of slope obtained with a regression corresponds to a decrease in the series obtained with RMA).

By considering the intercepts (

Table 5. Correlation matrix reporting associations among the intercepts from different regression approaches considered in this study (namely: OLS, RMA, MA, R, W, B). Progressively darker green colour is associated with a higher correlation coefficient. * Significant as p < 0.01.

	OLS	RMA	MA	R	PW	W
RMA	0.84 *
MA	0.27	0.74
R	0.44	0.23	−0.13
PW	0.94 *	0.83 *	0.31	0.46
W	0.59	0.60	−0.02	0.62	0.71 *
B	0.61	0.63	−0.01	0.63 *	0.71 *	0.98 *

), the degree of associations already highlighted for slopes was further confirmed with the exception of OLS–W and OLS–B (here not significant as p > 0.01). It should be noted that almost all correlations were slightly lower than the corresponding ones from the slopes.

The hierarchic cluster analysis (see Figure S5 in Supplementary Materials) among the slope series from the several regression approaches (reported as different branches composing the dendrogram) demonstrated that MA was associated with none of the other regression methods, while two main group of pairs were clearly separated: the first is represented by OLS–PW while the second by W–B. The aforementioned first and second group were associated with each other while longer branches further linked them to R and RMA.

With reference to intercepts, the dendrogram confirmed the nonassociation of MA with the other regression approaches. Moreover, the two closest series were still those of OLS and PW, which were in turn associated to W. Contrary to the case of slopes, B series was strictly associated to RMA.

4.3. Comparison between the Differences in the Slopes and Intercepts with Catchment Characteristics and Statistics

By taking into account only δ¹⁸O–δ²H regressions from surface water, the Pearson rank correlation matrix (we recall that Pearson assesses linear relationships) comparing differences in slopes with catchment characteristics (see

Table 6. Matrix of the Pearson (in grey) and Spearman (in green) rank correlation (values as r and r_s, respectively; r value evidenced in grey while r_s in green) between the differences in the slopes and the selected catchment characteristics considered for the 9 rivers. Relationship between differences in the slopes and coefficient of determination R² from δ¹⁸O–δ²H linear regressions are also reported. R² values from regressions were calculated starting from signed values of differences in slopes. * Significant as p < 0.01.

Descriptor	OLS–RMA		OLS–MA		OLS–R		OLS–PW		OLS–W		OLS–B
Hmin (m asl)	0.03	−0.03	0.05	−0.03	0.37	0.57	0.00	−0.01	−0.41	−0.61 *	−0.40	−0.65 *
A (kmq)	−0.19	−0.12	−0.22	−0.12	−0.02	−0.01	−0.04	−0.12	0.35	0.15	0.56	0.32
Hmax (m asl)	−0.41	−0.42	−0.35	−0.42	0.04	−0.01	−0.02	−0.12	0.08	0.12	0.313	0.05
Hmean (m asl)	−0.05	−0.28	−0.02	−0.28	0.38	0.28	−0.09	−0.18	−0.14	−0.19	−0.10	−0.19
q (Ls−1/km2)	−0.57	−0.30	−0.57	−0.31	0.05	0.11	0.04	0.05	0.18	−0.05	0.12	−0.03
q95(Ls−1/km2)	−0.14	−0.40	−0.10	−0.40	0.24	0.05	−0.21	−0.29	−0.04	−0.10	0.00	−0.04
P (mm)	0.09	0.01	0.08	0.01	0.00	−0.01	0.00	0.01	−0.15	−0.11	−0.14	−0.02
F (km)	−0.11	−0.05	−0.13	−0.05	−0.23	−0.25	−0.01	−0.02	0.26	0.34	0.42	0.56
Fyw (%)	0.41	0.70 *	0.32	0.70 *	−0.09	−0.10	0.00	0.03	0.00	0.00	0.00	0.01
R2	0.98 *	0.98 *	0.96 *	0.96 *	0.01	0.02	0.10	0.18	−0.46	−0.03	−0.32	−0.06
δ18O range	0.16	0.45	0.13	0.46	−0.38	−0.27	−0.04	0.00	0.05	0.29	0.15	0.19
δ2H range	0.00	0.03	−0.26	0.03	0.25	−0.28	−0.34	−0.34	0.25	0.46	0.36	0.45
n°of samples	0.04	0.03	0.04	0.03	0.17	0.12	−0.35	−0.27	−0.03	0.00	−0.03	−0.01

) did not provide significant (p < 0.01) correlations. If the Spearman rank correlation matrix (i.e., assessment of nonlinear relationships) was considered, we found positive and significant (p < 0.01) correlations with F_yw (Slope_OLS–RMA and Slope_OLS–MA) while negative ones with H_min (Slope_OLS–W, Slope_OLS–B).

In the case of the Pearson rank correlation matrix applied to differences in intercepts, we did not find significant correlations (see

Table 7. Matrix of the Pearson (in grey) and Spearman (in green) rank correlation (values as r and r_s, respectively; r value evidenced in grey while r_s in green) between the differences in the intercepts and the selected catchment characteristics considered for the 9 rivers. Relationship between differences in the slopes and coefficient of determination R² from δ¹⁸O–δ²H linear regressions are also reported. R² values from regressions were calculated starting from signed values of differences in intercepts. * Significant as p < 0.01.

Descriptor	OLS–RMA		OLS–MA		OLS–R		OLS–PW		OLS–W		OLS–B
Hmin (m asl)	0.04	−0.03	0.05	−0.03	0.36	0.56	0.00	−0.01	−0.40	−0.53	−0.48	−0.62 *
A (kmq)	−0.19	−0.12	−0.22	−0.12	−0.02	−0.01	−0.04	−0.07	0.35	0.14	0.42	0.28
Hmax (m asl)	−0.38	−0.42	−0.35	−0.42	0.04	−0.01	−0.02	−0.11	0.10	0.16	0.04	0.06
Hmean (m asl)	−0.04	−0.28	−0.02	−0.28	0.36	0.28	−0.08	−0.20	−0.12	−0.14	−0.18	−0.18
q (Ls−1/km2)	−0.53	−0.30	−0.56	−0.30	0.05	0.11	0.04	0.06	0.20	−0.05	0.10	−0.03
q95(Ls−1/km2)	−0.12	−0.40	−0.08	−0.40	0.24	0.05	−0.20	−0.28	−0.02	−0.06	−0.05	−0.03
P (mm)	0.10	0.01	0.09	0.01	0.00	−0.01	0.00	0.01	−0.13	−0.12	−0.08	−0.02
F (km)	−0.32	−0.05	−0.12	−0.05	−0.22	−0.25	−0.01	0.00	0.27	0.30	0.46	0.53
Fyw (%)	0.36	0.70 *	0.30	0.70 *	−0.09	−0.10	0.00	0.01	0.01	0.00	0.00	0.00
R2	0.98 *	0.98 *	0.96 *	0.96 *	0.01	0.00	0.02	0.15	−0.57	−0.12	−0.67 *	−0.21
δ18O range	0.13	0.45	0.12	0.45	−0.41	−0.27	−0.04	0.01	0.03	0.27	0.14	0.18
δ2H range	0.00	0.03	0.00	0.03	−0.30	−0.28	−0.34	−0.23	0.22	0.45	0.34	0.46
n°of samples	0.04	0.03	0.04	0.03	0.15	0.12	0.31	0.27	0.03	0.00	0.03	0.01

). On the contrary, and as in the case of slopes, significant (p < 0.01) nonlinear relationships were found for F_yw (Slope_OLS–RMA and Slope_OLS–MA) and H_min (negative correlation for Slope_OLS–B).

In all cases, the largest statistical performance (again significant as p < 0.01) was found for coefficient of determination R² from δ¹⁸O–δ²H linear regressions (Slope_OLS–RMA and Slope_OLS–MA; Intercept_OLS–RMA and Intercept_OLS–MA).

5. Discussion

We did not find significant discrepancies in the slopes and intercepts computed by the different regression methods in the case of precipitation data. On the contrary, marked variations were detected in the case of river water and groundwater (both from springs and wells in lowland aquifers) obtained using specific methods. Among others, such discrepancies were somehow reduced in the case of OLS, RMA, and PW. Because of different values of river and spring discharges and the corresponding changes in isotopic content of water during the year, weighting procedures (W, B) were characterised by diverse values of slopes and intercepts rather than the aforementioned OLS, RMA, and PW. Moreover, and as highlighted by both heat maps (Figure S3 in Supplementary Materials) and correlation matrices (Table 4 and Table 5) and dendrograms (Figure S5 in Supplementary Materials), slopes and intercepts from the MA and R approaches were not comparable to others (nonsignificant statistical associations). Regardless of water type, the aforementioned discrepancies were promoted when δ¹⁸O–δ²H regressions were characterised by weak statistical performances (low values of R² and larger values of standard deviations of the estimates). With reference to rivers, the weak statistical performances were linked to the presence of outliers and/or serial correlation of the residuals violating the stationary assumption of OLS, MA, and RMA approaches.

The investigation carried out on the data solely from rivers highlighted that the magnitude of the differences in the slopes and intercepts was related in all cases (with the exception of R and PW) to the coefficient of determination R² characterising δ¹⁸O–δ²H linear regressions. The largest values of Pearson coefficients (see Table 6 and Table 7) led us to consider R² as the main causal factor for such differences in slopes and intercepts.

In particular, the larger the correlations between δ¹⁸O and δ²H, the smaller the differences among slopes and intercepts detected by RMA, MA, W, and B within the specific sampling point (river, well, or spring). This is in agreement with the results reported by [15] and corroborated the hypothesis that statistical performance of the regression was the main driver of these slope and intercept variations. In any case, despite finding highly statistical significance with R² in our investigated dataset, no relations between differences of slopes (and intercepts) and the ranges in δ¹⁸O (and δ²H) along with the number of samples composing the dataset were noticed, thus indicating that extreme values of δ¹⁸O (and δ²H) were not significant causal factors.

With reference to RMA and MA, the Spearman rank correlation matrices involving differences in slopes and intercepts and catchment descriptors allowed us to find a significant nonlinear association with F_yw (we recall here that F_yw is the percentage of water younger than 2–3 months). In both cases (RMA and MA) the association (reported also as plots in Figure 2) indicated that the magnitude of the differences in the slopes and in the intercepts decreased along with the quota of young water.

This means that rivers showing low values of F_yw are likely to be more affected by differences in slopes and intercepts computed by different regression approaches. By examining the plots reported in Figure 2, it can be evidenced that nonlinearity is driven by two catchments (namely, the Secchia River “5” and Panaro River “6”). As already anticipated in Section 2, these two rivers (“5,6”) were the only ones characterised by nival–pluvial discharges due to the melting of the snow cover in the upper part of the catchments during the spring months. Moreover, [7] stated that there were evidences of sublimation in several water samples collected in rivers “5,6” from January 2017 to April 2017. This was further confirmed by the remarkable number of snowfall events that occurred between December 2013 and April 2014 over the highest part of the catchments “5,6”. Such events have allowed the consequent snowpack development alternating with partial snowmelt for a snow water equivalent higher than 600 mm [21].

There, we recall that sublimation occurring during sunny days can modify the former isotopic composition of the superficial snow layers, allowing the release of a vapour phase from the solid skeleton to the atmosphere. In this case, the final snow cover does not preserve the isotopic composition of the original snowfall from which it was derived [44,45], a fact that also led differences in slopes and intercepts from δ¹⁸O–δ²H regressions to be enhanced. In detail, sublimation acting on a snow cover can lead to an enrichment of heaviest isotopes (such as ¹⁸O and ²H) within the solid skeleton and can induce a similar δ¹⁸O–δ²H pattern of that charactering the residual liquid subjected to evaporation (slope decrease of the δ¹⁸O–δ²H alignments for snowpack samples if compared to the water meteoric line (see for field and experimental studies: [46,47])). In particular, the slope decrease can be much more intense if the only late-season snowpack samples are considered (a value of 3.7 was found by [48]).

In case the two rivers “5,6” are removed from the analysis, it is still possible to confirm such alignments, although linear, between y F_yw and differences in slope and intercept pairs. In this sense, such relations, still identifying an inverse association between differences in the slopes (and in the intercepts) and quotas of young water, may also be related to other hydrological processes taking place at the catchment scale. As already pointed out by [7], by checking both slopes (river water showed slightly lower values than those characterising rainwater; see Figure S6 in Supplementary Materials) and intercepts (that were negative compared to those from rainwater), all the river water considered underwent evaporation/evapotranspiration processes prior to their infiltration towards the aquifer.

Albeit to a lesser extent, these variations also affected groundwater from wells (which were also fed by streambed dispersion and therefore by water isotopes already modified in the river water; see [11] and low-yield springs (these potentially characterised by pre-infiltrative modification as slopes from δ¹⁸O–δ²H alignments were slightly lower than those obtained from precipitation water; see Figure S7 in Supplementary Materials)). On the contrary, the nonvariability of slopes and intercepts observed in the different δ¹⁸O–δ²H alignments from precipitation was somehow expected, as these waters were unlikely to be affected by evaporative/sublimation processes once they had entered the rain gauges.

With reference to the effective role of young water fraction F_yw in influencing the differences in intercepts and slopes, we believe that further efforts have to be made for such catchments from the hilly part of the northern Italian Apennines and dominated by higher quotas of young water (i.e., F_yw > 25%; not analysed in this study for lack of isotopic data). The latter are characterised by intermittent discharge and wide outcrops of low permeable soils and bedrocks (prevailing clayey and marls materials; G_C and G_M in Figure 1). Further investigations should be isotopic-based in order to verify also the role of pre-infiltrative evaporation in isotopic deviation and, above all, in the change of slopes and intercepts.

Following [7], these processes were likely to be promoted in the clay-rich bedrock, where water molecules composing the soil moisture were slowed in percolation and thus kinetic fractionation processes were enhanced.

However, we can provide some preliminary recommendations for use of the different regression approaches for the four water types (precipitation, surface water, groundwater from wells, and low-yield springs) from the northern Italian Apennines:

(i)

In the case of δ¹⁸O–δ²H alignments from precipitation, and as no remarkable discrepancies were detected among the several investigated methods, the OLS approach should be preferred.

(ii)

For precipitation and surface water, slopes and intercepts from the two weighting procedures W and B were similar. Moreover, there was no evidence of remarkable changes among results obtained from such weighting procedures with those from unweighted OLS. The latter confirms the convenience of using the OLS approach even if, during the year, rainfall or discharge amounts (and isotopic content too) are different between the seasons.

(iii)

For surface water and groundwater, the MA and R approaches should not be used in any case as they seem to provide unrealistic values for both slopes and intercepts. The reason has to be searched in the fact that these two approaches are more sensitive to the statistical performance of the regressions (i.e., standard deviations), especially if outliers are present. MA demonstrated to be more sensitive to the statistical performance of the regressions (i.e., standard deviations), especially if outliers are present. Although the R approach was selected to verify its behaviour in the case of outliers, it did not induce improvements in standardised residuals. Moreover, it was also demonstrated that kinetic fractionation processes acting on these water types lead to increase the differences in slopes and intercepts (see, for instance, relationships between differences in intercepts and slopes with young water fraction Fyw). Slopes and intercepts from OLS and PW were the closest, with lower standard deviations sometimes associated to PW regressions. In addition, and with reference to the surface water, PW results were not affected by the kinetic fractionation processes (see Table 6 and Table 7).

(iv)

Surface water may be affected by nonstationary processes induced by both nonmultivariate normality and serial correlations of the residuals. Thus, prior to carrying out OLS regression on δ¹⁸O–δ²H data from surface water (and groundwater), the presences of outliers, heteroscedasticity, and autocorrelation must be carefully detected by means of both conventional statistical tests and inspections of standardised residuals. In the case of outliers, their importance on the whole data series composing the regression should be evaluated (as an instance, in the case of δ¹⁸O–δ²H pairs from surface water collected during the late summer through the beginning of the autumn period, the strong reduction in discharges may induce their removal from the dataset prior to carrying out the regression). In the event of dealing with time series of stable isotopes affected by autocorrelation, we believe it is convenient to use the PW approach, which, in our case, has proven to solve the serial correlations of residuals.

6. Conclusions

We presented the comparison of five different regression approaches applied to δ¹⁸O–δ²H data from four different water types collected in the northern Italian Apennines. We found that all the tested approaches converged towards similar values of slopes and intercepts for only stable water isotopes from precipitation. Conversely, differences in slopes and intercepts from surface water and groundwater (collected from wells and low-yield springs) were often significant and related to the robustness of the regressions (i.e., standard deviations of the estimates) and their sensitiveness to outliers and autocorrelation. Moreover, and with reference to surface water, we found evidence of a relationship between young water fraction and the magnitudes in differences of slopes and intercepts, suggesting the control of kinetic fractionation processes (mainly related to sublimation acting on snow cover and, secondary, to active pre-infiltrative evaporation and evapotranspiration processes) on such discrepancies. These results allowed us to provide some recommendations for hydrological and hydrogeological studies involving δ¹⁸O–δ²H from the abovementioned water types collected in the northern Italian Apennines. Firstly, as no discrepancies were noticed between slopes and intercepts from all the methods applied to precipitation, the OLS approach is preferred. Secondly, and with reference to the other water types (surface water and groundwater from wells and springs), we warmly suggest carrying out conventional statistical tests coupled with inspection of standardised residuals for a preliminary check on the presence of outliers and autocorrelation phenomena. In the case of managing outliers, the MA and R approaches should be avoided as they are more sensitive to the statistical performance of the regressions and often provide unrealistic values of both slopes and intercepts. Thirdly, for surface water and groundwater, the OLS and PW approaches still showed the highest degree of robustness and produced the closest values of slopes and intercepts, thus resulting as the methods preferable for δ¹⁸O–δ²H regressions. PW would be more reliable in the presence of serial correlations of the residuals (which, in our case, often affected surface water). In the case of managing outliers, the possibility of removing them will have to be considered (as an example in the case of δ¹⁸O–δ²H pairs from marked low-flow periods).

Lastly, despite the presence of marked differences in the amounts of rainfall and their isotopic contents during the year, the convenience of using weighing approaches before applying OLS was not found.