Reassessing Hydrological Processes That Control Stable Isotope Tracers in Groundwater of the Atacama Desert (Northern Chile)

A collection of 514 stable isotope water samples from the Atacama Desert is being reassessed geostatistically. The evaluation reveals that adjacent Andean catchments can exhibit distinct δ18O and δ2H value ranges in meteoric waters, despite similar sample altitudes of up to 4000 m above sea level (a.s.l.). It is proposed that the individual topographic features of each catchment at the western Andean Precordillera either inhibit or facilitate vapor mixing processes of easterly and westerly air masses with different isotopic compositions. This process likely causes catchment-specific isotope value ranges in precipitations (between −7‰ and −19‰ δ18O) that are being consistently reflected in the isotope values of groundwater and surface waters of these catchments. Further, due to evaporation-driven isotopic fractionation and subsurface water mixing, isotope samples of the regional Pampa del Tamarugal Aquifer plot collectively parallel to the local meteoric water line. Besides, there is no evidence for hydrothermal isotopic water-rock interactions. Overall, the observed catchment-dependent isotope characteristics allow for using δ18O and δ2H as tracers to delineate regionally distinct groundwater compartments and associated recharge areas. In this context, δ18O, δ2H and 3H data of shallow groundwater at three alluvial fans challenge the established idea of recharge from alluvial fans after flash floods.


Introduction
In the hyper-arid environment of the Atacama Desert in northern Chile, the regional Pampa del Tamarugal (PdT) Aquifer forms a strategic and vital source of groundwater, being the largest recognized aquifer in northern Chile [1,2].Its resources are used in the mining industry, for agricultural purposes or human consumption [1].Remarkable for groundwater resources of the PdT Aquifer is the wide range of δ 18 O values that span from −13‰ to −6‰ [3][4][5][6].To explain this observation different conceptual hydrogeological models were proposed for the study area, among others a deep interbasin fracture flow of groundwater from the Andes to the Atacama Desert [3,5,[7][8][9][10][11].The encountered conceptual uncertainties impede the development of a sound hydrogeological model of the PdT Aquifer to support long-term groundwater management measures [11].Recharge areas and groundwater inflow are not sufficiently understood [11].At the same time, a reliable hydrogeological model is urgently needed to support water management because current groundwater production amounts from the PdT Aquifer account for ~4000 L/s (instantaneous) and groundwater resources are being persistently overexploited as in many other regions of northern Chile [2,12,13].
Stable isotopes in water of oxygen ( 16 O and 18 O) and hydrogen ( 1 H and 2 H) are typically inert and conservative in mixing relationships and used as a tracer to investigate hydrogeological and hydrological processes including groundwater recharge and groundwater-surface-water interaction [14,15].Accordingly, the deuterium excess (D ex ) parameter of water samples can be used to detect relative evaporative enrichment effects particularly in arid regions [14].In this regard, the isotopic composition of precipitation is the initial isotopic signal of derived meteoric waters (surface, vadose zone or groundwater).The subsequent change of their isotopic composition by either evaporation, mixing or interaction with the lithosphere can be utilized to trace these processes [14].
Based on a geostatistical and chemical assessment of a dataset of 514 stable isotope measurements, which includes data of previous publications as well as unpublished data, the presented study aims at deriving new insights into hydrological processes that control the stable isotope chemistry from the arid Andes to the PdT.It is demonstrated that a consistent understanding of the regions stable isotope hydrology, in conjunction with 3 H data, allows for clarifying some of the encountered conceptual hydrogeological uncertainties.
The findings yield a new explanation for the long-lived question why stable isotope values from the PdT Aquifer cluster parallel to the local meteoric water line (LMWL) and make it possible to differentiate different recharge areas.The results have general implications for the isotope hydrology of Atacama basins and challenge the established idea of a substantial alluvial fan recharge by occasional flooding events in the open desert.

Study Area
The PdT is a nonmarine, intramassive forearc basin enclosed by the Precordillera of the central Andes to the east and the Chilean Coastal Cordillera to the west (Figure 1A,B).It is situated between altitudes of ~900-1800 m a.s.l. and forms a part of the Atacama Desert.The PdT evolved together with the Altiplano Plateau uplift from the Oligocene to present over a time span of about 30 Ma [16].The sedimentary fill of the PdT is of a lens-like shape and reaches ~1000-1700 m thickness [17].It consists of horizontally deposited nonmarine Oligocene to modern clastic sediments intercalated with extrusive volcanic deposits (e.g., ignimbrites).The aquifer of the PdT is formed by the Middle Miocene to Quaternary sediment fill with a maximum depth of ~300 m b.g.l.[11,18].The PdT Aquifer exhibits various strata with a total thickness of 25-150 m [2,18].A longitudinal cross-section of the aquifer was elaborated by Rojas et al. [11].Mean corrected groundwater residence times in the PdT Aquifer likely reach a few thousand years [3,10].
In the study area, the overall tectonic setting of the basin is characterized by N-S, NNE-SSW and NNW-SSE striking west-vergent and east-vergent (basically blind) reverse faults [17].One of the major faults is the ~65 km north-south striking Longacho Hinge, passing the localities Pica and Puquio de Nunez [17,19].
East-west discharging ephemeral stream networks flow through narrow ravines (Quebradas) downslope from parts of the Altiplano and Precordillera to the relatively flat sedimentary basin of the PdT.Catchments in the area correspond to the Quebradas (from north to south) Aroma, Tarapacá, Quipisca, Juan de Morales, Quisma, Chacarillas and Ramada (Figure 1B).
Occasional storm events in the Andean Cordillera result in flash-floods of varying intensity which can reach the alluvial fans that spread into the PdT (return period ~4 years) [20,21].Along the Andean Precordillera, several springs can be found (Figure 1B) of which the most prominent are the thermal Pica springs (at the Pica Oasis) located at ~1300 m a.s.l.(Figure 1B) [10,22].
The general groundwater flow direction of the PdT Aquifer is south-south-west (Figure 1B) [1,2,6,18].Groundwater emerges to the surface at the basins western margin where the magmatic rocks of the Coastal Cordillera constitute an impermeable barrier [16].
The general groundwater flow direction of the PdT Aquifer is south-south-west (Figure 1B) [1,2,6,18].Groundwater emerges to the surface at the basins western margin where the magmatic rocks of the Coastal Cordillera constitute an impermeable barrier [16].1B), (B) Topographical map of the study area and sample types (elevation profiles are discussed in Section 5.2).Water table contour lines indicate the regional groundwater flow regime of the PdT Aquifer and are based on Jica [2].Arrows indicate resulting groundwater flow directions.Not all springs in the region are displayed.
These emerging groundwaters can be found near-surface at salt flats in the PdT, from where they are evaporated (Pintados-and Bellavista-Salar, Figure 1B) [16,17].
A characteristic of stable isotope data from the PdT Aquifer is that δ 18 O-2 H-values of groundwater plot parallel to the local meteoric water line and δ 18 O ranges from −13‰ to −6‰ [3,5,6,23].
Overall, there are two empirical local evaporation lines (LEL) established for the study area.One LEL was derived from water samples of an unsaturated soil profile at the river outlet of the Huatacondo River in the PdT (20 km south of the Ramada basin).In this case, the LEL showed a slope of ~4 [24].A second LEL was derived based on data from the local salt flats [6].Again, a slope of ~4 was detected.Potential evaporation in the PdT during austral summer months when more than 80% of yearly precipitation falls [25], is ~250 to ~350 mm/month [26].Precipitations of more than 20 mm/a fall only above ~2500 m a.s.l.[25].Vapor masses that contribute to precipitation in the Andes can originate from east, north-east and west-north-west (most significant sources: Amazon basin and Atlantic Ocean) [23,27,28].However, the PdT itself which is situated in the Andean rain shadow did probably not experience notable amounts of precipitation in the Late Quaternary [29,30].Nevertheless, the Altiplano area went through stages of significant variations in precipitation amounts during the last 20 ka (Central Andean Pluvial Event) [31][32][33].Wetter stages in the Altiplano indirectly influenced These emerging groundwaters can be found near-surface at salt flats in the PdT, from where they are evaporated (Pintados-and Bellavista-Salar, Figure 1B) [16,17].
A characteristic of stable isotope data from the PdT Aquifer is that δ 18 O-2 H-values of groundwater plot parallel to the local meteoric water line and δ 18 O ranges from −13‰ to −6‰ [3,5,6,23].
Overall, there are two empirical local evaporation lines (LEL) established for the study area.One LEL was derived from water samples of an unsaturated soil profile at the river outlet of the Huatacondo River in the PdT (20 km south of the Ramada basin).In this case, the LEL showed a slope of ~4 [24].A second LEL was derived based on data from the local salt flats [6].Again, a slope of ~4 was detected.Potential evaporation in the PdT during austral summer months when more than 80% of yearly precipitation falls [25], is ~250 to ~350 mm/month [26].Precipitations of more than 20 mm/a fall only above ~2500 m a.s.l.[25].Vapor masses that contribute to precipitation in the Andes can originate from east, north-east and west-north-west (most significant sources: Amazon basin and Atlantic Ocean) [23,27,28].However, the PdT itself which is situated in the Andean rain shadow did probably not experience notable amounts of precipitation in the Late Quaternary [29,30].Nevertheless, the Altiplano area went through stages of significant variations in precipitation amounts during the last 20 ka (Central Andean Pluvial Event) [31][32][33].Wetter stages in the Altiplano indirectly influenced water availability in the PdT, mainly triggered by stream discharge from the Quebradas into the basin [29,30].At 21 • S latitude pluvial phases in the Altiplano can be correlated with the dating of ancient riparian vegetation in the PdT for periods between 17.6-14.2ka BP, 12.1-11.4ka BP and 2.5-2 ka BP [29,34].These periods mark likely prominent regional groundwater recharge events [29].

