**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**

#### **Tiziano Boschetti 1,\*, José Cifuentes 2, Paola Iacumin <sup>1</sup> and Enricomaria Selmo <sup>1</sup>**


Received: 13 March 2019; Accepted: 11 April 2019; Published: 16 April 2019

**Abstract:** In this study, a revision of the previously published data on hydrogen (2H/ 1H) and oxygen (18O/ 16O) 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: <sup>δ</sup>2H <sup>=</sup> [(7.93 <sup>±</sup> 0.15) <sup>δ</sup>18O] + [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 <sup>δ</sup>2H = <sup>±</sup>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.

**Keywords:** precipitation; stable isotope ratios; local meteoric water line; amount-weighted mean; linear regression; confidence; prediction and generalized intervals

#### **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 (δ18O, the independent variable *x*) and hydrogen (δ2H, 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 δ18O and δ2H 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 δ18O and δ2H 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–16].

**Figure 1.** (**A**) Localization of the area of this study. (**B**) Morphostructural domains (shades of brown) and administrative regions (modified from [17]). In blue: Pampa de Tamarugal aquifer. (**C**) Morpho-tectonic units of the central Andes (modified from [18]). SBS: Santa Barbara System; SP: Sierras Pampeanas. The "Modern Forearc" coincides with the Precordillera and the Coastal Cordillera [19].

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 δ18O 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 δ18O signal of precipitation is influenced by a combination of factors ([24]). The δ18O -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 <sup>δ</sup>18O from <sup>−</sup>2 to <sup>−</sup>8%, but <sup>δ</sup>18O becomes more depleted toward the west as vapor is lifted across the Eastern Cordillera of the Andes ([26–29]). The dominant Altiplano summer rain has <sup>δ</sup>18O values from <sup>−</sup>8% to <sup>−</sup>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, [31–39]).

**Table 1.** Slope and intercept values of the previously published northern Chile meteoric water lines [31–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, <sup>δ</sup>2H <sup>=</sup> [(7.64 <sup>±</sup> 0.39) <sup>δ</sup>18O] <sup>+</sup> [9.44 <sup>±</sup> 2.32] and <sup>δ</sup>2H <sup>=</sup> [(8.04 <sup>±</sup> 0.45) <sup>δ</sup>18O] <sup>+</sup> [9.98 <sup>±</sup> 2.59] by PWLSR and RMA methods [6,7], respectively, both applied to all the available data (*N* = 40).



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–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 δ18O and δ2H data pair per station (Table 2; Supplementary File S1).


*Water* **2019** , *11*, 791

respectively)

 of the

**Table 2.**

Amount-weighted

 oxygen and hydrogen stable isotope ratios (δ18O and δ2H in permil versus the standard mean ocean water SMOW,

only the standard deviation related to spectrometric

File S2) this sample.

 measurement

(Supplementary

 File S1). The regression has been calculated both including (Table 3) and excluding

(Supplementary

Assuming xi to be the measurement of the stable isotope ratios expressed as δ18O and δ2H in % versus SMOW (standard mean ocean water)—determined on a specific precipitation amount pi (mm) collected over a specific time period—the standard deviation related to the amount weighted mean xw [4] of the meteoric stations with more than one accumulation periods of precipitation (N*a*.*p*. > 1) was calculated as follows [44]:

$$\sigma\_{\mathbf{W}} = \sqrt{\frac{\sum\_{i=1}^{\mathrm{N}\_{a,p.}} \mathbf{p}\_i (\mathbf{x}\_i - \overline{\mathbf{x}}\_{\mathbf{W}})^2}{\sum\_{\begin{subarray}{c} \mathrm{N}\_{a,p.} \\ \mathrm{N}\_{a,p.} \end{subarray}}}} \tag{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 σw, as described above, and a standard deviation related to the mass spectrometric measurements σms, which differs depending on the publication (Supplementary File S1). One of the easiest ways to calculate the combined standard deviation σcsd involves foreseeing the combination of σ<sup>w</sup> and σms following the law of propagation of uncertainty, that is, the square root of the sum-of-the-squares [45–47]:

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

The combined standard deviations, associated with the δ18O and δ2H 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 σ<sup>w</sup> 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 σ<sup>w</sup> was calculated among the years in which more than 70% of precipitation was analyzed for a given isotope composition [4]. Following this, the σcsd was calculated through using the σms declared in the analytical section of the publication dataset [32] or elsewhere [25]. With this, the σcsd on δ18O and δ2H 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–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–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–60], snow [43,58,59,61–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.

**Table 3.** Regression results obtained by different methods applied to the data of Table 2. (**a**) Crawford et al.'s non-weighted (OLSR, MA, RMA) and amount weighted (PWLSR, PWRMA, PWMA) regressions [6,7]. (**b**) Ordinary Least Squares Regression (OLSR) and error-in-variables regression (EIV) results by different codes [52,53,55,71]. (**c**) Comparison between regression's slopes obtained by EIV and PWLR, PWRMA, and PWRMA.


(**a**) regression results according to Crawford et al. 2014

Local Meteoric Water Line Freeware©—Australian Nuclear Science and Technology Organisation. OLSR: Ordinary Least Square Regression; RMA: Reduced Major Axis; MA: Major Axis; PWLSR: Precipitation Weighted Least Square Regression; PWRMA: Precipitation Weighted Reduced Major Axis. PWMA: Precipitation Weighted Major Axis. s.e.: standard error; rmSSEav: root mean sum of squared error average; Comparison between OLSR and the other regression models listed in (a): t-value: Student's *t*-test value; *p*: probability distribution by Microsoft®Office Excel (d.f. = N − 2).

(**b**) EIV regression results compared with OLSR


OLSR: here calculated by Microsoft®Office Excel—2013 (GOF parameter from Origin®Pro). EIV-a: Origin®Pro; EIV-b: SigmaPlot©; EIV-c: BFSL; EIV-d: Cantrell 2008; EIV-e: Real Statistics Using Excel© with a mean variances ratio of λ = Xvar/Yvar = 0.01913. C.I.: Confidence Intervals; r: Pearson's correlation coefficient; GOF: Goodness of Fit; rmSSE: root mean sum of squared error (this parameter is called "standard error" of regression in Microsoft®Excel's OLSR; in other regression codes, it could be easily calculated by the square root of the GOF parameter).

(**c**) EIV regression results compared with PWLSR, PWRMA, and PWMA


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 λ = δ18Ovar/δ2Hvar = 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):

