Local Meteoric Water Line of Northern Chile (18° S–30° S): An Application of Error-in-Variables Regression to the Oxygen and Hydrogen Stable Isotope Ratio of Precipitation

Abstract

In this study, a revision of the previously published data on hydrogen (²H/¹H) and oxygen (¹⁸O/¹⁶O) stable isotope ratio of precipitation in northern Chile is presented. Using the amount-weighted mean data and the combined standard deviation (related to both the weighted mean calculation and the spectrometric measurement), the equation of the local meteoric line calculated by error-in-variables regression is as follows: Northern Chile EIV-LMWL: δ²H = [(7.93 ± 0.15) δ¹⁸O] + [12.3 ± 2.1]. The slope is similar to that obtained by ordinary least square regression or other types of regression methods, whether weighted or not (e.g., reduced major axis or major axis) by the amount of precipitation. However, the error-in-variables regression is more accurate and suitable than ordinary least square regression (and other types of regression models) where statistical assumptions (i.e., no measurement errors in the x-axis) are violated. A generalized interval of δ²H = ±13.1‰ is also proposed to be used with the local meteoric line. This combines the confidence and prediction intervals around the regression line and appears to be a valid tool for distinguishing outliers or water samples with an isotope composition significantly different from local precipitation. The applicative examples for the Pampa del Tamarugal aquifer system, snow samples and the local geothermal waters are discussed.

1. Introduction

To fit a model to data pairs, the linear least squares regression is by far the most widely used modeling method (e.g., [1]). The ordinary least squares regression (OLSR) minimizes the sum of the squared vertical distances between the y data values and the corresponding y values on the fitted line (the predictions). The OLSR design assumes that there is no variation in the independent variable (x) because it is controlled by the researcher. This is also known as ‘Model I’ regression [2]. In contrast, in ‘Model II’ regression, both the response and the explanatory variables of the model are random (i.e., not controlled by the researcher), and there are errors associated with the measurements of both x and y [2]. Among the ‘Model II’ regression, major axis regression (MA) assumes equal variance on both variables and minimizes the orthogonal (perpendicular) distances from the data points to the fitted line. To the best of our knowledge, one of the first suggestions for applying orthogonal regression to the precipitation amount weighted stable isotope ratios of the oxygen (δ¹⁸O, the independent variable x) and hydrogen (δ²H, the dependent variable y) of water molecules in order to determine a better global meteoric water line (GMWL) dates back to the 1980s [3]. Ten years later, the International Atomic Energy Agency (IAEA) proposed the use of the reduced major axis (RMA) regression for calculation of the local meteoric water lines (LMWL) [4]. However, in this latter publication, RMA was wrongly defined as an orthogonal regression. In fact, RMA minimizes the sum of the areas (thus using both vertical and horizontal distances of the data points from the resulting line) rather than the least squares sum of the squared vertical distances, as in OLSR [5]. Meanwhile, the local meteoric water lines in the GNIP dataset of IAEA are calculated through precipitation amount-weighted least squares regression (PWLSR) [6,7]. This regression approach is considered to be the most suitable for representing the isotope composition of groundwater because these are recharged by important rainfall events [6,7]. Alternatively, the precipitation weighted RMA (PWRMA) appears to be most suitable for coastal, island, and Mediterranean sites [6]. Another proposed approach for the calculation of LMWL was the generalized least squares regression model (GENLS), an error-in-variables regression (EIV) that considers the combined standard deviation of the δ¹⁸O and δ²H values [8]. Within the statistical literature, synonyms of this approach are “Deming” [9] or “error-in-variables” [10] regression, the former often distinguished as simple (when the data errors are constant among all measurements for each of the two variables) and general (different data error at each observation) [11]. It should be noted that in the case where the variance ratio is equal to 1, Deming regression is equivalent to orthogonal regression. From a general statistical and chemometrics point of view, comparisons between the different approaches can be found in the existing literature [12]. In contrast, in the specific case of water isotope geochemistry and the related meteoric water line, a comparison between error-in-variables and other regression methods has never been presented or discussed. Moreover, despite the advantages identified by previous publications of the alternative regression methods [4,6,7,13], most of the studies on the stable isotope ratio of the waters still use OLSR to calculate the LMWLs. While this does not necessarily mean that the use of OLSR is wrong [7], it has been shown, especially in the hydrology study of arid regions, that the differences between the OLSR approach and other regression approaches can be significant, particularly in terms of the slope of the LMWL [7].