Groundwater Sampling
Prior to groundwater sampling, a temperature log of the well fluid was recorded with a temperature probe head.Subsequently, depth-specific groundwater sampling was carried out using a Robertson Geologging Gas Sampler.The gas sampler contains a sealed sample chamber with a moveable piston and motor-actuated valve that is led down into the well with a partial vacuum inside to retrieve uncontaminated samples of well fluids at specific depths.In order not to disturb the natural flow conditions in the well, the speed of the probe heads was regulated not to exceed 5 m/min.All other groundwater samples (from earlier studies) represent a mixture of water over the entire wells screened profile due to pumping.Wells in the PdT are typically 50-200 m b.g.l.deep with a water level typically at 20-70 m b.g.l.

Isotope Sample Analysis
For this study, stable isotope ratios of oxygen ( 18 O/ 16 O) and hydrogen ( 2 H/ 1 H) in H 2 O in water samples were measured with a PICARRO L1102-i isotope analyzer.The L1102-i is based on the WS-CRDS (wavelength-scanned cavity ring down spectroscopy) technique [35].Measurements were calibrated by the application of linear regression of the analyses of IAEA calibration material VSMOW, VSLAP and GISP.The stable isotope ratios of oxygen and hydrogen are expressed in the conventional delta notation (δ 18 O, δ 2 H) in per mil (‰) versus VSMOW [36].The reproducibility of replicate measurements is better than 0.1‰ for oxygen and 0.5‰ for hydrogen.
Stable isotopes from previously published studies were typically analyzed by mass spectrometry.Samples from studies published before the year 2000 give stable isotope ratios of oxygen and hydrogen in delta notation (δ 18 O, δ 2 H) in per mil (‰) versus SMOW.However, the reference standards SMOW and VSMOW are defined as equivalent to each other [37].
In this study, tritium data of four shallow wells at alluvial fans are presented.The data originates from an internal investigation of the mining company Compañía Minera Doña Inés de Collahuasi (CMDIC) (Santiago, Chile) [38].Tritium ( 3 H) measurements were carried out by Gas Proportional Counting at the Rosenstiel Laboratory of the University of Miami, which provides analytical errors of ±0.09 TU (tritium units).
The deuterium excess (D ex ) parameter of a water sample is calculated according to the following equation [14]:

Isotope Data
Earlier published studies include papers and reports about the isotopic content in precipitation, river water, and spring-and groundwater (Table 1, Figure 1B).The data is available attached to this article (Supplementary Materials File S1).
A geostatistical evaluation was conducted on the isotope data from spring, river and groundwater, together with unpublished data from an internal study of the mining company Compañía Minera Doña Inés de Collahuasi [38] and data collected in this study.The report by Salazar et al. [4] includes all isotope data collected by the Chilean water directorate (DGA) that were published by Magaritz et al. [5], Aravena and Suzuki [39] and Aravena [6].The precision of the GPS-coordinates of the data is restricted by the accuracy of the GPS-instruments of each report or study but was found to be of good to sufficient accurateness for a regional scale assessment.Samples of surface water from the salt flats Salar de Pintados and Bellavista and Salar del Huasco-strongly affected by evaporation-were excluded from the dataset because they would distort the assessment of regional groundwater conditions.In the PdT dataset, 170 samples are concentrated around the Pica area, where dozens to hundreds of wells can be found and several sample campaigns were carried out.The remaining 196 samples are distributed across the PdT over an area of ~9000 km 2 , including the slope of the Andean Precordillera.

Geostatistical Interpolation
The δ 2 H and δ 18 O data was processed by the geostatistical kriging method by the ArcGIS geostatistical analyst [38].Kriging allows for deriving spatial interpolation models of georeferenced point data and helps thereby to identify spatial trends in the data distribution [42].To find the most accurate model for the investigated cases an iterative, semi-automated procedure was followed, including data trend analysis, trend removal, semi-variogram modeling, automated model optimization and anisotropy correction.Models that are presented in this study represent those models that performed best when tested by cross-validation, with the lowest root mean square error (RMSE) and average standard deviation (ASD).Both chosen models (δ 2 H and δ 18 O) apply the nugget effect which allows kriging to deviate from exact measured values for yielding a lower RMSE and ASD globally [42].Allowing the nugget effect was necessary primarily because for some sample locations several measurements were available that deviated slightly from each other.This deviation was either due to sampling at different dates, different sample types (e.g., river and groundwater) or, for example, due to samples from areas that are locally affected by reinfiltrated irrigation water, which is influenced by additional evaporation (Section 4.1).The chosen models are therefore optimized for a regional trend assessment.However, all generated candidate models demonstrated similar regional trends.Model parameters of the chosen kriging models can be found in Supplementary Materials File S2.The final models were masked by the extension of the Tarapacá, Quipisca, Juan de Morales, Quisma, Chacarillas, Salar del Huasco and Laguna Lagunillas basin for enhanced visualization.
Trendlines of stable isotope data were derived by linear regression.Statistical visualization of data was done by Tukey box plots [43].

Mathematical Calculation of the Intersection of the Local Meteoric Water Line and a Given Local Evaporation Line
To calculate the point of intersection of a specific local evaporation line-with a known slope of approximately four (Section 2) and containing a sample point S-and the local meteoric water line, the following equations [23] are needed: LMWL: δ 2 H p = 7.8 × δ 18 O p + 9 (2) δ 2 H p and δ 18 O p are the isotopic values of regional precipitations, δ 2 H s and δ 18 O s are the isotopic values of a given meteoric water sample S and n s is the respective y-intercept of the specific local evaporation line.For calculating the intersection point (δ 2 H i and δ 18 O i ) where the deuterium value of the specific local evaporation line is equal to the coordinates δ 2 H p and δ 18 O p , the following equation is finally entertained: The resulting value δ 18 O i is the initial δ 18 O value in precipitation from which the isotopic composition of sample S evolved by evaporation-driven fractionation.
Equation (4) will be fed with data of the kriging models for δ 18 O and δ 2 H to yield a spatial distribution of the initial isotopic composition of meteoric waters in the study area.