$$\text{Northerm Chile EIV-LMWL: } \delta^2 \mathcal{H} = \left[ \left( 7.93 \pm 0.15 \right) \delta^{18} \mathcal{O} \right] + \left[ 12.3 \pm 2.1 \right] \tag{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, <sup>δ</sup>18O < <sup>−</sup>12% and <sup>δ</sup>2H < <sup>−</sup>80%) shows numerous monitored stations with wider standard deviations—especially in terms of hydrogen composition δ2H—in comparison with the higher sector of the line (Figure 2A).

**Figure 2.** (**A**) δ2H and δ18O amount-weighted mean of the samples (open circles; Table 2) along with the combined standard deviation on hydrogen and oxygen isotope composition (*y* and *x* bars, respectively); error-in-variables local meteoric water line (EIV-LMWL solid line) and the generalized interval (dashed lines; <sup>δ</sup>2H = <sup>±</sup>13.1% up and down the fitted values on the line). (**B**) Pampa de Tamarugal aquifer (thick-orange line, [69]), snow and penitentes (dark and light blue diamonds, respectively [41,43,59,61,63–68]) and hydrothermal groundwater (red triangles, [70]; PL: Pampa Lirima, TT: Torta de Tocorpuri) samples compared with the lines described in 'A'. Arrows explain the isotope effects in hydrothermal groundwater [70,74]. *y* and *x* bars in snow and penitentes depict the combined sigma related to the volume weighted mean (samples from Cerro Tapado [67], Supplementary File S5). Light-blue diamonds without error bars are single samples of the penitentes from Parinacota volcano ([64], Supplementary File S5).

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–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 <sup>δ</sup>2H (9.1 <sup>±</sup> 4.8%, excluding the station #9 with N*a*.*p.* <sup>=</sup> 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–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 δ2H 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 <sup>δ</sup>2H 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 δ2H 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–68]); and (iii) the geothermal waters samples [70].

In terms of the first groundwater sample type (i), in previous studies, the clustering of the δ18O and δ2H 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 (δ2Hvar/δ18Ovar >> 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, <sup>δ</sup>2H = <sup>±</sup>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.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4441/11/4/791/s1: Supplementary File S1—Table 2 raw data. Dark blue-highlighted bold values are the data shown in Table 2. Light blue-highlighted values are the results of partial calculations or additional parameters. Other colors depict parameters from the previously published tables (e.g., yellow-highlighted annual weighted mean in stations 1–9, 10–11, and 12–24) or warnings (grey-highlighted discharged data in stations 29–31 and 32; orange-highlighted measurement error in all the sheets). Supplementary File S2—Regression results as in Table 3 but excluding the sample #9 (Na.p. = 1) from Table 2 dataset. Supplementary File S3—Intervals of the meteoric water line calculated using Table 2 dataset by BivRegBLS [77] and Real Statistics Using Excel© [71]. In the former, code admits only Na.p. > 1 stations, whereas the latter use constant variance ratio (λ = δ18Ovar/δ2Hvar = 0.019). Supplementary File S4—δ18O and δ2H of single precipitation without amount data [28,30,36,56–60,68,79]; OLSR results and prediction band calculated on the isotope composition of the single precipitation dataset. Supplementary File S5—Sheet 1: δ18O and δ2H of single snow samples (collected from snow cover or by precipitation sampler; [43,58,59,61–66,68]). Sheet 2: volume weighted mean and combined standard deviation on snow and penitentes [67].

**Author Contributions:** Study design, data search, data and interpretation, writing the manuscript—original draft preparation, communicating with the journal, T.B.; data search, components for discussion on results and conclusions, writing—review and editing, J.C.; components for discussion on results and conclusions, writing—review and editing, P.I. and E.S.

**Funding:** This research received no external funding.

**Acknowledgments:** We are indebted to Bernard Francq, Université Catholique de Louvain—Belgium, who gave a great deal of assistance in the R code and for a critical revision of a first version of the manuscript. Special thanks to Javier Uribe and Jorge Gironás, Pontificia Universidad Católica de Chile, for the raw data on precipitation collected at Salar del Huasco basin. We extend acknowledgments and thanks to all the authors who published isotope data of precipitation in the studied area. Finally, we would like to thank the Editor, Polona Vreˇca, and two anonymous reviewers for their helpful comments and suggestions.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