In this manuscript, we compare the regression lines obtained from OLSR, RMA, MA, PWLSR, PWRMA, and PWMA with that from the error-in-variables (EIV) approach, which is very similar to the method followed for GENLS [8]. Here, the existing data on δ¹⁸O and δ²H values of precipitation in northern Chile were chosen. This brief study does not attempt to be mathematically or statistically exhaustive and the reader is advised to refer to the referenced literature for more details. In fact, the main purpose here is to provide the results of alternative regression methods applied to stable isotope ratios of precipitation and to provide guidance and advice when the obtained local meteoric water lines are employed for the interpretation of the isotope data of groundwater.

2. Materials and Methods

Northern Chile extends southwards from approximately 18° S to 30° S and includes the two macro-regions of Norte Grande (Regions XV, I and II; Figure 1B) and part of Norte Chico (Regions III and IV; Figure 1B) [14,15,16].

This presents a narrow strip of land (~300 km) between the Andean ridge and the Pacific coast [17]. The latitudinal range of this northern part of Chile aligns with a major extension of the arid desert climate (BW climate zone [20,21]), followed by a tundra climate (ET climate zone [20,21]), and a semi-arid climate in a minor extension [14,21,22]. The Atacama Desert, the driest place on earth with mean precipitation as low as 10 mm per decade, and the higher elevations in the Andes are striking examples of the BW and the ET climates, respectively. The Amazonian influence in the South American monsoon system, with concentrated summer rainfall or dry winters, affects the climate in this area from the north boundary to Nevado Ojos del Salado [22], a stratovolcano with an altitude of 6893 m, the highest point in Chile and the second highest outside of Asia.

The previous studies on isotope composition of meteoric water revealed a general increase in δ¹⁸O on the Andean plateau (the “Altiplano”, Figure 1C) from north to south, concomitant with an increase in aridity and decrease in convective moistening (amount effect; [23]). Stable isotope and seasonal precipitation patterns suggest an eastern provenance of the vast majority of moisture that falls as precipitation across the Andean Plateau and Western Cordillera (Figure 1C), with Pacific-derived moisture contributing a minor amount at low elevations near the coast ([23,24]). However, over most regions, the δ¹⁸O signal of precipitation is influenced by a combination of factors ([24]). The δ¹⁸O -depleted values observed in the high altitude area (i.e., the Andean plateau) were related to processes that affect the air masses that (i) originated over the Atlantic Ocean, (ii) cross the Amazon Basin (continental effect), (iii) ascended the Andes (altitude effect) and (iv) precipitated (convective effect) in the Andean plateau ([25]). In particular, over the eastern Andes, precipitation at low elevations has δ¹⁸O from −2 to −8‰, but δ¹⁸O becomes more depleted toward the west as vapor is lifted across the Eastern Cordillera of the Andes ([26,27,28,29]). The dominant Altiplano summer rain has δ¹⁸O values from −8‰ to −15‰, with the expectation that precipitation in the driest places (e.g., Atamaca desert) should be moderately to strongly depleted at all elevations up to approximately −18‰ ([28,30]).