Geostatistical Assessment of Stable Isotope Data from the Andean Altiplano to the PdT Aquifer
The kriging method is applied to assess spatial trends in the georeferenced δ 2 H and δ 18 O datasets that include for both parameters 435 single measurements (Figure 2A,B, Section 3.4) (26 samples from the Salar del Huasco basin were not included due to a lack of coordinates, Supplementary Materials File S1).The presented interpolation models are those that exhibit the lowest RMSE and ASD.The derived δ 2 H model exhibits an RMSE of 6.82‰ and an ASD of 5.92‰.The δ 18 O model exhibits an RMSE of 0.89‰ and an ASD of 0.80‰.Reports on the model configuration and respective cross-validation plots can be found in Supplementary Materials File S2.
The interpolation yields a clear spatial correlation that is almost identical for both parameters.Distinct meteoric water compartments (CP) within the study area are being visualized that exhibit different isotope value characteristics and correspond to the topographic watersheds in the region (Figure 2A,B).The statistical attributes of samples within the different compartments are summarized in Figure 3A,B.Recently it was demonstrated that the corresponding groundwater is recharged at the Altos de Pica recharge area (Figure 2B) [10].A slight increase of δ 2 H and δ 18 O values can be found in some samples of the local and shallow Pica aquifer, which is influenced by reinfiltrated irrigation and spring water affected by evaporation, which causes an additional isotopic fractionation [3,22,44].
Data of CP 5 (24 samples) shows a first and third quartile of δ 2 H values at −90‰ and −73‰.For δ 18 O respective values are −10.6‰and −8.1‰.It lies in the drainage area of the Quebradas Chacarillas and Ramada.Some groundwater samples of this compartment and catchment, which were taken in the Andean Altiplano, reach δ 2 H values of up to −95‰.Also here an isotopic value increase with decreasing altitude can be observed.
CP 6 corresponds to the Altiplano basins Salar del Huasco and Laguna Lagunillas (91 samples of the Salar del Huasco basin and three from the Laguna Lagunillas basin).The first and third quartile of δ 2 H values plot at −102‰ and −96‰ respectively.δ 18 O values range dominantly between −13.3‰ and −12.1‰.
The catchment-dependent isotopic characteristics are also reflected in the respective D ex values as illustrated by Tukey box plots (Figure 3C).Extreme values outside the lower and upper bound of the box plots can be considered outliers [43] (green boxes cover the data range that plots within the first and third quartiles).While CP 1 and CP 5 exhibit relatively low D ex values in majority between 1‰ and −7‰, CP 3 exhibits highest D ex values with 2-4‰.The first and third quartile of D ex data of all sample groups plot below the first and third quartile of D ex values in regional precipitations (Figure 3C).
influenced by reinfiltrated irrigation and spring water affected by evaporation, which causes an additional isotopic fractionation [3,22,44].
Data of CP 5 (24 samples) shows a first and third quartile of δ 2 H values at −90‰ and −73‰.For δ 18 O respective values are −10.6‰and −8.1‰.It lies in the drainage area of the Quebradas Chacarillas and Ramada.Some groundwater samples of this compartment and catchment, which were taken in the Andean Altiplano, reach δ 2 H values of up to −95‰.Also here an isotopic value increase with decreasing altitude can be observed.The catchment-dependent isotopic characteristics are also reflected in the respective Dex values as illustrated by Tukey box plots (Figure 3C).Extreme values outside the lower and upper bound of the box plots can be considered outliers [43] (green boxes cover the data range that plots within the first and third quartiles).While CP 1 and CP 5 exhibit relatively low Dex values in majority between 1‰ and −7‰, CP 3 exhibits highest Dex values with 2-4‰.The first and third quartile of Dex data of all sample groups plot below the first and third quartile of Dex values in regional precipitations (Figure 3C).

δ 2 H-δ 18 O-Charts of All Stable Isotope Data
Charting all PdT samples together in one δ 2 H-δ 18 O-plot confirms results of earlier studies that found that stable isotope data from the PdT plots parallel to the LMWL (Figure 4A) [3,5,6].The resulting trend line of all PdT samples shows a slope of 7.4 (R 2 = 0.89) which is almost equal to the slope of the established LMWL of 7.8 [23].
In comparison to that, Figure 4B summarizes a collection of stable isotope data (91 samples) from the Salar del Huasco basin (spring, river and groundwater).The Salar del Huasco basin abuts on watersheds draining into the PdT (Figure 1B).Contrary to the trend line of samples of the PdT, the trend line of data from the Salar del Huasco basin exhibits a slope of 3.9 (R 2 = 0.61), which is in accordance with empirically established LELs for the study area (Section 2).

δ 2 H-δ 18 O-Charts of All Stable Isotope Data
Charting all PdT samples together in one δ 2 H-δ 18 O-plot confirms results of earlier studies that found that stable isotope data from the PdT plots parallel to the LMWL (Figure 4A) [3,5,6].The resulting trend line of all PdT samples shows a slope of 7.4 (R 2 = 0.89) which is almost equal to the slope of the established LMWL of 7.8 [23].
In comparison to that, Figure 4B summarizes a collection of stable isotope data (91 samples) from the Salar del Huasco basin (spring, river and groundwater).The Salar del Huasco basin abuts on watersheds draining into the PdT (Figure 1B).Contrary to the trend line of samples of the PdT, the trend line of data from the Salar del Huasco basin exhibits a slope of 3.9 (R 2 = 0.61), which is in accordance with empirically established LELs for the study area (Section 2).

δ 18 O against Sample Altitude for Samples of Compartments Two, Three and Five
Compartments two, three and five that correspond to the catchments Tarapacá/Quipisca, Juan de Morales and Chacarillas (Section 4.1), provide enough Precordilleran slope samples to assess the relationship between δ 18 O values and the sample altitude (Figure 5A-C).As an apparent tendency, an increase of the δ 18 O values with decreasing elevation can be observed for CP 5 and CP 2. CP 3 demonstrates no apparent increasing trend.However, available electrical conductivity data indicates generally an increase in salinity along the flow path with decreasing altitude.Isotope samples from the PdT Aquifer downstream of the respective compartments cluster typically around the upper end of the δ 18 O values for each catchment.Nevertheless, the isotopic composition of groundwater below the alluvial fan of the Juan de Morales catchment does not correlate with the isotopic composition of its meteoric water upstream (Figures 2A and 5C).Also, the isotopic values of groundwater below the alluvial fan of the Chacarillas river cannot be correlated with isotopically enriched flooding water at its alluvial fan but correspond to the isotopic composition of CP 4 and CP 5 upstream (further discussed in Section 5.3).Floodwater was sampled on the day of flooding (data from Aravena et al. [24]).

δ 18 O against Sample Altitude for Samples of Compartments Two, Three and Five
Compartments two, three and five that correspond to the catchments Tarapacá/Quipisca, Juan de Morales and Chacarillas (Section 4.1), provide enough Precordilleran slope samples to assess the relationship between δ 18 O values and the sample altitude (Figure 5A-C).As an apparent tendency, an increase of the δ 18 O values with decreasing elevation can be observed for CP 5 and CP 2. CP 3 demonstrates no apparent increasing trend.However, available electrical conductivity data indicates generally an increase in salinity along the flow path with decreasing altitude.Isotope samples from the PdT Aquifer downstream of the respective compartments cluster typically around the upper end of the δ 18 O values for each catchment.Nevertheless, the isotopic composition of groundwater below the alluvial fan of the Juan de Morales catchment does not correlate with the isotopic composition of its meteoric water upstream (Figures 2A and 5C).Also, the isotopic values of groundwater below the alluvial fan of the Chacarillas river cannot be correlated with isotopically enriched flooding water at its alluvial fan but correspond to the isotopic composition of CP 4 and CP 5 upstream (further discussed in Section 5.3).Floodwater was sampled on the day of flooding (data from Aravena et al. [24]).

Time Series of Stable Isotope Samples from Spring Sites between 1967 and 2014
The selected spring sites are the only ones in the study area that were sampled several times during the last 50 years.Figure 6

Depth-Specific Stable Isotope Samples from the PdT Aquifer
Depth-specific stable isotope samples were taken at wells J5, LC2, J8 and JF during October 2015 (location displayed in Figure 2B, data presented in Table 2).Before sampling, temperature logs were recorded at the respective wells (Figure 7).While wells J5, J8 and JF are cased and exhibit screen sections [2] well LC2 is an open borehole.The measurements at wells J5, JF and LC2 indicate that at the same locations stable isotope compositions of groundwater stayed almost constant over depth, with a tendency to slightly increase in salinity.Only at well J8 a marked difference between two sample depths can be identified.While the sample J8.2, taken at 170 m b.g.l., shows a δ 18 O value of −10.9‰ with a conductivity of 875 µS/cm, the sample J8.1 (60 m b.g.l.) demonstrates −9.3‰ δ 18 O and 997 µS/cm.Also, the Dex value differs by 4.6‰.

Time Series of Stable Isotope Samples from Spring Sites between 1967 and 2014
The selected spring sites are the only ones in the study area that were sampled several times during the last 50 years.Figure 6

Depth-Specific Stable Isotope Samples from the PdT Aquifer
Depth-specific stable isotope samples were taken at wells J5, LC2, J8 and JF during October 2015 (location displayed in Figure 2B, data presented in Table 2).Before sampling, temperature logs were recorded at the respective wells (Figure 7).While wells J5, J8 and JF are cased and exhibit screen sections [2] well LC2 is an open borehole.The measurements at wells J5, JF and LC2 indicate that at the same locations stable isotope compositions of groundwater stayed almost constant over depth, with a tendency to slightly increase in salinity.Only at well J8 a marked difference between two sample depths can be identified.While the sample J8.2, taken at 170 m b.g.l., shows a δ 18 O value of −10.9‰ with a conductivity of 875 µS/cm, the sample J8.1 (60 m b.g.l.) demonstrates −9.3‰ δ 18 O and 997 µS/cm.Also, the D ex value differs by 4.6‰.

Time Series of Stable Isotope Samples from Spring Sites between 1967 and 2014
The selected spring sites are the only ones in the study area that were sampled several times during the last 50 years.Figure 6

Depth-Specific Stable Isotope Samples from the PdT Aquifer
Depth-specific stable isotope samples were taken at wells J5, LC2, J8 and JF during October 2015 (location displayed in Figure 2B, data presented in Table 2).Before sampling, temperature logs were recorded at the respective wells (Figure 7).While wells J5, J8 and JF are cased and exhibit screen sections [2] well LC2 is an open borehole.The measurements at wells J5, JF and LC2 indicate that at the same locations stable isotope compositions of groundwater stayed almost constant over depth, with a tendency to slightly increase in salinity.Only at well J8 a marked difference between two sample depths can be identified.While the sample J8.2, taken at 170 m b.g.l., shows a δ 18 O value of −10.9‰ with a conductivity of 875 µS/cm, the sample J8.1 (60 m b.g.l.) demonstrates −9.3‰ δ 18 O and 997 µS/cm.Also, the Dex value differs by 4.6‰.Besides, the temperature log indicates that in the case of wells J8 and LC2 a temperature gradient increase occurs at round about 860 m a.s.l.At well J5 a continuous temperature increase is observed.At well JF the gradient appears to decrease slightly from ~900 m a.s.l.on.Noteworthy is the shift of the temperature profiles depending on the earlier delineated aquifer compartments (Section 4.1).

3 H Samples from Shallow Alluvial Fan Groundwater in the PdT
3 H values of groundwater of four shallow wells situated at alluvial fans (af) show that no significant amounts of 3 H can be detected (Table 3 and Figure 2).A sample from the 54 m deep well C1 (water level at ~20 m b.g.l.) at the western part of the Chacarillas fan in CP 4 shows no apparent tritium activity.Also, two samples of the Juan de Morales fan were taken that is situated in CP 2. The samples show tritium values ≤0.09TU.A sample from the Quipisca alluvial fan (located in CP 2, Figure 2A) demonstrates likewise values of ≤0.09 TU.Besides, the temperature log indicates that in the case of wells J8 and LC2 a temperature gradient increase occurs at round about 860 m a.s.l.At well J5 a continuous temperature increase is observed.At well JF the gradient appears to decrease slightly from ~900 m a.s.l.on.Noteworthy is the shift of the temperature profiles depending on the earlier delineated aquifer compartments (Section 4.1).

3 H Samples from Shallow Alluvial Fan Groundwater in the PdT
3 H values of groundwater of four shallow wells situated at alluvial fans (af) show that no significant amounts of 3 H can be detected (Table 3 and Figure 2).A sample from the 54 m deep well C1 (water level at ~20 m b.g.l.) at the western part of the Chacarillas fan in CP 4 shows no apparent tritium activity.Also, two samples of the Juan de Morales fan were taken that is situated in CP 2. The samples show tritium values ≤0.09TU.A sample from the Quipisca alluvial fan (located in CP 2, Figure 2A) demonstrates likewise values of ≤0.09 TU.

Hydrological Processes Controlling Stable Isotope Tracers in the Salar del Huasco Basin (Altiplano)
A significant hydrological process that alters stable isotope values of groundwater in arid regions is isotopic fractionation by evaporation that affects meteoric water before infiltrating into the underground [14,15].It causes the δ 18 O-δ 2 H composition of groundwater to deviate from the trend of the LMWL-following a LEL and leads typically to an increase of the δ 18 O value and hence a decrease of the D ex value.Also, D ex values in groundwater that are significantly lowered compared to D ex values of regional precipitations indicate an isotopic fractionation by evaporation prior to recharge [14].
The effect of isotopic fractionation by evaporation can be observed well for data from the Salar del Huasco basin (Figure 4B).The trend line of respective samples exhibits a slope of 3.9 and is thereby in accordance with the slope of established empirical LEL in the region showing a value of ~4 (Section 2).Hence, the general distribution of stable isotope data of spring, river and groundwater from the Salar del Huasco basin reflects the effect of an isotopic enrichment induced by evaporation of meteoric water prior to infiltration and is satisfyingly explained by this process.The dominance of isotopic fractionation by evaporation is confirmed by the lowered D ex values of the dataset compared to values of regional (winter and summer) precipitations (Figure 3C, CP 6).However, Figure 4B also demonstrates a certain fuzziness of the data clustering in a corridor above and below the calculated local evaporation trend line, which is quantified by a moderately good regression coefficient (R 2 ) of 0.61.It is proposed that the cause of this fuzziness is the natural variation of the initial δ 18 O-δ 2 H composition of precipitations that fall in the Salar del Huasco basin.These variations can be influenced by amount, altitude and continental effects, that are known to occur in the Andean Altiplano [8, 23,27].Accordingly, the first and third quartile of δ 18 O values of sixteen weighted rainwater samples of the Salar del Huasco basin plot between −17‰ and −13‰ (based on data of Fritz et al. [3], Aravena et al. [23] and Uribe et al. [8]).In accord with that, the intersections of the estimated upper and lower bound of the identified LEL with the LMWL mark a δ 18 O value range between −17‰ and −14‰ (Figure 4B).That is because waters with isotope values along the upper and lower bounds must have evolved from rainwater departing from the LMWL along an evaporation line with a slope of ~4.
An additional effect that probably influences this isotope data are groundwater mixing processes occurring between the highest and the lowest isotope values in the respective aquifer (conservative mixing within red circle in Figure 4B).However, the most extreme isotope values that were used to delineate the predominant isotope value range of precipitations can be considered to be (almost) unaffected by mixing processes.That is because the datasets suggest that margin values represent the extreme ends of the investigated group and hence they cannot be the result of a substantial mixing.Therefore, mixing processes within the respective meteoric water compartment are negligible for estimating the predominant value range of feeding precipitations.
Another factor that must be taken into account when assessing stable isotope data in the Andean area are interactions with the lithosphere under hydrothermal conditions [45,46].Although previous isotopic studies in the Andean Altiplano suggest that isotopic exchange processes can occur in arid basins associated with hot springs [45,46], isotopic interactions with the lithosphere cannot be observed for data from the Salar del Huasco basin.Such an effect would cause a limited set of samples -particularly of thermally influenced springs-to deviate from the evaporative trend observed for the basins meteoric water.This is not the case.Overall, the Salar del Huasco basin exhibits only slightly thermally influenced springs (~22-25 • C) but no hot springs [8,10].

Hydrothermal Water-Rock Interactions
Jayne et al. [9] discussed the possibility that isotopic anomalies encountered for waters close to the Juan de Morales basin (CP 3), might be caused by hydrothermal water-rock interactions.In fact, the Juan de Morales basin exhibits some hot springs (water temperature ~50  2B and 5C).Despite the absence of a hydrothermal imprint on these cold springs, the respective waters show likewise relatively high δ 18 O values that are in perfect accord with values measured at nearby hot springs.Hence, it is not likely that hydrothermal isotope exchange processes affect the isotopic composition of groundwater in CP 3 despite a possible circulation along deep faults into fractured segments.Also, thermal waters at Pica which circulate along disruption zones within the Andean slope to depths of ~900 m b.g.l.[10], show no anomalous deviation in their isotopic composition from surrounding waters in CP 4 (Figure 2A,B and Figure 6).Overall there is, therefore, no evidence for the assumption that hydrothermal water-rock interactions affect the isotopic composition of groundwater west of the Andean Precordillera.Hence, there must be other effects involved that cause catchment-dependent isotopic value ranges.