During the last forty years, several meteoric water lines calculated through OLSR regression were proposed to describe the best fit for the oxygen and hydrogen stable isotope ratios of the precipitation in northern Chile (

Table 1. Slope and intercept values of the previously published northern Chile meteoric water lines [31,32,33,34,35,36,37,38,39]. All the shown data were calculated by OLSR. The IAEA/GNIP data pairs are referred to as the “La Serena” station (Figure 1B). For this station, other two equations are also proposed on the same database, δ²H = [(7.64 ± 0.39) δ¹⁸O] + [9.44 ± 2.32] and δ²H = [(8.04 ± 0.45) δ¹⁸O] + [9.98 ± 2.59] by PWLSR and RMA methods [6,7], respectively, both applied to all the available data (N = 40).

	Slope	±	s.e.	Intercept	±	s.e.	N	R2	Reference	Year	Ref. [#]
	7.8		-	10.3		-	29	-	Fritz et al.	1981	[35]
	8.2	±	0.17	16.6	±	2.2	39	0.98	Chaffaut et al.	1998	[32,33]
	7.8		-	9.7		-	129	-	Aravena et al.	1999	[31]
	7.9		-	14		-	-	-	Herrera et al.	2006	[36]
	7.7		-	9.6		-	-	0.91	Squeo et al.	2006	[38]
	7.95		-	14.9		-	-	0.99	DGA	2015	[34]
	8.07		-	13.5		-	23	0.98	Troncoso et al.	2012	[39]
	7.52	±	0.46	7.18	±	2.64	40	0.87	IAEA/WHO	2015	[37]
mean	7.87			11.97
std.de.	0.21			3.22

-: data not calculated or not declared; s.e.: standard error.

, [31,32,33,34,35,36,37,38,39]).

However, not all the studies show the oxygen and hydrogen stable isotope ratios of the water molecule along with precipitation amount [31,32,35,40,41,42]. “La Serena” is the northern Chile station (Region IV, Figure 1B) monitored by the Global Network of Isotopes in Precipitation (GNIP dataset) from 1988 until 2015 [37].

For this study, a dataset of 32 stations has been constructed through collecting previously published data, considering only the isotope values of precipitation with the measured amount and recalculating the amount-weighted mean in order to have only one δ¹⁸O and δ²H data pair per station (

Table 2. Amount-weighted oxygen and hydrogen stable isotope ratios (δ¹⁸O and δ²H in permil versus the standard mean ocean water SMOW, respectively) of the previously published precipitation samples collected in northern Chile [31,32,35,37,40,41,42,43]. Calculations are described in the Supplementary File S1. Samples were mainly collected in Region I (from #1 to #27), Region II (from #28 to #31) and Region IV (#32) (Figure 1B).

ID #	Location	Coordinates		Elevation Meters	δ18O (a.w.m. ± c.s.)	δ2H (a.w.m. ± c.s.)	Na.p.	Precipitation Amount (Mean Value) Millimeters	References	Ref. [#]
		Latitude	Longitude		‰ vs. SMOW	‰ vs. SMOW
1	Apacheta Tapa	−19.5891	−68.9940	4350	−17.48 ± 1.14	−122.3 ± 8.9	4	73.73	Fritz et al. 1981	[35]
2	Chusmiza	−19.7060	−69.1906	3360	−7.51 ± 0.22	−48.9 ± 1.0	3	103.93	Fritz et al. 1981	[35]
3	Alto Mocha	−19.8181	−69.2987	2590	−7.13 ± 0.31	−43.7 ± 1.5	2	33.90	Fritz et al. 1981	[35]
4	Mocha	−19.8436	−69.2840	2200	−4.91 ± 0.90	−29.2 ± 6.7	2	39.65	Fritz et al. 1981	[35]
5	Collacagua	−20.0703	−68.8259	3915	−15.00 ± 1.22	−108.1 ± 8.6	7	82.23	Fritz et al. 1981	[35]
6	Huasco	−20.3144	−68.8065	3800	−17.13 ± 0.65	−127.7 ± 9.3	3	34.63	Fritz et al. 1981	[35]
7	Indio Muerto	−20.3610	−68.9615	4135	−18.27 ± 0.35	−129.9 ± 3.4	4	53.05	Fritz et al. 1981	[35]
8	Tambillos	−20.4672	−69.1415	3300	−8.02 ± 0.27	−50.1 ± 2.1	3	77.13	Fritz et al. 1981	[35]
9	Apacheta Mama *	−20.4734	−69.1580	4115	−19.90 ± 0.15	−146.0 ± 1.0	1	77.00	Fritz et al. 1981	[35]
10	Pumire	−19.0955	−69.1108	4200	−2.25 ± 0.78	−11.5 ± 3.1	3	10.00	Salazar et al. 1998; Aravena et al. 1989	[41,42]
11	Coposa	−20.7089	−68.6942	3460	−16.66 ± 2.01	−119.8 ± 16.2	9	9.17	Salazar et al. 1998; Aravena et al. 1989	[41,42]
12	Collaguasi	−20.9833	−68.7000	4250	−16.96 ± 1.99	−122.3 ± 16.0	5	13.30	Aravena et al. 1999	[31]
13	Ujina	−20.9833	−68.6000	4200	−11.98 ± 1.57	−83.3 ± 12.0	8	18.79	Aravena et al. 1999	[31]
14	Puchuldiza	−19.4000	−68.9500	4150	−14.55 ± 0.93	−98.4 ± 11.0	15	11.88	Aravena et al. 1999	[31]
15	Pampa Lirima	−19.8200	−68.9000	4100	−16.97 ± 1.23	−121.6 ± 9.7	27	9.86	Aravena et al. 1999	[31]
16	Colchane	−19.7003	−68.8833	3965	−14.95 ± 2.28	−106.3 ± 18.2	4	51.90	Aravena et al. 1999	[31]
17	Collacagua	−20.0333	−68.8500	3990	−12.68 ± 1.34	−88.3 ± 10.5	14	15.04	Aravena et al. 1999	[31]
18	Cancosa	−19.8500	−68.6000	3800	−14.00 ± 2.45	−96.4 ± 20.7	5	28.50	Aravena et al. 1999	[31]
19	Huaytane	−19.5433	−68.6042	3720	−14.10 ± 1.48	−100.9 ± 12.1	5	32.20	Aravena et al. 1999	[31]
20	Copaquire	−20.9305	−68.8922	3490	−12.47 ± 0.98	−84.6 ± 7.8	16	5.75	Aravena et al. 1999	[31]
21	Poroma	−19.8667	−69.1833	2880	−5.00 ± 0.66	−28.1 ± 6.0	9	8.01	Aravena et al. 1999	[31]
22	Parca	−20.0167	−68.8500	2570	−6.64 ± 0.90	−40.9 ± 7.8	4	4.88	Aravena et al. 1999	[31]
23	Huatacondo	−20.9167	−69.0500	2460	−8.26 ± 1.26	−53.8 ± 9.8	12	6.21	Aravena et al. 1999	[31]
24	Camina	−19.3116	−69.4299	2380	−7.02 ± 1.00	−42.8 ± 8.1	4	11.50	Aravena et al. 1999	[31]
25	Sillillica	−20.1738	−68.7377	4270	−15.96 ± 1.07	−110.9 ± 7.9	3	95.30	Uribe et al. 2015	[40]
26	Altos del Huasco	−20.3221	−68.9022	3784	−14.82 ± 1.57	−104.8 ± 12.3	3	51.80	Uribe et al. 2015	[40]
27	Diablo Marca	−20.0505	−68.9962	4603	−18.43 ± 1.59	−128.8 ± 11.8	3	78.60	Uribe et al. 2015	[40]
28	Quisquiro	−23.2100	−67.2500	4260	−13.40 ± 2.12	−91.1 ± 13.3	36	14.30	Chaffaut 1998	[32]
29	Tuyajto	−23.9435	−67.5921	4040	−6.23 ± 0.50	−32.4 ± 3.1	2	13.30	Herrera et al. 2016	[43]
30	Pampa Colorada	−23.8588	−67.4774	4426	−8.80 ± 0.86	−51.3 ± 7.6	3	28.50	Herrera et al. 2016	[43]
31	Aguas Calientes-3	−23.9171	−67.6971	3900	−6.00 ± 0.51	−31.9 ± 4.5	3	21.67	Herrera et al. 2016	[43]
32	La Serena	−29.8981	−71.2425	142	−5.78 ± 1.32	−35.9 ± 9.9	34	43.10	IAEA/WHO 2015	[37]

a.w.m. ± c.s.: amount-weighted mean ± combined sigma (see methods paragraph for details); N_a.p.: Number of accumulation periods of precipitation. This value coincides with the total number of the available isotope data pairs per station (Supplementary File S1); * sample with only one accumulation period (N_a.p. = 1). In this case, the combined sigma (c.s.) contain only the standard deviation related to spectrometric measurement (Supplementary File S1). The regression has been calculated both including (Table 3) and excluding (Supplementary File S2) this sample.

; Supplementary File S1).

Assuming ${\mathrm{x}}_{\mathrm{i}}$ to be the measurement of the stable isotope ratios expressed as δ¹⁸O and δ²H in ‰ versus SMOW (standard mean ocean water)—determined on a specific precipitation amount ${\mathrm{p}}_{\mathrm{i}}$ (mm) collected over a specific time period—the standard deviation related to the amount weighted mean ${\overline{\mathrm{x}}}_{\mathrm{w}}$ [4] of the meteoric stations with more than one accumulation periods of precipitation (${\mathrm{N}}_{a.p.}$ > 1) was calculated as follows [44]:

$${\mathsf{\sigma }}_{\mathrm{w}}=\sqrt{\frac{{\displaystyle \sum _{\mathrm{i}=1}^{{\mathrm{N}}_{a.p.}}{\mathrm{p}}_{\mathrm{i}}{\left({\mathrm{x}}_{\mathrm{i}}-{\overline{\mathrm{x}}}_{\mathrm{w}}\right)}^2}}{\frac{\left({\mathrm{N}}_{a.p.}-1\right){\displaystyle \sum _{\mathrm{i}=1}^{{\mathrm{N}}_{a.p.}}{\mathrm{p}}_{\mathrm{i}}}}{{\mathrm{N}}_{a.p.}}}}$$

(1)

It should be noted that the accumulation periods, which coincides with the total number of the isotope data pairs and precipitation amounts available per each station, are not uniformly distributed through the year and not equally long for most of the stations considered (Supplementary File S1).

In the existing literature, there are several definitions of ‘combined standard deviation’ or ‘combined uncertainty.’ Here, we deal with two different ‘uncertainty’ sources [8,45], that is, a standard deviation related to the amount-weighted mean calculation ${\mathsf{\sigma }}_{\mathrm{w}}$, as described above, and a standard deviation related to the mass spectrometric measurements ${\mathsf{\sigma }}_{\mathrm{ms}}$, which differs depending on the publication (Supplementary File S1). One of the easiest ways to calculate the combined standard deviation ${\mathsf{\sigma }}_{\mathrm{csd}}$ involves foreseeing the combination of ${\mathsf{\sigma }}_{\mathrm{w}}$ and ${\mathsf{\sigma }}_{\mathrm{ms}}$ following the law of propagation of uncertainty, that is, the square root of the sum-of-the-squares [45,46,47]:

$${\mathsf{\sigma }}_{\mathrm{csd}}=\sqrt{{\left({\mathsf{\sigma }}_{\mathrm{w}}\right)}^2+{\left({\mathsf{\sigma }}_{\mathrm{ms}}\right)}^2}$$

(2)

The combined standard deviations, associated with the δ¹⁸O and δ²H amount-weighted mean data in 30 of the total 32 stations, were used to calculate the local meteoric water line through EIV regression. In two stations, Quisquiro [32] and La Serena [37], the calculation of the standard deviation of the weighted mean ${\mathsf{\sigma }}_{\mathrm{w}}$ was performed differently for two reasons. In the case of Quisquiro, the complete dataset is unavailable; therefore, the standard deviation was calculated among the amount-weighted data corresponding to the extreme seasons (winter and summer, [32]). Meanwhile, in the case of La Serena, the ${\mathsf{\sigma }}_{\mathrm{w}}$ was calculated among the years in which more than 70% of precipitation was analyzed for a given isotope composition [4]. Following this, the ${\mathsf{\sigma }}_{\mathrm{csd}}$ was calculated through using the ${\mathsf{\sigma }}_{\mathrm{ms}}$ declared in the analytical section of the publication dataset [32] or elsewhere [25]. With this, the ${\mathsf{\sigma }}_{\mathrm{csd}}$ on δ¹⁸O and δ²H is within the mean of the other stations.

The calculation of the EIV regression through different codes employ similar methods based on York and Williamson’s works [48,49,50,51]. Generally speaking, the process of searching for best straight-line parameters requires an iterative procedure in which each new slope estimate is used to modify a weight term, which combines the reciprocal variance of the x and y variables, in a series of convergent fitting cycles [11,52,53,54,55]. The results of this method are compared with those obtained by OLSR, RMA, MA, PWLSR, PWRMA, and PWMA [6]. In accordance with Crawford et al. [6], the Student’s t-test was used to check the difference between the two regression slopes. Furthermore, the obtained EIV meteoric water line was also compared with previously published isotope data on precipitation without amount data [28,30,36,56,57,58,59,60], snow [43,58,59,61,62,63,64,65,66,67,68], groundwater [69], and geothermal water [70] samples collected up to this point in northern Chile.

3. Results

The parameters of the northern Chile meteoric water lines obtained through the different software and regression methods are listed in Table 3.

The values of slope and intercept through EIV appear to be slightly higher than those obtained by other methods that do not take into account the combined standard deviation of the oxygen and hydrogen measurements. However, there is no difference with other regression methods, aside from the results from Cantrell’s bivariate worksheet (EIV-d, Table 3b), which was tested on atmospheric chemical data [54]. The regression results of EIV-d method include: (i) a lower standard error in comparison with those from other codes used for EIV calculation; and (ii) a higher slope when compared to PWLSR, but not significantly different (p > 0.05). The difference between the slopes are further smoothed when the sample #9 is removed (N_a.p. = 1; Supplementary File S2). It should be noted that lower standard errors on slope and intercept similar to those of Cantrell’s bivariate worksheet are obtained in the other EIV models introducing a scale factor in the maximum likelihood method, which is estimated by the square root of the reduced chi-square statistic that results from the fit analysis [11,53]. Lower than half standard error are also obtained by “Deming regression” with a constant variance ratio of λ = δ¹⁸O_var/δ²H_var = 0.019 (EIV-e, Table 3b), that is the mean ratio of the squared combined sigma (Table 2). Therefore, we consider as most representative of the samples the regression results with the standard errors on slope and intercept calculated without scale factor and with combined sigma for each station (EIV results from “a” to “c” in Table 3b):

Northern Chile EIV-LMWL: δ²H = [(7.93 ± 0.15) δ¹⁸O] + [12.3 ± 2.1]

(3)

The upper and lower confidence intervals at 95% probability of the slope and intercept in the EIV regressions (results from “a” to “c” in Table 3b) include all the previous published values obtained by OLSR (Table 1) and all the other slopes and intercepts of the meteoric water lines calculated by different regression methods (Table 3a). This is clear and straightforward, in particular looking to the standard errors of slope and intercept, which are higher in EIV regression results.

4. Discussion

In isotope hydrology applications, it is often desirable to know whether a water sample from the surface (river, lake, or lagoon) or from an aquifer is shifted in a statistically significant way in relation to the local meteoric water line. This is particularly true in the following: (i) arid zones, where fractionation due to evaporation produces isotopically enriched residual waters; and (ii) geothermal areas, where hot water samples could be isotopically enriched due to water-rock interaction. Northern Chile offers an ideal location for this kind of study because both of these conditions exist here [65,72].

The calculated confidence intervals obtained through EIV regression represent the scatter of the weighted mean data around the fitted line, disregarding the combined standard deviation (Table 3). In OLSR, prediction intervals for y at a specific x value estimate the bands around the fitted line in which the future observation will fall, with a certain probability (usually 95% as confidence interval), given what has already been observed with the sample dataset employed for the regression (t-distribution multiplied by the standard error of the fit). The OLSR prediction interval can be quickly calculated by using commercial or free codes (e.g., [71,73]). However, OLSR does not take into account the measurement error in x when predicting y, either in the combined standard deviations of the data in the fittings or the prediction interval calculations. Moreover, an officially recognized method for the calculation of prediction intervals related to EIV regression still does not exist.

An approximate prediction area around the EIV local meteoric water line could be traced through using an empirical-graphical method, joining the outer limits of the bars representing the combined standard deviations around the mean of the isotope data of precipitation. However, its irregular shape, related to the non-homogeneous values of the combined standard deviations of the regressed stations, necessitates an ad-hoc correction. Specifically speaking, the lower sector of the line (approximately, δ¹⁸O < −12‰ and δ²H < −80‰) shows numerous monitored stations with wider standard deviations—especially in terms of hydrogen composition δ²H—in comparison with the higher sector of the line (Figure 2A).

The wider combined standard deviation of the samples that fall within the lower sector is mainly due to the irregular/inhomogeneous frequency of the sampling and to the different accumulation periods of precipitation in the same station (mixing of monthly and seasonal data). In contrast, the stations on the higher sector of the line, which are not depleted by the “altitudinal effect” (lower isotope values at a higher elevation, [74,75]), show a minor combined standard deviation. This is probably due to the fact that only four stations for rainwater are located at an altitude of lower than 2500 m, and monitored by a more regular frequency of sampling (e.g., La Serena station).

An innovative method has recently been proposed for the calculation of a “generalized interval,” which combines the confidence and predictive concepts in EIV [76,77,78]. Using this approach, the obtained interval forcing a mean on 12 values per station (qx = qy = 12; i.e., as supposing a sampling of precipitation extended to 12 months) is 8.9 ± 4.2‰ up and down the fitted line, which is fairly similar to the mean combined standard deviation of δ²H (9.1 ± 4.8‰, excluding the station #9 with N_a.p. = 1, Table 2; Supplementary File S3). In contrast, the prediction intervals calculated through the jackknife method [71]—a statistical approach which simulate the resampling of the collected data taking into account their variance—and through Deming regression [71]—which necessitate a constant variance ratio (constant λ [71], calculated on the variances ratio of the 32 samples)—give a prediction interval of ±1.53‰, which is narrower than the combined standard deviation (Supplementary File S3).

Considering that the OLSR and EIV meteoric water lines of northern Chile are not different in terms of slope, the stable isotope ratio data on 115 precipitation samples available from the existing literature (Supplementary File S4; [28,30,36,56,57,58,59,60,68,79]) could help with checking the effectiveness of prediction and generalized intervals. While all these isotope data do not report the amount of precipitation and are therefore not included in the EIV regression calculation, most of them are clustered between or extended over the enriched and depleted sector of the meteoric water line, thus completing the station’s gap. The upper and lower 95% prediction intervals on the δ²H values related to the OLSR regression of this latter dataset are shifted ±13.4‰ up and down respectively on the regressed line (Supplementary File S4). This latter value is not significantly different, neither from the upper/lower δ²H combined standard deviation values of the weighted-amount dataset (±13.8‰), nor from the upper/lower extremes of the generalized interval (±13.1‰). This is particularly true if we take into account that the unweighted rainfall dataset pertains to single rain events, meaning it is quite normal that their predicted δ²H values are shifted towards the higher value of the combined standard deviation of the weighted rainfall dataset.

To verify the effectiveness of this approach, we compared the EIV regression and the above described generalized interval with the isotope composition of the following samples from northern Chile: (i) the groundwater from Pampa del Tamarugal [69] (Figure 1B); (ii) the snow samples (Supplementary File S5; [43,58,59,61,62,63,64,65,66,67,68]); and (iii) the geothermal waters samples [70].

In terms of the first groundwater sample type (i), in previous studies, the clustering of the δ¹⁸O and δ²H values of the groundwater from Pampa del Tamarugal falling on the right side of the meteoric water line has led several authors to consider this aquifer system as isotopically different from local precipitation [35,80]. This has been attributed to the following: (i) a recharge in climatic conditions different from the current ones [35]; and (ii) to the evaporation that affects the precipitation in the unsaturated zone of the recharge area [80]. A recent re-evaluation of the isotope data published up to that point shows that the isotope composition of the groundwater from that area is effectively parallel to the meteoric water line, while it was noted that the previous hypothesis on climate variation must be reconsidered [69]. However, the “uncertainty wings” around the meteoric water line have not been traced in any of the above-described cases.

In Figure 2B it is shown how the recent calculated regression line (in orange) representing the groundwater in Pampa del Tamarugal [69] falls within the generalized intervals of the EIV meteoric water line. This demonstrates that the mean composition of the groundwater from this area is not significantly different from the precipitation. This does not exclude the fact that some groundwater samples could be significantly affected by evaporation during recharge. Indeed, the right side of the generalized interval of the EIV regression is the location of the enriched precipitation samples associated with the evaporation during rainfall. Furthermore, the regression line of the groundwater samples should have a proper prediction interval, which certainly does not cover exactly that of the rainwater on that side. However, it is beyond the scope of this study to show also the prediction interval related to groundwater regression.

In terms of the snow samples, all but two of the 61 snow samples from northern Chile fall within the area between the generalized intervals determined in this study, confirming that this area could be approximated as an uncertainty ribbon around the EIV regression (Figure 2B). It should be noted that this dataset includes both samples from precipitation samplers and the Andean snow cover. The two samples out of this area are one outlier (left side of the meteoric water line, Figure 2B) and one isotopically enriched sample from the so-called “penitents” of the Parinacota volcano ([64], Figure 1B—Region I, Supplementary File S5). The penitentes are an annual phenomenon of snow sublimation related to solar radiation and arid climates; consequently, the isotope composition of these samples could be enriched, in particular, if sampled from surface layers [64,67,81]. The volume weighted mean and the related combined standard deviation calculated on the isotope composition of the snow samples from the Cerro Tapado glacier ([67], Figure 1B—Region IV, Supplementary File S5) fall exactly on the EIV-LMWL’s extension beyond the most negative value and within the generalized interval (Figure 2B). As expected, the penitentes from the same area are enriched in comparison with the mean isotope composition of the local snow due to sublimation (Figure 2B, [67]). Moreover, although their mean values fall close to EIV-LMWL, the combined sigma of the penitentes can be higher than the generalized interval (up to ±16‰; Figure 2B, Supplementary File S5).

Finally, geothermal water samples from a recent study were plotted with the EIV regression line and the generalized interval, showing more clearly the distinction between the meteoric waters and the hydrothermal water samples that can be isotopically enriched as a result of high-temperature isotopic fractionations between water–rock, mixing with andesitic water or water–gas interactions [70,74] (Figure 2B). Specifically speaking, the Torta de Tocorpuri (TT in Figure 2B) and Pampa Lirima (PL in Figure 2B) samples are not different from meteoric waters. According to the Tassi et al. [70], the limited temperature and involvement of the so-called “andesitic water” does not produce a relevant shift from the meteoric composition in these two sample groups (Figure 2B). In particular, the PL sample group lies precisely between Pampa del Tamarugal line and the EIV local meteoric water line (Figure 2B).

5. Conclusions

The regression analysis based on the error-in-variables approach (EIV) produces a meteoric water line with a slope similar to those obtained with ordinary least squares regression or other regression approaches.

This most likely occurs because the variance ratio (i.e., the square of the combined standard deviations) is high (δ²H_var/δ¹⁸O_var >> 3; [2]).

Further studies that focus on other areas and which use a greater number of amount-weighted data pairs (for example, the GNIP/WMO dataset [37]) could be useful in terms of verifying this.

The present study offers an EIV local meteoric water line that summarizes all the previously published isotope data of precipitation from northern Chile. Moreover, the use of the innovative generalized interval (locally, δ²H = ±13.1‰ up and down the fitted values on the line), which is a good compromise between the confidence and the prediction intervals, appears to be a useful tool for distinguishing significantly different data from precipitation.

The case of the groundwater system of Pampa del Tamarugal, snow samples and the local geothermal waters is emblematic. An efficient application of the generalized interval or prediction band could be extended to waters that have been subjected to evaporative fractionation (surface water or precipitation samples from desert climate areas), while it could be used to distinguish the outliers.