Isotopic Fractionation by Evaporation of Meteoric Waters and Initial Isotopic Composition of Feeding Rainwater
It is evident that spring-and groundwater from the PdT and Andean Precordillera must be strongly affected by an isotopic fractionation caused by evaporation.Apart from the high potential evaporation in the region (Section 2), this is demonstrated by two facts.Firstly, D ex values of meteoric waters west of the Andean Precordillera are consistently lower than in regional precipitations (Figure 3C).Secondly, there is evidence of a continuous δ 18 O enrichment with decreasing altitude along at least two catchments in the region (Figure 7A,B).These two arguments allow for the conclusion that discharging spring and river water from Andean slope catchments is affected by isotopic fractionation caused by an ongoing evaporation on its way downstream and should evolve dominantly along an evaporation line as in the case of the Salar del Huasco basin (Section 5.1).
When reviewing the different datasets of the spatial compartments of Section 4.1 (Figures 2A  and 3), the isotopic data of CP 1-5 plot indeed in clusters similar to the data of the Salar del Huasco basin (Figure 3A,B).Analogous to the procedure applied to data from the Salar del Huasco basin-by estimating an upper and lower bound for evaporation-driven fractionation-the predominant isotope value range of feeding precipitations of CP 1-5 can be estimated.Instead of following a graphical procedure as in the case of the Salar del Huasco basin (Figure 3B), the assessment can be carried out mathematically based on the data provided by the kriging models of δ 18 O and δ 2 H (Figure 2A,B, Section 3.5).
The assessment proposes that rainwater falling in the Juan de Morales basin exhibits predominantly δ 18 O values between −11 and −7‰ (Figure 8A).On the contrary precipitations in CP 6 and CP 5 would exhibit dominantly δ 18   This is emphasized by Figure 8A that displays the distribution of the theoretical initial δ 18 O value of meteoric waters prior to evaporation-driven fractionation, when neglecting eventual groundwater mixing processes (mixing occurs primarily in the PdT-Aquifer) (Section 3.5).Anomalies in the initial δ 18 O values of meteoric waters are associated with high topographic elevations (~above 5000 m a.s.l.) that form barriers against easterly air masses with depleted δ 18 O values (Figure 8A).Further, Figure 8B depicts the contrast in δ 18 O values in precipitations east and west of the Andean mountain range, which is caused by different vapor sources situated either east-northeast or west-northwest of the Andean mountain range.Similar observations of a rapid change in vapor isotope values when comparing windward and leeward sides of mountain ranges were also made, for example, in the Canadian Rocky Mountains and the Himalaya Mountains [14,47,48].
In this understanding, the Juan de Morales basin exhibits a mountainous barrier against the windward side of easterly trade winds, with maximum elevations above 5000 m a.s.l. and so does the western limit of the Tarapacá basin.This condition probably leads to a rainout at their eastern margins, due to the adiabatically cooling of rising easterly air masses.Eventually crossing moistures could mix with predominantly north-western air masses and yield relatively elevated Finally, it is crucial to understand why rainfall above the Quisma basin shows partially relatively depleted δ 18 O values.Vapors that arrive from east-northeast on the mountainous barrier of the Juan de Morels basin will be depleted on its leeward side due to altitude and amount effects (rainout).Altitude effects are observed to occur in the Salar del Huasco basin with a δ 18 O depletion of approximately −0.64‰/100 m [8].These depleted vapor masses probably arrive occasionally with northeastern winds at the adjacent Altos de Pica mountain ridge of the Quisma basin and precipitate (Figures 2B and 8A,B).
While altitude, amount and continental effects affecting regional precipitations in the study area were described before [23,27,49], the importance of local topographic features on the inhibition or facilitation of vapor mixing processes was disregarded so far when assessing stable isotopes of groundwater in the Atacama Desert.Overall, the mentioned processes (catchment-dependent isotopic composition of rainwater, evaporation and eventual groundwater mixing) would cause the samples to shift from the LMWL collectively but at the same time deviate from the respective LEL.
Consequently, the reason for the different δ 18 O value ranges in precipitations for different catchments will be discussed.

Catchment-dependent δ 18 O Value Ranges of Precipitations: Topographic Controls on Air Mass and Vapor Mixing
Previous authors speculated about the reasons why collectively examined δ 18 O-δ 2 H samples from the PdT Aquifer show such a wide value range and plot parallel to the LMWL (Figure 4A) [5,6].Reasons proposed for this observation were different types of precipitation (snow and rain) or a change in climatic conditions (paleo-rain) [5,6,23], as well as assumptions about recharge areas at different altitudes (altitude effect) [6].There are two arguments to disregard the hitherto proposed reasons.
The first argument is that data from the Salar del Huasco basin of the nearby Andean Altiplano consistently plot along an expected LEL (Figure 4B).Apparently, this water must have been recharged under the same climatic conditions as the water from the PdT Aquifer because it is known that pluvial phases in the Altiplano foster recharge in the PdT through discharging rivers (Section 2).To assume that different types of precipitation or a variation in climatic conditions would have changed or affected the isotopic composition of rainfall in northern Chile, to be similar to waters from the PdT, is therefore not justified.
The second argument is that meanwhile it is known that altitude effects are limited to the Andean Altiplano.Such altitude effects develop with the elevation increase along the eastern Andean flank and are consequently linked to abundant summer rainfalls with mainly easterly and north-easterly vapor sources (originating from the Amazon basin and the Atlantic Ocean) [23,27].Typical altitude effects do not apply to vapor masses rising along the western flank of the Andes with dominantly west-north-western vapor sources [23,27].Hence, different recharge altitudes cannot explain the spatially strongly varying isotope values in the PdT that exist up to sample altitudes above 4000 m a.s.l.
It is argued that the main reason for the observed variation of δ 18 O and δ 2 H values in the different compartments of the PdT Aquifer (Figure 2A) are predominant long-term stable isotope value ranges in precipitations that aliment the relevant recharge areas (as discussed in Sections 5.1 and 5.2.2).These δ 18 O-δ 2 H value ranges are likely controlled by topographic characteristics of the catchments in the transition zone from the Altiplano to the PdT that can show elevation differences of up to 1200 m (Figure 1B).This is emphasized by Figure 8A that displays the distribution of the theoretical initial δ 18 O value of meteoric waters prior to evaporation-driven fractionation, when neglecting eventual groundwater mixing processes (mixing occurs primarily in the PdT-Aquifer) (Section 3.5).Anomalies in the initial δ 18 O values of meteoric waters are associated with high topographic elevations (~above 5000 m a.s.l.) that form barriers against easterly air masses with depleted δ 18 O values (Figure 8A).Further, Figure 8B depicts the contrast in δ 18 O values in precipitations east and west of the Andean mountain range, which is caused by different vapor sources situated either east-northeast or west-northwest of the Andean mountain range.Similar observations of a rapid change in vapor isotope values when comparing windward and leeward sides of mountain ranges were also made, for example, in the Canadian Rocky Mountains and the Himalaya Mountains [14,47,48].
In this understanding, the Juan de Morales basin exhibits a mountainous barrier against the windward side of easterly trade winds, with maximum elevations above 5000 m a.s.l. and so does the western limit of the Tarapacá basin.This condition probably leads to a rainout at their eastern margins, due to the adiabatically cooling of rising easterly air masses.Eventually crossing moistures could mix with predominantly north-western air masses and yield relatively elevated δ 18  Finally, it is crucial to understand why rainfall above the Quisma basin shows partially relatively depleted δ 18 O values.Vapors that arrive from east-northeast on the mountainous barrier of the Juan de Morels basin will be depleted on its leeward side due to altitude and amount effects (rainout).Altitude effects are observed to occur in the Salar del Huasco basin with a δ 18 O depletion of approximately −0.64‰/100 m [8].These depleted vapor masses probably arrive occasionally with northeastern winds at the adjacent Altos de Pica mountain ridge of the Quisma basin and precipitate (Figures 2B  and 8A,B).
While altitude, amount and continental effects affecting regional precipitations in the study area were described before [23,27,49], the importance of local topographic features on the inhibition or facilitation of vapor mixing processes was disregarded so far when assessing stable isotopes of groundwater in the Atacama Desert.
Overall, it is this phenomenon of locally varying isotope value ranges in precipitations that allows for using stable isotopes in the PdT as a useful tracer for groundwater flow mapping.The influence of Andean topographic features on the movement of air masses and vapor mixing are very likely of relevance for similar basins of the Atacama Desert where basin-dependent pronounced topographic differences occur.

Kriging Models, Regional Groundwater Flow Regime and Identified Recharge Areas
The correlation structure that is being mapped by the kriging models (Figure 2A,B) is that of the spatial distribution of the isotopic composition of meteoric waters in the region.The respective RMSE and ASD of both models (δ 2 H model: RMSE = 6.82‰,ASD = 5.92‰ and δ 18 O model: RMSE = 0.89‰, ASD = 0.80‰) are relatively low and reasonable when considering that the dataset consists of different kinds of meteoric water (surface and groundwater) and includes samples of campaigns spanning over a period of ~50 years.Particularly spring sites were sampled several times during that period.δ 18 O measurements at these springs show an ASD of 0.54‰ (Figure 5), which is quite close to the ASD of the respective kriging model (ASD of respective δ 2 H values is 4.2‰).However, the low ASD of isotopic measurements at springs also demonstrates that the isotopic composition of groundwater at different localities is not exposed to prominent (seasonal) variations.That is reasonable when taking into account that groundwater in the Atacama Desert exhibits typically corrected residence times of a few thousand years (e.g., Pica springs and Puquio de Nunez spring [3,10]).Long residence times damp signals of slightly varying seasonal variations due to mixing within the respective reservoir.Hence, extreme isotopic values-as associated with depleted rains during storm events-are smoothed out in the given spring data (as far as observed).Applying this insight to other groundwater measurements in the region implies that the observed isotopic conditions-as visualized by the kriging models-reflect long-term averaged conditions due to subsurface mixing within each compartment.However, based on the given dataset alone it is impossible to quantify or further distinguish these mixing processes.
The interpolations of δ 2 H and δ 18 O data confirm each other in their spatial extent, which allows for identifying six meteoric water compartments that are apparently physically related to different catchments in the study area (Section 4.1, Figure 2A,B).Due to the catchment-dependent isotope value ranges (Section 5.2) the spatial trends in δ 2 H and δ 18 O data can be used to trace the regional groundwater flow regime and indirectly derive conclusions about the recharge areas of respective waters.The fact that the geostatistically modeled spatial molding of CP 1-5 is in high accordance with the regional groundwater flow direction-as indicated by regional water table contour lines-further substantiates the validity of the identified spatial trends (Figures 1B and 2A,B).Besides, a preliminary comparison of electrical conductivity measurements for CP 2, 3 and 5 (Figure 5), also confirms the delineated groundwater flow regime by a continuous increase of the salinity along the flow path.Also, an earlier report on major ions of the central PdT Aquifer distinguished very similar aquifer compartments based on chemical groundwater types [50].
Hence, groundwater in CP 1 likely originates from the Aroma basin, groundwater in CP 2 is recharged at the Tarapacá and Quipisca basins, and groundwater from CP 3 infiltrates in the Juan de Morales basin.Groundwater in CP 4 is recharged in the Quisma basin (in particular at the Altos de Pica area [10]) and groundwater associated with CP 5 percolates in the Chacarillas catchment.Finally, close to the western margin of the PdT where near-surface water prevails at salt flats, further evaporation processes probably cause waters to enrich isotopically increasingly as reflected by the respective kriging models (Figure 2A,B).These enriched brines might affect groundwater close to salt flats isotopically by either density-driven recirculations or modern groundwater pumping from production wells (change of natural flow conditions).
A first preliminary observation of the distinct patterns in the distribution of δ 2 H values in the PdT (based on less than 25 samples and without any quantitative method applied) was made by Fritz et al. [3].Nevertheless, the results presented in the current study yield for the first-time a statistical assessment of the isotopic conditions based on an extensive stable isotope dataset.
Depth-specific samples indicate that the geostatistical trend of isotope values in the region is consistent over varying depths of the aquifer (up to ~250 m b.g.l.) despite the existence of different aquifer storeys.This observation is reasonable when considering as one of the primary recharge mechanisms an inflow from slope infiltrations of streambeds together with water mixing along the river course and subsurface flow path (check graphical abstract).The temperature profiles of wells J8 and LC2 indicate that very likely different aquifer storeys were sampled because of a change in the temperature gradient (Figure 6).However, due to the constantly rising temperature profiles at wells J5 and JF, here vertical flows could occur and lead to water mixing within the wells.Nevertheless, the deviation of a deep and shallow sample from J8 could indicate a mixing of water from CP 4 and 5 (Figure 7).Such a mixing, in turn, would be in accordance with the spatial distribution of observed regional groundwater compartments.Based on a binary isotopic mixing model-incorporating water from the Puquio de Nunez area and groundwater from the fan apex of the Chacarillas fan-the shallow sample at well J8 consists of a water mix of 50:50 from both reservoirs.The deeper sample at J8 is almost identical to groundwater from the Puquio de Nunez area.
The spatial extension of the identified regional isotopic CP 4 also demonstrates that the groundwater which is recharged at Altos de Pica through rainfall is passing the Longacho Hinge (conservatively treated as an impermeable barrier; [11,18]) and finally discharging into the PdT Aquifer (Figure 2A).A recent geophysical survey (gravity and TEM measurements) supports that the bedrock of the Longacho Hinge dips section-wise deeper into the underground, facilitating pathways for a shallow inflow from the Andean slope through sedimentary units into the PdT Aquifer [44].

Recharge Mechanisms
The isotopic assessment allows for some conclusions regarding recharge mechanisms of the PdT Aquifer.An indication that can be derived from the spatial isotope distribution is that alluvial fan recharge in response to flooding is probably negligible for the Juan de Morales and the Chacarillas alluvial fans.It was proven that water, flooding the Chacarillas alluvial fan after storm events, is rapidly enriching in δ 18 O by ~1.2‰/day due to evaporation, and is relatively enriched with values above −4‰ δ 18 O (Figure 5A) [24].Furthermore, it is known that evaporation affects percolating waters up to 2 m b.g.l. in the Atacama Desert [51].
Infiltrating floodwater caused by storms should, therefore, exhibit significantly higher δ 18 O and δ 2 H values when compared to water recharged along the river course and fan apex.As can be seen from Figure 5A, floodwater at the Chacarillas fan recently sampled after flooding, can isotopically not be correlated with shallow groundwater at the well J8 (water level at well J8 is ~39 m b.g.l.).On the contrary, as visualized in Figures 2A and 5A,C the isotopic composition of water below the alluvial fans of the Chacarillas (CP 5) and Juan de Morales (CP 3) catchments demonstrates lower δ 18 O and δ 2 H values than exhibits the water of the relevant basins upstream.This finding is very likely due to groundwater flowing in from adjacent catchments as also indicated by the groundwater level contour lines (Figure 1B).The water of CP 2 can be found below the Juan de Morales alluvial fan and water of CP 4 below the Chacarillas alluvial fan (Figure 2A).The inflow from other basins is confirmed by a preliminary comparison of respective electrical conductivities of groundwater (Figure 5).Water below the Juan de Morales fan, for example, corresponds with values of 1500-3500 µS/cm to resources in CP 2. In this context, δ 18 O-data and electrical conductivity measurements propose that shallow groundwater at well J8 (Chacarillas alluvial fan) is the result of a mixing of a recharge provided from the lower stream segment of the Chacarillas river and the Puquio de Nunez area (CP 4) (Figure 5A).
Due to the apparent absence of the isotopic imprint of fractionated floodwater on shallow groundwater around alluvial fans, flash floods might not be in the long-term a significant source of recharge.The absence of 3 H in shallow groundwater at alluvial fans (<0.09TU) emphasizes this conclusion because modern rainwater exhibits 3 H values of 3 to 10 TU in the region [23].The respective shallow groundwater must, therefore, show residence times above 50 years.This finding implies that modern floodwater may be lost in majority to evaporation.Accordingly, earlier detected 3 H amounts in groundwater of the PdT were always associated with sites of agricultural activity and were attributed to isotope exchange with the atmosphere during watering and irrigation-return flow [3].However, a more detailed study of the different alluvial fans could provide more evidence for the presented conclusion.Recharge conditions may also vary from fan to fan.
Nevertheless, the understanding that at the Chacarillas alluvial fan there is no significant recharge occurring from floodwater is in contradiction with Houston [21].He argues-based on the interpretation of a hydraulic head rise at well J8 (total depth ~200 m b.g.l.)-that a regular in situ recharge from the alluvial fan takes place, with a minimum average annual recharge equivalent to around 200 L/s.However, the phenomenon of the hydraulic head rise at well J8 could be explained solely by a short-term hydraulic head propagation (pressure) induced by recharge occurring at the sloping stream bed and fan apex of the Chacarillas river (as depicted in the graphical abstract).Such a continuous propagation of pressure signals is controlled by the aquifers hydraulic diffusivity and was recently demonstrated to occur in the Quisma basin [10,52].Due to the low specific storage of confined aquifers, which receive a focused recharge at the basins margin (see graphical abstract), the pressure signal could propagate within weeks in the PdT Aquifer.This signal propagation could cause a slightly delayed water level rise at well J8 in response to distant recharge events, which can be easily misinterpreted as an in situ recharge (a hydraulic head time series analysis to demonstrate this effect exceeds the framework of this article but may be subject to a later study).
Another factor that supports the absence of alluvial fan recharge is the presence of thick (>20 m) unsaturated and clayey sediment layers at the alluvial fans of the Juan de Morales and Chacarillas basins.These unsaturated top layers were identified earlier by remote sensing imagery, drilling profiles and extensive TEM measurements carried out in the PdT [2,53].Clayey top sediments show resistivities in the magnitude of 700-1000 Ωm compared to ~20 Ωm of the saturated aquifer section and ~100 Ωm for sediments of the vadose zone [2].A regular percolation from alluvial fans (4-year return period of floods, Section 2) would need to result in thick vadose zones connecting clayey surface sediments at alluvial fans with the aquifer, which is not the case.This argumentation line yields another indication towards the conclusion that no regular and significant recharge occurs from the discussed alluvial fans.
Despite the vital recharge area Altos de Pica (~4000 m a.s.l., Section 5.3) in the Quisma basin (CP 4), it is proposed that the principal recharge facilitated into the PdT Aquifer occurs along the lower segments of rivers that discharge into the sedimentary basin.That is derived from the relation between catchment-specific δ 18 O values and their sample altitudes (Figure 5).For CP 2 and 5 the δ 18 O values of groundwater from the PdT Aquifer match with the δ 18 O values of river and spring water of the corresponding catchments at altitudes between ~2200 m and 1000 m a.s.l.(Figure 5A).This matching implies that the water was recharged at these altitudes because isotopic fractionation by evaporation stopped at some point downstream.Major north-south striking faults in the region could reduce a prominent lateral groundwater inflow from higher elevations which is supported by the presence of several springs along the Andean Precordillera (Figure 1B, Section 2, [10,17]).
However, due to a lacking evaporative enrichment trend in data of CP 3 (Figure 5C), here, recharge areas providing higher fractions of lateral groundwater inflow to the PdT Aquifer could include streambed sections at higher elevations (Figure 5C and graphical abstract).That is underlined by the fact that respective waters show highest D ex values for the region (Figure 3C), which implies a minor fractionation by evaporation along the river course and hence possibly higher recharge elevations.
Another factor that influences the D ex value is whether a basins geologic characteristics facilitate the opportunity for interflow to take place along the streambed (reduced evaporation of river water).The importance of interflow is demonstrated by the fact that D ex values from the Chacarillas basin are the lowest for the region (Figure 3C) and its host formations are in great parts impermeable Jurassic and Cretaceous quartzous sandstones (Chacarilla formation) or marine sedimentary rocks (Majala formation) [19].Whereas streams in the Juan de Morales basin run mainly above Oligocene to Miocene sedimentary rocks (Altos de Pica formation) with thin layers of ignimbrites (Tambillo ignimbrite) that tend to weather easily [10,54].
However, overall it is known that a focused recharge occurring beneath ephemeral streams is in many arid environments the principal recharge mechanism [55,56].The derived conclusions imply that diffuse recharge at higher elevation is-in the mentioned cases-of minor importance.The presence of thick dry deposits (up to 500 m thickness and 4000 m a.s.l.[10,19,54]) along the Andean Precordillera, that form ravines through which rivers discharge, supports that because a diffuse recharge would demand of infiltrating rainfall fractions to penetrate these dry layers.That is not the case.Hence, focused recharge along ephemeral rivers fed by direct runoff is much more likely.Only the Altos de Pica area appears to facilitate geological conditions for a substantial diffuse recharge along the Andean Precordillera to the PdT Aquifer, due to its relatively plainly deposited and well-fractured ignimbrites [8,10,57].However, further research is recommended to substantiate the elaborated conclusions.

Conclusions
Earlier introduced assumptions about why collectively examined groundwater isotope samples from the Pampa del Tamarugal (PdT) Aquifer plot parallel to the local meteoric water line need to be reconsidered.A correlation of isotope data from the PdT with data from the nearby Andean Altiplano excludes a variation of isotope values in precipitation due to a change of regional climatic conditions.Other mechanisms are proposed to explain the isotope characteristics of meteoric water in the PdT.These mechanisms comprise rainout (amount effects) and altitude effects at the Andean eastern windward side and a mixing of air masses with different isotopic compositions at the transition zone from Precordillera to Altiplano (western and eastern vapor sources).It is proposed that exceptionally high topographic elevations can function as barriers for eastern air masses and inhibit the mentioned vapor mixing processes.The described topographically controlled effects could lead to characteristic isotope value ranges in precipitations of relevant basins.Hence, there is evidence that rainwater precipitating over PdT catchments can show distinct δ 18 O value spectra that in total cover a range from −19 to −7‰.Once precipitated, meteoric waters undergo an isotopic evolution that involves evaporation-driven fractionation and subsurface groundwater mixing.There is, however, no evidence for hydrothermal isotopic water-rock interactions in the investigated cases.Overall, it is likely that similar mechanisms control isotope values of groundwater in other basins of the Atacama Desert.
Due to the observed basin-specific isotope characteristics in the PdT, it is possible to use δ 18 O and δ 2 H as tracers to delineate six major aquifer compartments.These compartments are consistent with flow directions indicated by regional groundwater level contour lines.The recharge areas of the mentioned compartments can be spatially correlated with different Precordilleran slope catchments.Furthermore, there is isotopic evidence that water recharged at Altos the Pica (~4000 m a.s.l.)-in the Quisma basin-passes the regional Longacho Hinge and provides a recharge into the PdT Aquifer.
However, isotopically (δ 18 O, δ 2 H, 3 H) there is no proof that regular flash floods reaching the alluvial fans of the Juan de Morales and Chacarillas catchments contribute significantly to groundwater recharge.In this context, there are indications that the main recharge facilitated by the Chacarillas and Tarapacá rivers is infiltrating along their lower stream segments between ~2200 and 1000 m a.s.l. as focused recharge.

Figure 1 .
Figure 1.(A) View on South America (black rectangle represents the extent of Figure 1B), (B) Topographical map of the study area and sample types (elevation profiles are discussed in Section 5.2).Water table contour lines indicate the regional groundwater flow regime of the PdT Aquifer and are based on Jica [2].Arrows indicate resulting groundwater flow directions.Not all springs in the region are displayed.

Figure 1 .
Figure 1.(A) View on South America (black rectangle represents the extent of Figure 1B), (B) Topographical map of the study area and sample types (elevation profiles are discussed in Section 5.2).Water table contour lines indicate the regional groundwater flow regime of the PdT Aquifer and are based on Jica [2].Arrows indicate resulting groundwater flow directions.Not all springs in the region are displayed.

Figure 2 .
Figure 2. (A) Kriging model based on δ 18 O values in river, ground and spring water and yielded meteoric water compartments 1-6 (limits marked by dashed lines), (B) Kriging model based on δ 2 H data of river, ground and spring water and implied groundwater flow directions.Note that the extension of the alluvial fans at CP 3 and CP 5 is not corresponding with the respective groundwater flow direction (Figure 1B).CP 1 (18 samples) represents western marginal samples where the first and third quartiles of δ 2 H values cover a range of −69-60‰ and for δ 18 O values plot between −8.7 and −7.5.The water of CP 1 appears to originate from the Aroma basin.CP 2 (48 samples) is associated with the drainage area of the Quebradas Tarapacá and Quipisca.δ 2 H data range dominantly between −82‰ and −75‰ but decrease with higher altitudes.The first and third quartile of δ 18 O values ranges from −10.5‰ to −9.9‰.Due to a lack of samples from the Quipisca basin, the Quipisca basin cannot be isotopically distinguished from the Tarapacá basin.Therefore, both basins are summarized into one compartment (CP 2).CP 3 (54 samples) forms a distinct plume of higher δ 2 H and δ 18 O values which is physically related to the Juan de Morales River.This group marks δ 2 H values in majority higher than −57‰.Its lower bound is formed by the δ 2 H value −69‰ and the δ 18 O value −9‰.CP 4 (223 samples) shows low δ 2 H values with down to −113‰ and δ 18 O values down to ~−13.5.It is associated with the Quisma basin.Recently it was demonstrated that the corresponding groundwater is recharged at the Altos de Pica recharge area (Figure2B)[10].A slight increase of δ 2 H and δ 18 O values can be found in some samples of the local and shallow Pica aquifer, which is

Figure 2 .
Figure 2. (A) Kriging model based on δ 18 O values in river, ground and spring water and yielded meteoric water compartments 1-6 (limits marked by dashed lines), (B) Kriging model based on δ 2 H data of river, ground and spring water and implied groundwater flow directions.Note that the extension of the alluvial fans at CP 3 and CP 5 is not corresponding with the respective groundwater flow direction (Figure 1B).
CP 6 corresponds to the Altiplano basins Salar del Huasco and Laguna Lagunillas (91 samples of the Salar del Huasco basin and three from the Laguna Lagunillas basin).The first and third quartile of δ 2 H values plot at −102‰ and −96‰ respectively.δ 18 O values range dominantly between −13.3‰ and −12.1‰.

Figure 3 .
Figure 3. (A) Box plot of δ 2 H data by compartment, (B) Box plots of δ 18 O data by compartment, (C) Box plot of Dex values by compartment, lower values are an indication for a higher isotope fractionation by evaporation (higher evaporative loss) (no significant difference in Dex values of precipitations of north-western and east-north-eastern air masses).

Figure 3 .
Figure 3. (A) Box plot of δ 2 H data by compartment, (B) Box plots of δ 18 O data by compartment, (C) Box plot of D ex values by compartment, lower values are an indication for a higher isotope fractionation by evaporation (higher evaporative loss) (no significant difference in D ex values of precipitations of north-western and east-north-eastern air masses).

Figure 4 .
Figure 4. (A) δ 18 O against δ 2 H (all 367 samples of the PdT) and resulting trend line, (B) Collection of 94 river, groundwater and spring samples from the Salar del Huasco basin and resulting trend line (local evaporation line) (P = precipitation).The predominant δ 18 O value range in P of the Salar del Huasco basin (~−17-14‰) is derived by the intersections of the lower and upper bound of the local evaporation line with the local meteoric water line (LMWL).The stable isotope value range in P for the Salar del Huasco basin varies due to amount, altitude and continental effects.The first and third quartiles of δ 18 O data of weighted rainwater in the Salar del Huasco basin plots consistently at −17‰ and −13‰.

Figure 4 .
Figure 4. (A) δ 18 O against δ 2 H (all 367 samples of the PdT) and resulting trend line, (B) Collection of 94 river, groundwater and spring samples from the Salar del Huasco basin and resulting trend line (local evaporation line) (P = precipitation).The predominant δ 18 O value range in P of the Salar del Huasco basin (~−17-14‰) is derived by the intersections of the lower and upper bound of the local evaporation line with the local meteoric water line (LMWL).The stable isotope value range in P for the Salar del Huasco basin varies due to amount, altitude and continental effects.The first and third quartiles of δ 18 O data of weighted rainwater in the Salar del Huasco basin plots consistently at −17‰ and −13‰.
depicts the δ 18 O values of the springs for different sampling campaigns from 1967 to 2014.The temporal resolution indicates that stable isotope values did not vary significantly over the last decades (average standard deviation for all springs: δ 18 Oasd = 0.54‰, δ 2 Hasd = 4.2‰).No strong long-term or seasonal variations can be observed.The spring site with the highest δ 18 O variation are the Mamina springs in CP 3.

Figure 6 .
Figure 6.Time series of δ 18 O measurements at different spring sites in the study area and average standard deviation in ‰VSMOW (1967-2014).
depicts the δ 18 O values of the springs for different sampling campaigns from 1967 to 2014.The temporal resolution indicates that stable isotope values did not vary significantly over the last decades (average standard deviation for all springs: δ 18 O asd = 0.54‰, δ 2 H asd = 4.2‰).No strong long-term or seasonal variations can be observed.The spring site with the highest δ 18 O variation are the Mamina springs in CP 3.
depicts the δ 18 O values of the springs for different sampling campaigns from 1967 to 2014.The temporal resolution indicates that stable isotope values did not vary significantly over the last decades (average standard deviation for all springs: δ 18 Oasd = 0.54‰, δ 2 Hasd = 4.2‰).No strong long-term or seasonal variations can be observed.The spring site with the highest δ 18 O variation are the Mamina springs in CP 3.

Figure 6 .
Figure 6.Time series of δ 18 O measurements at different spring sites in the study area and average standard deviation in ‰VSMOW (1967-2014).

Figure 6 .
Figure 6.Time series of δ 18 O measurements at different spring sites in the study area and average standard deviation in ‰VSMOW (1967-2014).

Figure 7 .
Figure 7. Temperature logs of wells J5, J8, JF and LC2 and depth of depth-specific stable isotope samples (blue arrows mark sample depth of depth-specific samples, italic numbers indicate measured δ 18 O values).

Figure 7 .
Figure 7. Temperature logs of wells J5, J8, JF and LC2 and depth of depth-specific stable isotope samples (blue arrows mark sample depth of depth-specific samples, italic numbers indicate measured δ 18 O values).
O values from ~−17‰ to −13‰ (in accordance with the assessment in Section 5.1) and in CP 4 from ~−19‰ to −13‰.The precipitations alimenting meteoric waters in CP 2 would show δ 18 O values of dominantly ~−15 to −11‰.The respective signals lose in intensity downstream-below the rainfall limit of ~2500 m a.s.l.-where groundwater mixing processes become more significant and hence respective values are being averaged (Figure 8A).

Hydrology 2018, 5 , 3 15 of 22 Figure 8 .
Figure 8. (A) Approximate initial δ 18 O values of meteoric waters in the study area prior to evaporation-driven fractionation.The map is based on the interpolation data of the kriging models for δ 18 O and δ 2 H (Figure 2A,B), arrows indicate suggested air mass contributions, (B) Schematic sketch of proposed topographic controls of isotope values in precipitations contributing to recharge in the PdT (isotope data of rain originates from Fritz et al. [3], Uribe et al. [8], and Aravena et al. [23]).

Figure 8 .
Figure 8. (A) Approximate initial δ 18 O values of meteoric waters in the study area prior to evaporation-driven fractionation.The map is based on the interpolation data of the kriging models for δ 18 O and δ 2 H (Figure 2A,B), arrows indicate suggested air mass contributions, (B) Schematic sketch of proposed topographic controls of isotope values in precipitations contributing to recharge in the PdT (isotope data of rain originates from Fritz et al. [3], Uribe et al. [8], and Aravena et al. [23]).

Table 2 .
Depth-specific samples taken for this study (locations displayed in Figure2B).

Table 2 .
Depth-specific samples taken for this study (locations displayed in Figure2B).

Table 3 .
Isotope data of shallow groundwater at alluvial fans (alluvial fan = af, EC = electrical conductivity) (position of wells displayed in Figure2B).