The provenance and persistence of the perennial Río Loa in the Atacama Desert: links between crustal processes and surface hydrology

1 Introduction

The perennial flow of the Ríos Loa and Salado crossing the Atacama Desert, considered the driest place on earth, is anomalous (Figure 1). There are no other rivers that cross the Atacama and reach the ocean between 19 and 27^oS, a distance of over 700 km.

FIGURE 1. The Río Loa catchment. (A) extent, location and physical features. The DEM shows the principal sedimentary basins as near flat surfaces. The yellow contour is the 120 mm/a isohyet. PP Puente Posada, CB Calama Basin, TB Turi Basin, SP San Pedro Basin, LB Linzor Basin, ET El Tatio, PG Pastos Grandes Caldera, LC Laguna Colorada. (B) Río Loa catchment, surface water divide, principal abstractions and flow gauging/sample locations. Red dots represent major springs, pink polygons minor springs (a- Angostura, b- Laguna Inca Coya, c- Río Salado, d- La Teca, e− Chitor) and green triangles rain gauge stations.

The Río Loa system is a rich ecological resource within the desert and has acted as a corridor for Paleo-Indian cultures since at least 13,000 years BP (Aldenderfer, 1989; Núñez and Grosjean, 2003; Moreno et al., 2009; Latorre et al., 2013; Pestle et al., 2015; de Porras et al., 2021). Since the beginning of the 20th century water demand has increased as much as ten-fold as a consequence of mining, agricultural and municipal demands on the river and groundwater resources within its catchment. The resource is now in danger of being so exploited as to cease functioning at some point in the future. Sustainable management of these resources needs to be based on detailed knowledge of the system as a hydrological entity. Whilst considerable attention has been paid to the groundwater resources of the Río Loa catchment, little has been paid to the surface water system, except for the pioneering isotopic work of Magaritz et al. (1989), Aravena and Suzuki (1990) and Aravena (1995).

Surface water runoff is initiated by precipitation, but over 90% of annual precipitation in the catchment of the Río Loa occurs at high altitude in the summer months of December to February (DJF) (Garreaud et al., 2009; Espinoza et al., 2020; Arias et al., 2021; Houston and Latorre, 2022). The river is maintained throughout the rest of the year by interflow and baseflow which originate as transient release from the vadose zone and long-term discharge from groundwater respectively (Horton, 1933; Barnes, 1939; Rorabaugh, 1964). For any gauging station, only active aquifer storage above the elevation of the gauging station contributes to baseflow; once the groundwater gradient decays to this putative level, flow ceases. However, this does not preclude deeper groundwater entering the hydrological cycle, where topographic and flow gradients are steep (Tóth, 2009; Frisbee et al., 2017; Somers and McKenzie, 2020), nor the possibility that the hydraulic basement as currently envisaged (Figure 2) (Naranjo and Pascoff et al., 1982; Herrera et al., 2021) may be permeable on long timescales, as a result of diffusivities greater than zero, and as changes in baselevel, gradients or recharge rates occur (Manning and Ingebritsen, 1999; Ingebritsen and Manning, 2002; Diamond et al., 2018). Crucially, it is hypothesized that lithospheric sources such as geothermal, magma volatile release and/or subduction dehydration may enter the hydrological cycle associated with arc volcanism (Aravena and Suzuki, 1990; Godfrey et al., 2019; Godfrey et al., 2021). Several aquifers have been identified within the catchment (Figure 2): high altitude intra-arc, structurally and volcanically controlled aquifers (Houston, 2009; Urrutia et al., 2019), Pre-Andean Calama and Turi Basin volcano-sedimentary aquifers at 2,300–3,200 m (Houston, 2004; Jordan et al., 2015), and forearc hosted aquifers of the Central Depression (JICA, 1995; Houston et al., 2001; Rojas et al., 2010; Samuel et al., 2019; Viguier et al., 2019; Jordan et al., 2022).

FIGURE 2. Catchment geology. (A) schematic hydrostratigraphic chart for major forearc and Pre-Andean basins in the Río Loa catchment. The aquifer conditions for Quillagua-Llamara and Puente Posada are based on outcrop lithology, geological mapping and inferred hydrogeological conditions. The Calama and Turi Basins are based on geological mapping, drilling, geophysical logging and well test data. Generalized river incision in relation to formations for each area is shown. Published sources are given for each basin. (B) lithology and structure of formations hosting the upper, phreatic aquifers within the surface water catchment of the Ríos Loa and Salado. Q- Quillagua, PP- Puente Posada, CB- Calama Basin, TB- Turi Basin. AFZ- Atacama Fault Zone, PCFZ- Pre-Cordillera Fault Zone,WF- West Fissure Fault Zone, LF- Loa Fault, AFS- Ayquina Fault System. ACL- Antofagasta-Calama Lineament, COT- Calama-Olcapato-El Toro Lineament, LCL- Lipez-Coranzuli Lineament. Pink line- Calama Basin monocline. Red triangles-volcanoes.

Additionally, several recent studies (Jordan et al., 2002; Corenthal, 2016; Munk et al., 2018; Moran et al., 2019; Boutt et al., 2021) have suggested that so-called orogenic groundwater (here termed far-field) flow that is sourced from beyond the surface water catchment boundaries, contributes significant fluxes to both surface and groundwater along the western slope of the Andes. The most extensively reported example is to the south, in the Salar de Atacama where both the water and salt fluxes within the catchment appear to be out of balance unless sources outside the catchment are taken into consideration (Moran et al., 2019; Boutt et al., 2021). To the north, Magaritz et al. (1990), Aravena (1995), Jayne et al. (2016) and Alvarez-Campos, et al. (2022) have all suggested that a proportion of the groundwater in the Pampa Tamarugal (Central Depression) originates from the high Andes as infiltrating surface water, transmitted via deep seated fissures open to groundwater flow. Elsewhere, smaller scale cross-catchment boundary flow has been widely reported (Anderson et al., 2002; Montgomery et al., 2003; Houston, 2006c; Rissmann et al., 2015).

The present study is based on a field program during 1999–2001 that measured twenty two surface water flows, and seven springs at different times of the year, at selected locations throughout the whole catchment, concurrently obtaining water samples for major ion and isotope analysis. Daily precipitation and river flow data for the period 1976–2001 were obtained from the Dirección General de Aguas (DGA) for three key locations, in order to place the study data in context, whilst groundwater and geothermal data from the literature are also used to provide background and are referenced in the text. The study integrates physical, chemical and isotopic aspects in order to explain the hydraulic and hydrogeochemical processes throughout the catchment. In particular, the origins of the system in the volcanic arc are examined in relation to lithospheric processes. The results have wider implications for fluvial activity in arid regions, subduction zones worldwide, and specifically for the sustainable management of water resources in northern Chile.

2 The Río Loa catchment

2.1 Morphological development

These hydrological systems have their origin in the uplift of the Andes. Figure 1A is a DEM of the catchment showing its major morphological features. Shallow slopes of <2° signify basins with volcano-sedimentary infill hosting the major aquifers throughout the catchment, whilst slopes of 2–4° are generally indicative of colluvium and alluvial fans where significant recharge takes place under favourable conditions (Houston, 2002). The thalweg of the Río Loa is around 450 km long from its catchment boundary (maximum elevation ∼6,000 m) to the ocean encompassing approximately 33,570 km². The river presents three segments, each with coherent gradients reflecting successive phases of tectonic uplift (Hoke et al., 2007). At the upper bound of each segment lie three deeply incised canyons, 110–170 m deep in the Western Cordillera, 160–190 m deep through the Pre-Cordillera and 450–490 m deep through the Coastal Range. Each canyon ends with a significant knickpoint that appears to have stalled in its upstream migration where it encounters resistant bedrock. The upstream lip of the gorges is the effective base level for active groundwater storage in the basins above, which are up to 1 km or more deeper than the canyon lip (Jordan et al., 2015; Jordan et al., 2022), suggesting that there may be considerable stagnant or slowly overturning storage within them. Localised dewatering of these aquifers orthogonal to the river occurred during incision distorting adjacent groundwater flow nets (Houston, 2006c).

The forearc and Pre-Andean basins were infilled with clastic sediments, together with repeated pyroclastic deposits, and progressively overlain by carbonate, sulphate and halite evaporites between the Eocene and Late Miocene (De Silva, 1989; Lamb, 2016; Garzione et al., 2017; Horton, 2018; Jordan et al., 2022). The Calama Basin was largely endoreic during this period (May et al., 1999; Blanco et al., 2003; Mpodozis et al., 2005), but beginning 10–7 Ma, sedimentation spilled over the Pre-Cordillera through the Calama Gap to the Central Depression with deposition of the Opache Formation (May et al., 1999; May et al., 2005; Hoke et al., 2007; Jordan et al., 2010; Jordan et al., 2015). However, surface water flow beyond Calama to the Central Depression was limited, as demonstrated by shallow gradients (<1°) and transient ponding in the Calama gap (Houston et al., 2001; Jordan et al., 2010; de Wet et al., 2015; de Wet et al., 2020). Uplift of the Western Cordillera due to tectonic activity and tilting took place during the Middle Miocene through the Pliocene (Hartley and Evenstar, 2010; Jordan et al., 2010; Evenstar et al., 2016) leading to incision of the Río Toconce on the western slope of the Cordillera between 3.9 and 1.1 Ma (Houston, 2006c) and the paired Ríos Loa and San Salvador as marginal drainage channels through the structural gap in the Pre-Cordillera and fully connecting the Calama Basin with the forearc hydraulically post ca. 3 Ma (May et al., 1999; Hartley and Evenstar, 2010).

The Central Depression between Quillagua and the Salar de Llamara was also endoreic from the Eocene, with a significant Plio-Pleistocene evaporite sequence overlying lacustrine carbonates and clastics. Accommodation space in the basin was filled by this sequence, which, coupled with tectonic activity caused breaching of the Coastal Range at its lowest point allowing the Río Loa to incise the canyon that led to a new (oceanic) baselevel after ca. 0.3–0.2 Ma (Sáez et al., 1999; Sáez et al., 2012; Mohren et al., 2020; Jordan et al., 2022).

2.2 Geology and aquifers

The geology of these basins has been studied extensively (Ferraris, 1978; Skarmeta and Marinović, 1981; Naranjo and Pascoff, 1982; Marinović and Lahsen, 1984; Jensen, 1992; Sáez et al, 1992; Medina et al., 2012; Duhart et al., 2018; Vásquez et al., 2018; May et al., 1999; Houston, 2004; May et al., 2005; Blanco et al., 2009; Jordan et al., 2015; Jordan et al., 2022) so that a hydrostratigraphic chart can be presented (Figure 2A), in addition to a sketch map of the lithology hosting the upper, phreatic aquifers of the catchment (Figure 2B). Here the geology is only briefly reviewed with respect to the main aquifers.

At high elevations volcanic edifices and ignimbrite sheets of Miocene to Pleistocene age have steep slopes giving rise to rapid runoff as well as infiltration into a series of perched aquifers consisting of unconsolidated ashes, tuffs and blocky lava interbeds. Within the Western Cordillera a number of aerially small intra-arc sedimentary basins exist at elevations between 3,700 and 4,500 m. Despite their apparently limited size these basins and associated porous volcanic edifices represent intra-arc aquifers which are of particular interest since they receive significant recharge under present day conditions (e.g., Linzor Basin - Houston, 2009), and are a major source of supply for nearby copper mines (e.g., San Pedro-Inacaliri Basin).

Mesozoic rocks with low matrix permeability dominate the Pre-Cordillera. Between the high Cordillera and Pre-Cordillera a series of moderately sized Pre-Andean basins occur: the Calama and Turi Basins. Houston (2004) identified upper and lower aquifers in the Calama Basin separated by an intermediate aquiclude-aquitard, and recent studies have confirmed this (Blanco and Tomlinson, 2009; Tomlinson et al., 2010; Jordan et al., 2015; Herrera et al., 2021). The upper aquifer in the Turi and Calama Basins is hosted by the Toconce, Chiquinaputo and Opache Formations, largely fluvio-clastic and palustrine-lacustrine carbonate sediments (Figure 2). The lower aquifer in the Calama Basin is hosted by the coarse to fine-grained clastic sediments of the Jalquinche, Lasana, Yalqui and possibly the Calama Formations, and by an unnamed formation (assumed to be an extension of the Jalquinche Formation) into the Turi Basin. The aquitard separating the upper and lower aquifers in the Calama Basin is a mixture of welded ignimbrite (principally the Sifon Ignimbrite) and clay-rich facies of the Jalquinche Formation (Houston, 2004; Jordan et al., 2015; Herrera et al., 2021).

Between the Pre-Cordillera and Coastal Range lies the extensive Central Depression (Pampa Tamarugal) with several, serially connected basins rising northward (and hence draining south toward the Río Loa) from 700 to 1,300 m, as a result of tectonic activity on cross-strike faults (JICA, 1995; Rojas et al., 2010; Samuel et al., 2019; Viguier et al., 2019; Jordan et al., 2022). The Quillagua Formation, largely comprised of palustrine carbonates together with diatom-rich silts and marls is equivalent to the upper aquifer of the Calama Basin, whereas the Quebrada Amarga Formation comprised of carbonates and clastic sediments (Sáez et al., 1999; Vásquez et al., 2018; Jordan et al., 2022) represents a semi-confined lower aquifer in the vicinity of Quillagua-Llamara, becoming unconfined in the vicinity of the Río Loa incised valley.

A number of structures conceivably impact groundwater, and potentially surface water flow by extrapolation (Figure 2B). A series of long-lived north-south along-strike faults, at various times in their history normal, reverse, and strike-slip occur sequentially from west to east: the Atacama Fault Zone bounding the fore-arc basins (Brown et al., 1993), the Pre-Cordillera Fault and the West Fissure, bounding the Pre-Andean Basins (Reutter et al., 1993; López et al., 2020), the Loa Fault (Jordan et al., 2015) and the Ayquina Fault system bounding the Turi Basin (Houston, 2004). Three deep seated, crustal lineaments cross the Calama Basin obliquely: The Antofagasta- Calama (NE-SW), and both the Calama-Olcapato-El Toro and the Lipez-Coranzuli lineaments (NW-SE) (Palacios et al., 2007; Lanza et al., 2013; Norini et al., 2013; Tibaldi et al., 2017). Finally, the Chiu Chiu monocline, a sinuous monocline dividing the northwest part of the Calama Basin from the southeast may also affect groundwater flow (Blanco and Tomlinson, 2009; Jordan et al., 2015).

2.3 Climate

Convective activity generated by the South American Summer Monsoon (SASM) brings low-latitude summer precipitation to the basin, while the Southern Hemisphere Westerlies (SHW) bring winter precipitation from high latitudes. Between these systems lies the Andean Dry Diagonal (ADD), a zone of minimal precipitation at the core of the Atacama Desert (Garreaud et al., 2003; Arias et al., 2021; Houston and Latorre, 2022). Precipitation is strongly influenced by the Andes which create a rainshadow on their western slopes: 3.5 mm at Calama (2,300 m), compared with ∼120 mm/a at 4,000 m (Houston and Hartley, 2003). The Río Loa catchment lies at the core of the Atacama Desert, largely to the north of the ADD where ca. 90% of rainfall occurs in the summer months (DJF) compared with >90% south of the ADD during winter (AMJJAS) (Garreaud et al., 2003; Espinoza et al., 2020; Arias et al., 2021; Houston and Latorre, 2022; Valdivielso et al., 2020). Although it is recognized that winter precipitation occurring as snow is often under recorded (Vuille and Ammann, 1997).

High-intensity long-duration convective precipitation is required to generate extreme flood runoff and groundwater recharge events, which occurs sporadically over the Western Cordillera and Pre-Cordillera (Houston, 2005; Houston, 2006b; Poveda et al., 2020). Such extreme events are augmented by quasi-periodic circulation modes: El Niño-Southern Oscillation, Atlantic sea surface temperature, Pacific Decadal Oscillation, and the Madden Julian Oscillation (Lenters and Cook, 1997; Takahashi et al., 2011; Boers et al., 2015; Valdés-Pineda, 2018; Sun et al., 2019; Espinoza et al., 2020; Vuille et al., 2000).

Of particular importance to past groundwater recharge are the wetter precipitation episodes that occurred during the Holocene and Late Pleistocene (de Porras et al., 2021; González-Pinilla et al., 2021). Wetter periods tend to have recurred at intervals of 10²–10³ years, although not always with the same intensity. The Central Andean Pluvial Event (CAPE; 15.9–7.5 ka BP) led to widespread perennial discharge in what are now ephemeral channels (Latorre et al., 2002; Houston, 2005; Nester et al., 2007; Gayo et al., 2012), and likely contributed to contemporaneous groundwater recharge events, the latest arguably in the Late Holocene (5.0–4.6, 3.2–2.1 and 1.4–0.7 ka BP) (Tully et al., 2019; González-Panilla et al., 2021).

3 Hydrology

3.1 Flow regime

The perennial river network is shown in Figure 1B. Three principal high altitude subcatchments; the upper Loa, San Pedro and Salado are all fed from aquifers during the dry season. A middle section transfers water from these catchments through the gap in the Pre-Cordillera, from where the lowest section crosses the hyper-arid forearc basin, finally to discharge into the Pacific Ocean. A number of tributaries are present, especially in the Salado catchment, as well as minor springs throughout which are inferred to represent bank storage return water (see Supplementary Material S1).

Using data from the Chilean national catchment dataset: CAMELS-CR2, (Alvarez-Garreton, et al., 2018), daily rainfall and stream flows for three key stations are plotted in Figure 3A for the years 1975 through 2002, and their details are given in Table 1. The rainfall station at Linzor for the Río Salado (at El Sifon), and Ujina for the Río Loa (above Lequena) are located near the respective upper catchment boundaries, at some distance from the stream gauging site but illustrative of catchment precipitation, whereas that for Yalquinche (Calama) is adjacent to the gauging site. Low flow for these stations is defined as the lowest daily flow recorded at the end of the dry season (usually November), and baseflow is defined as all flow that occurs after direct runoff and interflow have ceased and before the next flood.

FIGURE 3. Flow time series. (A) stream hydrographs and precipitation at three key locations within the Loa catchment, with time-invariant flow rate arrowed. (B) recession analysis for Río Salado at El Sifon (locality 17).

TABLE 1. Catchment characteristics, estimated recharge and aquifer properties for the middle and upper Ríos Loa and Salado derived from recession analysis of baseflow.

Summer rainfall events give rise to flood peaks with short lag times in all stations, indicating a rapidly responding flash flood regime typical of arid environments (see Supplementary Material S2). Lequena and El Sifon have flash flows correlated with their respective rainfall stations, but it is noteworthy that Calama flood events are mainly out of phase with local precipitation demonstrating that summer precipitation-runoff events in the Western Cordillera influence the river throughout its course (Houston and Hartley, 2003; Houston, 2006b; Houston and Latorre, 2022). The flood peak is always short-lived, usually no more than a few days after the rainfall peak, although there may be more than one event per season. Thereafter, interflow, the release of water taken into bank storage during intense rainfall and flood flow, dominates, and depending on the magnitude and intensity of the rainfall event, normally lasts no more than 1–2 months. After bank storage has been depleted the remainder of the dry season flow is supported by baseflow. During this period frequent low-amplitude, low-frequency increases in flow occur during the winter months April to September. The cause of these small (<0.1 m³/s) increases in flow are postulated to be due to one or more of the following: a) winter precipitation including snowmelt, b) reduced winter evaporation, or c) human intervention. The data to support a) or b) of these options is inconclusive. Human intervention (c) in the form of irrigation takeoff and return, as well as abstraction and reservoir release on a daily basis are a mostly unknown and unrecorded quantitative feature of streamflow. Nevertheless, Figure 3B shows that baseflow and recession can still be clearly identified.

During the period 1977–2002 low flows never drop below a threshold specific to each gauging station (Figure 3; Table 1), and the recession is asymptotic to the low flow threshold value and never to zero. The threshold is stable over time arguing against it being a precipitation related source which would be variable year-on-year. This time-invariant component of baseflow is theorized to have its source within the lithosphere from several possible causes: a) slab dehydration (Ingebritsen and Manning, 2002), b) magma volatile exsolution (e.g., Urann et al., 2022), c) hydrothermal discharge (Giggenbach, 1992), d) deep recirculation of meteoric water (Diamond et al., 2018), or e) far-field (extra catchment) recharge (Boutt et al., 2021). These possibilities are explored further below.

Based on the lowest flows in the upper reaches of the Loa and Salado subcatchments (localities 17 and 12), the combined time-invariant flow is 0.717 m³/s, only 0.03 m³/s (4%) less than that in the Loa at Calama (between localities 7 and 9). The difference may be accounted for by additions from the Río San Pedro, a left bank tributary of the Río Loa joining between localities 11 and 12, which has average low flows of <0.1 m³/s at the confluence. Critically, the lack of additional flow in the Pre-Andean basins signifies that the upper Ríos Loa and Salado supply virtually all low flow at Calama.

Data from October 1916 reported by Jordan et al. (2015) indicate that in the era before large scale abstraction took place, baseflow was 1.0–2.5 m³/s greater throughout the middle Loa (Calama Basin), and possibly as much as 3 m³/s greater in the upper Loa basin.

3.2 Recession analysis

Analysis of the recession curve allows for the estimation of groundwater storage depletion, or contribution to baseflow over time. Baseflow is derived from active storage within an aquifer and is controlled by the physical properties of the aquifers (Szilagyi et al., 1998; Wittenberg, 1999; Troch et al., 2013). Analysis of the recession curve assumes no significant anthropogenic interventions such as reservoirs, abstraction or return flows (Singh and Stall, 1971; Sujono et al., 2004; Bloomfield et al., 2011; Yao et al., 2021). Such interferences do occur at several locations on the Ríos Loa and Salado (Figure 1B), although during the period of analysis (1977–2002) the slope of the recession limbs remain virtually constant for each station, suggesting that artificial regulation is either short-lived or near-constant in time. With the possible exceptions of the small variations as noted above, a simple tripartite division of the receding limb into flood runoff, interflow, and groundwater baseflow can be made. For the purpose of recession analysis, the storage reservoir is assumed to be a single horizontal, uniform, homogeneous aquifer with no bedrock leakage (Hall, 1968; Tallaksen, 1995; Troch et al., 2013). Since the recession occurs over a long period dry season (>280 days) the initial critical time (change from interflow to baseflow) becomes less important (Anderson and Burt, 1980; Szilagyi et al., 1998; Sujono et al., 2004) and the recession in its simplest form (Barnes, 1939) is then given by:

where Q₀ is the initial flow, Q_t the flow at time t after the initial flow and before any further flood flow, and k is the recession constant. By integration this leads to the storage equation where S (the volume remaining in active storage after time t):

Furthermore, at long times the aquifer hydraulic diffusivity (D) can be calculated from the recession constant thus:

where L is the length of the perennial river from the gauging station to the divide and k is the recession constant (Rorabaugh, 1964; Szilagyi et al., 1998). The diffusivity (D) of the source aquifer is a function of its properties:

where K is the hydraulic conductivity, b the aquifer thickness and Sy is the effective porosity (more correctly the specific yield). Kb is also equal to T, the transmissivity. Finally, the magnitude of the year-on-year displacement of the recession curve at the moment when interflow ceases and baseflow initiates is equivalent to the annual recharge of the storage reservoir (Meyboom, 1961; Szilagyi et al., 1998; Chen and Lee, 2003).

Evapotranspiration losses are additional to measured baseflow. Above 3,000 m widespread Puna and high Andean vegetation belts reduce infiltration and enhance evapotranspiration, as do spring-fed wetlands (Vilagrán et al., 1981; Houston, 2006a; Suárez et al., 2023). In the middle and lower reaches, evaporation is limited to in-channel areas, and transpiration to riparian vegetation. Estimated evaporation rates based on elevation (Houston, 2006a) and isotopic fractionation provide an estimate (see Supplementary Material S3) that varies from 0.04 to 0.12 m³/s in the upper catchments and 0.21 m³/s in the Calama Basin.

Throughout the catchments, the respective hydrographs show no changing impact from abstraction during the period of analysis, so that no allowance for abstractions is made in the recession analysis.

The Río Salado at El Sifon is used as an example of the technique (Figure 3B) which is then applied to all three streamflow stations (Table 1). When plotted in semi-logarithmic form baseflow recession yields a straight line, with the slope equal to the recession constant (k). Single, near-parallel late-stage recession curves support the concept of an ideal single linear reservoir (aquifer) with annual replenishment. However, given the rather noisy data (see above), a graphical technique for the derivation of the recession constant k, and the moment of baseflow initiation has been favored (Szilagyi et al., 1998; Sujono et al., 2004; Stoelzle et al., 2020).

Mean values of k for each subcatchment are given in Table 1 and vary from 0.9988 through 0.9994 per day. Such consistently high values of k are indicative of slow recession from large active storage reservoirs. Unusually, each year’s recession decays to the same flow rate, and this is true for all three gauging stations. Typically, wet years should end on higher flows and vice versa. This limiting flow may be labelled as time-invariant flow, whilst all baseflow above this rate is considered time-variant. Mean time-variable baseflow for the upper Ríos Loa and Salado is 0.04 and 0.05 m³/s respectively, whilst the time-invariant components (baseflow minus time-variant flow) averages 0.57 and 0.32 m³/s. Unit baseflows are 0.96 × 10⁴ m³/a/km² for the upper Río Loa and 1.48 × 10⁴ m³/a/km² for the Salado. Unit baseflow in the Río Loa at Calama (0.36 m³/a/km²) is significantly lower because it is generated in only a small proportion of the total catchment (above ∼3,000 m).

Active storage at year (dry season) end is ∼28 × 10⁶ m³ for the upper catchments combined, compared with 32 × 10⁶ m³ for the whole catchment at Calama, indicating 88% is derived from intra-arc aquifers above ∼3,000 m. This large buffer against stream desiccation in the dry season, is interpreted as the reason for the perennial nature of the Loa-Salado river system.

Diffusivity is high for all active storage reservoirs, implying that permeability will be high, whereas specific yield is likely to be low suggesting semi-confined conditions (see Supplementary Material S4), and indicating coarse clastic or pyroclastic material possibly with significant fracture flow.

Annual precipitation recharge for each subcatchment based on the recession curve displacement method gives values between 10 and 15 mm/a above 3,000 m, and less than 5 mm/a at 2,300 m (Table 1). These estimates are based on total areas, yet it is probable that not all areas allow infiltration to the water table. Based on assumptions regarding suitable lithologies that might allow infiltration applied to existing geological maps (Marinović and Lahsen, 1984), conjectured recharge areas for each subcatchment are given in Table 1, which lead to actual precipitation recharge rates of 16–23 mm/a above 3,000 m.

3.3 Instantaneous flow transect

During October 1999 stream flows were measured throughout the Loa and Salado Basins, with the objective of identifying and quantifying gaining and loosing reaches that might indicate transfers from or to groundwater sources respectively. Temporary gauging sites were located along an east-west alignment from the catchment boundary in the Western Cordillera (Río Toconce), then due west along the Río Salado and Río Loa to the point where the latter eventually turns north and then west to the ocean. This section was chosen to give a typical hydraulic profile of the western Andean slope but does not give the full flow conditions for the whole catchment. Measurements were made using a float and the velocity-area method (Dingman, 1993) on suitable straight and uniform reaches of ∼3 m. Three to five repeat measurements at each site suggest that errors were <7%. Although these measurements took 5 days to cover the whole catchment, the data effectively give an instantaneous picture of the surface water balance. The measurement sites are located in Figure 1B and the results are shown in Table 2 and Figure 4. In order to discount tributary inflows to and licensed abstractions from the transect analysis, the adjusted flow at each site was calculated as:

where Q_A is the adjusted flow, Q_O is the measured flow, Q_T is the tributary flow which accounts for off-transect surface water gains, and Q_L is the licensed abstraction on-transect. There are three tributaries whose contributions need to be discounted (Figure 4), as well as licensed abstractions at Linzor (0.05 m³/s), and Toconce (0.35 m³/s) that need to be accounted. At several points along the transect, notably Toconce, Ayquina, Calama and Quillagua un-monitored agricultural takeoffs exist. These are not included in the analysis since a) their magnitude is not known with accuracy, and b) much of the water is returned via small-scale canalization. A repeat flow inventory was carried out in March 2000, when the hydrographs show that the flow regime was in its interflow phase.

TABLE 2. Instantaneous water budget for the Río Loa-Salado transect during baseflow and interflow conditions.

FIGURE 4. Instantaneous flow profile. (A) measured flows during baseflow for Loa-Salado Andean west slope transect. (B) calculated gains and losses for baseflow conditions. (C) Measured flow during interflow. (D) gains and losses during interflow. Note that whilst baseflow lasts up to 45 weeks pa, interflow only lasts 2–4 weeks pa.

A plot of the baseflow along the transect is shown in Figure 4 which indicates that a significant increase in flow (groundwater discharge) occurs at high altitude above Linzor, amounting to 0.43 m³/s in October 1999.

Not included in the transect but measured in October 1999, outflow from El Tatio to the south was measured at 0.19 m³/s and to the north, discharge was measured at 0.195 m³/s in the upper Loa (localities 12+13). Thus the combined baseflow from intra-arc aquifers to the Río Loa catchment based on the instantaneous water balance was 0.815 m³/s at the end of the dry season, compared with a mean of 0.972 m³/s throughout the year estimated from the recession analysis. These estimates of flow from the intra-arc aquifers reinforce their importance as reservoirs supporting baseflow throughout the Loa catchment.

Downstream from Linzor between 340 and 350 km distance from the ocean, baseflow loses 0.4 m³/s to unnamed formations (potentially an extension to the Jalquinche Formation of the Calama Basin) underlying the Sifon Ignimbrite which is breached by erosion in the Toconce gorge. A small gain of 0.05 m³/s occurs in the vicinity of Toconce (330–340 km from the ocean) due to drainage from the last downstream exposure of the Jalquinche Formation. On the transect, the Calama Basin phreatic aquifers receive leakage of 0.21 m³/s from the Río Salado in the vicinity of Chiu Chiu and again to the Río Loa of 0.43 m³/s at Cascad (locality 7). The last loss (0.05 m³/s) before the ocean occurs in the vicinity of Quillagua to the Quillagua Formation. Finally, below Quillagua there is a small gain (0.01 m³/s) from the Salar de Llamara aquifers (here the Quebrada Amarga Formation). In summary, during the dry season baseflow is largely derived from intra-arc aquifers, with a small contribution from the Pampa Tamarugal in the lowest reach. On the other hand the river recharges aquifers in the Pre-Andean and forearc basins.

The interflow regime in March 2000 shows higher flow rates as expected (Table 2; Figure 4), but the pattern of gains and losses show several important changes. Intra-arc aquifers are still providing significant (inter-) flow (0.67 m³/s), with an additional 0.24 m³/s downstream at 330 km. Leakage (0.37 m³/s) to the presumed Jalquinche Formation occurs between 332 and 342 km and again at 329 km, (0.29 m³/s). However, throughout the Calama Basin from the point where rivers exit the Cordillera to the passage through the Pre-Cordillera, instead of losses there are gains to the river, reaching a maximum of 1.11 m³/s at Cascad, indicating a reversal of flow when compared with the dry season. In the Central Depression, losses are again observed from the Río Loa to the Quillagua Formation, but with increased rates (0.35 m³/s) compared with the dry season.

4 Hydrogeochemistry

4.1 Sampling and analysis

During October 1999, 22 baseflow samples were taken from surface waters at the gauging sites discussed above, as well as other selected locations. Surface water samples were passed through a 0.45 μm filter in the field before being measured for pH, temperature and subsequently sealed in airtight 1 L bottles before transport to the laboratory for analysis. In addition, volume weighted samples of precipitation were collected from six sites during the period December 1999 to March 2000. These samples were collected by the DGA from selected rain gauge stations in the Cordillera. At each station the rain gauges were decanted immediately after storm events into small airtight bottles which were sealed and stored at ∼4°C until the end of the wet season when samples were aggregated for each station before being sent to the laboratory for analysis. The snow sample from Cerro Aucanquilcha is a bulk sample collected on 25 March 2000, 300 m above the snowline and 10 days after the storm event of 15 March that produced the snowfall. The snow sample from Cerro Cupo is also a bulk sample taken on 5 June 2000, 200 m above the snowline and 3 days after the storm event of 2 June from which the snow derived. Laboratory analyses were undertaken by ALS in Antofagasta, using ion chromatography for anions (Cl, SO₄, NO₃), volumetric titration for HCO₃, atomic absorption spectrography for cations (Na, K, Ca, Mg, As) and colorimeter for SiO₂. The results of the analyses are reported in Tables 3 for precipitation, and 4 for surface waters. Ion balance errors for the surface water sampled averaged 4.2% (±2.3%, 1σ). PHREEQC (Parkhurst and Appelo, 2013) was used to calculate P_CO2, and saturation indices (SI) for calcite, dolomite and gypsum (see Supplementary Material S5).

TABLE 3. Precipitation chemistry.

At selected sites, samples were taken for stable isotope, ultra low-level ³H, δ¹³C_DIC and ¹⁴C_DIC analysis at the GNS Science Laboratory in New Zealand. Field sampling followed the same protocols as above, whilst samples for AMS ¹⁴C_DIC determination were stabilized with HgCl₂. Samples were stored in the dark at ∼4°C until transport to the laboratory. Stable isotope measurements are with respect to VSMOW with error limits of ±1.0‰ for δ ²H and ±0.1‰ for δ ¹⁸O. Ultra-low level ³H determinations have a detection limit of 0.05 TU. In-situ production of ³H by fission of ²³⁸U and/or ⁶Li in magmatic rocks has been calculated at << 0.05 TU (Goff and McMurtry, 2000), so no allowance for these sources of ³H need be made. δ¹³C_DIC measurements are with respect to VPDB with a precision of 0.2‰. ¹⁴C_DIC results are normalized to δ¹³C_DIC at −25‰. The ¹⁴C_DIC detection limit is ±1 pmC and determinations have a precision of 1%–2%. The results are presented in Table 5.

4.2 Precipitation chemistry

The summer rainfall season of 1999–2000 was significantly wetter than average due to the impact of cold phase ENSO (La Niña) conditions that persisted from the end of 1998 through early 2001. The possibility exists therefore that the sample results are biased toward wetter than average conditions. Nevertheless, since there are few other major ion analyses available for precipitation in the region, they are presented here as an important aspect constraining the hydrogeochemical system.

TDS concentrations vary from 12 to 26 mg/L, excluding Ascotan with a maximum of 62 mg/L. The latter shows high Na, Ca, SO₄ and Cl concentrations indicating contamination from deflated gypsum and halite in the nearby salar. Snow is significantly lower, <3 mg/L and depleted in all ions as a result of leaching during snowmelt in the time between deposition and sampling (Hornbeck, 1986). Hence snow samples do not reflect the original composition of the precipitation. Although ionic concentrations are around 15-20 times greater than Amazonian precipitation (Faure, 1997), ion ratios indicate a greater affinity with Amazonian precipitation than with seawater. This is consistent with an easterly source expected during La Niña conditions (Cai et al., 2020; González-Pinilla et al., 2021).

4.3 River water chemistry

Chloride evolution through the catchment resolves into a dual process when Cl is plotted against distance from the ocean (Figure 5A). The first process is a gradual increase in concentration in the Río Loa and its tributaries, due to either evaporation of baseflow or dissolution of evaporitic minerals along the flow path (or both). The second process in the Río Salado and its tributaries, points to a gradual dilution of Cl from 7,760 mg/L at the exit to the El Tatio hydrothermal system, to 2,510 mg/L at locality 16 just before the confluence with the Río Loa. This process is interpreted to be the result of dilution of Cl-rich thermal water by incremental discharge from moderate to low-Cl intra-arc aquifers. The molar ratio of Na:Cl is close to 1 in the Pre-Andean and Forearc basins, suggesting the dissolution of halite, but subtle differences may be noted elsewhere (Figure 5B). Within the Western Cordillera excess Na (Na:Cl >1) is likely due to silicate weathering, whereas Na depletion (Na:Cl <1) is found in the hydrothermal waters. Downstream, below Quillagua, the molar Na:Cl ratio (0.701) is somewhat lower than sea water which may suggest that leakage of salts from the Salar de Llamara brine is more likely than dissolution of deposited marine aerosols (Rech et al., 2003; Wang et al., 2014; Voigt et al., 2020).

FIGURE 5. Hydrochemistry of the Ríos Loa-Salado surface waters. All ratios are molar. Dashed arrows show direction of flow. The symbol key applies to Figures 5–9, and the blue blocks locate the Coastal Range, Pre-Cordillera and Western Cordillera (see Figure 1). Between these structural blocks lie the major basins: the Central Depression or Pampa del Tamarugal, and the Calama-Turi Basins.

The Na normalized plot of HCO₃ to Ca ratio in Figure 5C indicates that the majority of baseflow samples plot close to the silicate weathering field (Gaillardet et al., 1999). There are two outliers: locality 1 next to the ocean plots close to the evaporite field, and locality 22 which is the outlet from the El Tatio hydrothermal field. No samples plot close to the carbonate field, so calcite dissolution appears to be unimportant. High silicate loads associated with rock weathering and geothermal activity are particularly noticeable in the headwaters of the Ríos Loa, San Pedro, Toconce and Hojalar as variously indicated by high Na, K, Ca and Mg to Cl ratios, low HCO₃:SiO₂ ratios and temperatures above ambient (Figures 5B,D,E,G,J,O respectively). Within the Western Cordillera the molar Ca:Na ratio varies widely between 0.04 and 0.49, but after the rivers descend below 3,000 m, the Ca:Na ratio remains stable at around 0.1, suggesting that silicate weathering is limited to higher elevations and volcanic rocks. Based on the Ca:Na ratios within the Cordillera, the principal source rock is rich in plagioclase (An_30-90) and associated ferromagnesian minerals (high K and Mg) consistent with the petrology of the widespread Miocene-Pleistocene lavas and pyroclastic rocks that range from andesites to dacites with a high percentage of unstable glass and some olivine (de Silva, 1989; Godoy et al., 2014; Godoy et al., 2017; González-Maurel et al., 2020).

The corresponding carbonate system is illustrated in Figures 5H,J pH increases from slightly acidic rain (pH = 6.5, Table 3) to slightly alkaline (8.0 ±0.4, 1 σ) throughout the flow regime, so that dissolved inorganic carbon (DIC) is virtually all HCO₃ (Table 4). The zone between 3,500 and 4,500 m has an extensive cover of cushion plants and tussock grasses that respire CO₂ (Houston, 2009), causing soil water to be enriched in HCO₃. As the rivers descend from this elevation, they lose HCO₃, which coupled with Ca released from silicate weathering and/or deflated gypsum incorporated from nearby salars in the presence of P_CO2 well above atmospheric in all samples, quickly proceeds towards calcite saturation (Figure 5I, Supplementary Material S5). Downstream throughout the Calama and Turi Basins the rivers maintain calcite saturation with travertine and calcareous tufa being deposited along their course as degassing of CO₂ takes place. Dolomite saturation follows a similar path to calcite and is supersaturated in rivers after they leave the Cordillera.

TABLE 4. Surface water chemistry compared with groundwater chemistry (sources indicated).

The headwaters of the Ríos Loa, San Pedro and Toconce present molar SO₄:Cl of ∼1 (Figure 5K) decreasing downstream to ∼0.1, suggesting high elevation volcanic S sources, and/or gypsum aerosol deposition from nearby salars, notably Ascotan. Conversely, the Río Salado headwaters have comparatively low SO₄ given that phase separation at depth within the El Tatio hydrothermal system allows for the production of peripheral acid-sulphate springs (Giggenbach, 1978). However, the same author reports flow rates amounting to 0.012 m³/s (95% total) for the mature Cl springs and only 0.0006 m³/s (5% total) for the SO₄ springs, which would explain the low SO₄ concentration in the Río Salado. Downstream from El Tatio, dilution of Cl increases the SO₄:Cl molar ratio. Once the rivers leave the Cordillera, the SO₄:Cl ratio maintains a constant 0.1 as far as Quillagua. Thereafter a rise in sulphate from 770 mg/L to 1750 mg/L is mirrored by increases in Ca, Na and Cl, coupled with a decrease in HCO₃. This change may be associated with the dissolution of marine aerosol deposition (Rech et al., 2003) or groundwater inflow (halite and gypsum derived brines) from the Salar de Llamara via the Quebrada Amarga (Rojas et al., 2010; Viguier et al., 2019). In the vicinity of Calama (locality 8) low-volume, high SO₄ and SiO₂ loads discharge to the San Salvador. The source of this input is likely associated with leachate from mining waste rock and tailings dams, some 3 km to the northeast of the Río San Salvador source.

The constant Na, K and SO₄ to Cl ratios downstream throughout the Pre-Andean and forearc basins, coupled with the increase in gypsum saturation (Figures 5B,D,K,L respectively) may be attributed to evaporative concentration, and argues against significant halite dissolution. Gypsum saturation is never sufficient for precipitation from surface waters anywhere in the catchment, unlike calcite and dolomite that are saturated throughout the Pre-Andean and forearc basins. Both Ca and Mg show depletions relative to Cl and Na through these basins, suggesting that they are being removed by calcite precipitation or that cation exchange is taking place. Since Mg:Ca ratios are everywhere below 1, calcite would be the mineral precipitating under current conditions (Tucker and Wright, 1990; Alonso-Zara and Tanner, 2010; deWet et al., 2015).

Nitrate is generally not a significant component of the flow system (Figure 5M), being everywhere less than 0.4 mg/L, except in the Río San Salvador where it reaches a maximum concentration of 65 mg/L close to the outskirts of Calama (locality 8). Its source is observed to be the sewage outlet from the town. Leakage from the Calama sewage collection system permeates the phreatic aquifers in the Opache Formation underlying the town; 19 mg/L at Ojos de Opache (locality 6), and 0.6 mg/L at Cascad (locality 7). At Chacance (locality 5) the Río San Salvador has a NO₃ load reduced to 6 mg/L, whilst the Río Loa (locality 4) has increased to 3.5 mg/L. Such convergence in concentration between the rivers implies mixing occurring in the passage through the Pre-Cordillera incised valley, potentially through the Batea or Calama Formations where they are exposed in the river beds that are ∼5 km apart above Chacance.

Arsenic is found throughout the surface water system varying from 0.1 to 29 mg/L, with the highest concentration from the outflow of El Tatio where Munoz-Sáez et al. (2018) report up to 47 mg/L for the upper geyser field. Low levels of SO₄ argue against sulphides as a source for the As. Alternative sources include leaching from ignimbrites (the Sifon contains 16 ppm - Álvarez-Amado et al., 2022, and the Puripicar 23–36 ppm - Youngman, 1984), as well as silicates and carbonates in the crust and mantle. Indeed, Tapia et al., 2020 suggest that As concentrations are positively correlated with crustal thickness, some 75 km under the Western Cordillera at this latitude (Sodoudi et al., 2011). Plots of As:Cl ratios against thalweg distance (Figure 5N) show that the impact of As from El Tatio reduces downstream in the Río Salado as dilution occurs, and subsequent mixing with the Río Loa. Of interest are the high As:Cl ratios found in the Toconce-Hojalar and San Pedro catchments which, together with elevated temperatures, are considered to be indicative of more extensive geothermal activity.

Silica concentrations are generally above 50 mg/L throughout the catchment, reaching more than 80 mg/L and up to 169 mg/L, in subcatchments above 3,000 m where water temperatures are above ambient, suggesting the widespread influence of geothermal activity. With the rather low Ca:Na ratios, alteration of albite-rich felspar and ferromagnesian minerals to clay minerals, especially smectite is indicated.

Samples taken in March 2000 are representative of interflow (Table 4). Interflow demonstrates similar processes and fluxes to baseflow samples taken in the previous October, the main difference being a uniform dilution of ionic concentration of 16%–24% (see Supplementary Material S6). This compares with ∼36% measured flow excess during interflow which may suggest that interflow chemistry includes a small admixture of water from groundwater storage, although the significantly lower HCO₃ suggests a shorter residence time for precipitation in contact with soil CO₂ in high elevation soils.

4.4 Isotopes

Stable isotope analyses are reported in Table 5. Summer precipitation in 1999/2000 gave δ¹⁸O values between −9 and −14‰, with δ²H varying between −60 and −100‰. A plot of δ²H vs. δ¹⁸O (Figure 6) shows that the local meteoric water line (LMWL) is given by:

which, considering the limited dataset, is indistinguishable from the summer LMWL of Valdivielso et al. (2020). The δ¹⁸O lapse rate for precipitation is −1.3‰/km, somewhat lower than that found by Aravena et al. (1999) and Bershaw et al. (2016), but closer to that found by Villablanca (2009), and Valdivielso et al. (2022) for summer precipitation (−1.5 to −2.8‰/km), supporting the concept that mesoscale convection differences due to topography create variations in precipitation δ (Gat, 2010; Fiorella et al., 2015; Scheihing et al., 2017; Valdivielso et al., 2022).

TABLE 5. Surface water and precipitation isotope analyses.

FIGURE 6. Stable isotopes. (A) δ²H vs. δ¹⁸O showing LMWL and LELs for the Ríos Loa and Salado. (B) stable isotope relations at El Tatio hydrothermal field and their relation to parental andesitic waters based on data sources listed. (C) ³H to d-excess crossplot with additional data from Cortecci et al. (2005) and projected data boxes from Giggenbach (1978) for an assumed ³H value.

Surface waters diverge from the LMWL along paths indicative of evaporitic enrichment with a slope (Δs) equal to 4.8. The evaporation lines (EL) have dissimilar departures from the LMWL at −14‰ δ¹⁸O for the Río Loa and −11‰ for the Río Salado. Based on the lapse rate, Río Loa waters originate from precipitation at around 6,000 m, whereas Río Salado waters originate from ∼4,000 m. This is broadly consistent with the local watershed elevation and the cloud base in the respective source areas. Turi Basin spring waters straddle the LMWL (Figure 6A) and the Salado EL showing a significant meteoric component in these spring waters.

Progressive downstream enrichment of δ¹⁸O indicates evaporitic fractionation during overland and in-channel flow. Based on a simplified Rayleigh equation:

where ΔɛO_w-v is the kinetic enrichment factor between water and vapor, f is the residual water fraction (Clark, 2015), and with enrichment assumed complete at −55 δ²H, the evaporative loss for the Río Loa approximates to 26%, and 16% for the Río Salado (see Supplementary Material S3). For comparison, Figure 6B shows data from Giggenbach (1978) and Cortecci et al. (2005) for the El Tatio hydrothermal field. Acid sulphate springs show the effect of mixing volatiles from a deep sources with meteoric water, whilst the Cl waters indicate variably sourced meteoric waters mixing with hot fluids from an andesitic reservoir, and stored in aquifers at varying depths (Giggenbach, 1978; Cortecci et al., 2005; Munoz-Saez et al., 2018).

Deuterium excess (d) is a function of source area humidity and is useful in differentiating source area and processes. It is expressed as:

In Figure 6C the plot of d against ³H (as an indicator of water age) reveals hydrothermal waters as old or ³H depleted. Modern water ³H >2 TU and d is greater than +10‰, whereas andesite waters are expected to have d around −100‰ with depleted ³H. From a highly negative d, increasing values are assumed to be the result of incremental volumetric mixing with meteoric waters: hot Cl waters evolve into moderate Cl+HCO₃ waters, with the HCO₃ component either the result of CO₂ degassing from the hydrothermal plumbing system and/or soil CO₂ generation in recharge waters of near surface reservoirs. There is considerable overlap in the ³H depleted 0 d region indicating widespread mixing of hydrothermal and meteoric water beyond the limits of the El Tatio geyser field, especially in the Turi samples.

DIC appears to be bimodal in the Cordillera with samples from both the Loa and Salado subcatchments below 7 mmol/L or above 10 mmol/L (Table 5). Below 3,000 m river DIC declines from 8 to 4 mmol/L. δ¹³C_DIC results are comparably bimodal in the Cordillera (Figure 7A); values between −6‰ and +6‰ could indicate a buried carbonate source but the preferred interpretation here is that these values suggest volcanic CO₂. In either case, degassing and the initiation of calcite precipitation in both the Ríos Loa and Salado occurs around 3,000 m. δ¹³C_DIC values close to −20‰ are associated with high altitude (>4,000 m) sources in equilibrium with soil CO₂. Farther downstream, in the forearc, δ¹³C_DIC values approximate to atmospheric CO₂ equilibrium at 8‰ reinforcing their identification as residual flood flow.

FIGURE 7. Carbon and tritium isotopes. (A) surface water δ¹³C_DIC plotted against distance from ocean. (B) ¹⁴C_DIC plot vs. distance from ocean. (C) relationship between ¹⁴C_DIC, and ³H indicates mixing between meteoric and hydrothermal water. (D) four groups distinguishable: (I) meteoric, (II) El Tatio hydrothermal derivatives, (III) San Pedro-Paniri-Toconce geothermal waters, and (IV) delayed flood water (see text for discussion). Additional data from Godfrey et al., 2021: + upper Río Loa river and springs, x Calama Basin groundwater. Red figures in (D) are ³H values in TU.

The presence of ¹⁴C in waters may originate from plant respiration, buried carbonates or volcanic volatiles, evolving via mixing and decay during transit. Given this complexity, no attempt is made to apply any age dating model with this data, although ¹⁴C_DIC is still invaluable as a tracer. From values around 50 pmC, ¹⁴C_DIC decreases downstream in the Western Cordillera (Figure 7B) to <10 pmC, but increases again throughout the Pre-Andean basins to ∼60 pmC, probably due to the decreased influence of dead volcanic carbon in these aquifers (Godfrey et al., 2021). Values above 100 pmC in the forearc confirm their interpretation as delayed flood waters.

Tritium values in precipitation vary between 2.6 and 4.7 TU, representing modern background levels in the central Andes (Houston, 2006c; Moran et al., 2019; Basaldúa et al., 2022). The downstream changes in ³H mirror those of ¹⁴C_DIC: at higher elevations in the Cordillera, ³H decreases towards the Calama Basin (<1 TU) and then increase again in the forearc to >1 TU (Table 5).

The presence of low levels of ¹⁴C_DIC are incompatible with measurable ³H in the same sample (Figure 7C), so that mixing between waters must be invoked. This is complicated by the decreasing then increasing values of ¹⁴C_DIC and ³H as baseflow moves downstream, with the inflexion point at approximately 3,000 m (Figures 7A,B).

4.5 Discrimination of water groups

Four water groups may be recognized based on location, chemistry and isotopic signatures (Figure 7D; Figure 8). Group I waters are found in the San Pedro, Hojalar and Toconce (and probably upper Loa) subcatchments associated with near surface aquifers above 4,000 m. The waters of group I (localities 18, 19, 20) are characterized by relatively dilute TDS (<600 mg/L), molar HCO₃:Cl ∼0.5, and d-excess between −10 and +10‰, with ³H 0.five to one TU, confirming they are predominately meteoric waters.

FIGURE 8. (A) distribution of low pmC, HCO₃:Cl, and temperature representing groups I-III. Also shown are the extent of the San Pedro-Linzor low resistivity zone (Mancini et al., 2019), and the El Tatio hydrothermal field (Giggenbach, 1978) and low resistivity zone (Letelier et al., 2021), the Cerro Pabellon geothermal field (Godoy et al., 2022), and the Sol de Manana hydrothermal field (Morata et al., 2020). (B) three component mixing diagram for groups I-III. See Figures 5, 6 for key.

Group II waters are exclusively from the Salado subcatchment, originating in the El Tatio hydrothermal system where HCO₃:Cl ∼0.01, SiO₂ > 100 mg/L, ³H <1 TU and d-excess <0‰ (locality 22). Río Salado samples show decreasing values of ¹⁴C_DIC from 50 to 10 pmC downstream, whilst also decreasing in Cl from 7,760 to 1820 mg/L along the flow path. Moderate pmC, together with high Cl and SiO₂ at locality 22 is interpreted as a mix of carbon dead thermal water (high Cl, SiO₂) (similar or close to a parental magmatic water), with meteoric water that has been in contact with CO₂ in the soil zone. Downstream (localities 17, 21), the fall in pmC and both Cl and SiO₂ suggests the dilution of thermal water with old groundwater discharge.

Group III waters (localities 11, 14, T1, T3) have not previously been identified as a discrete source. They are distinct, having the lowest ¹⁴C_DIC (<10 pmC) and ³H (<1 TU), coupled with high HCO₃:Cl (∼1) and low SiO₂ of ∼70–90 mg/L. These waters are concentrated either side of the Pleistocene San Pedro-Paniri-Toconce volcanic axis including the Turi Basin springs, so that it is hypothesized they originate from a discrete but diffuse ascent of Cl-rich magmatic water from the APMB (Altiplano-Puna Magma Body) with significant dilution from meteoric waters. Furthermore, their high HCO₃ may develop from conductively cooled thermal waters that lack any volatile release (Fournier, 1989). It is noteworthy that group III waters are closely associated with Volcan Paniri, La Poruña cinder cone and the Chao dacite dome, that have the most recent volcanic activity in the area (respectively 150ka, <100ka, and 56-110 ka) (Godoy et al., 2017; Godoy et al., 2018; González-Maurel et al., 2019; Marin et al., 2020).

Group IV waters, found in the lowest reaches of the Río Loa (localities 1, 2) have high TDS (8,000–12000 mg/L) as a result of evaporative concentration throughout the river’s course, as well as possible mixing/dissolution from Salar de Llamara brines and/or saline marine aerosol deposits (Figure 9) (Rech et al., 2003). Their ¹⁴C_DIC ≥100 pmC mark them out as flood waters delayed from the previous wet season.

FIGURE 9. Hydrochemical pathways and processes summarized using TDS vs. distance from the ocean (A) and ternary diagrams for specified zones. (B-E) Data based on Tables 3 and 4, as well as Giggenbach (1978) for El Tatio (see Figures 5, 6 for key to symbols; note expanded Ca, Mg, HCO₃, SO₄ in (B–D)).

The acid SO₄ springs that occur at El Tatio and might be considered a separate group, but they represent only ∼5% of the total outflow from the geyser field (Giggenbach, 1978), so it is not surprising that SO₄ is not a major component of the surface water system. Fournier (1989) suggested they may represent vapour dominated waters with little outflow.

5 Discussion

5.1 Baseflow characteristics and interaction with aquifers

A fundamental control on both streamflow and water chemistry is the annual cycle leading to periods of flood, interflow and baseflow, largely driven by summer (DJF) rainfall. Rising river levels in summer are followed by falling levels during autumn and winter (Mar to Nov). How far floodwater penetrates bank storage, the vadose zone or a contiguous aquifer, causing mixing of waters that have variably evolved, depends on a) the stage head difference between the river and the adjacent vadose zone, b) the duration of positive head difference, c) the degree of hydraulic connection and d) the diffusivity of the adjacent vadose zone or aquifer (Welch et al., 2013).

During flood flow (rising/high stage), there is a tendency for rivers to recharge aquifers throughout their length. During interflow (falling/moderate stage) flow reversal (back-drainage from aquifers) begins in the Cordillera and Pre-Cordillera, whilst the river is still leaking to aquifers in the Central Depression. Finally, during full baseflow (falling/low stage), flow continues to be supplied from intra-arc aquifers, whereas in the Pre-Andean basins and Central Depression river flow reverts to leakage into aquifers, albeit with small fluxes. This cycle of flow into and out of aquifers and/or the vadose zone leads to hysteresis, particularly during falling stages which include mixes from long-lived and evolved aquifer storage, indicating that piston flow may not be a good representation of groundwater transit times (Bethke and Johnson, 2002). As demonstrated by both hydrographs and the instantaneous flow transect, almost 90% of the flow in the Río Loa in the Calama Basin originates from the intra-arc aquifers located above 3,000 m. Throughout the annual cycle, intra-arc aquifers maintain the perennial nature of the flow regime.

Hydrogeochemical processes operating within the catchment may be summarized using TDS vs. location within the catchment, together with cation-anion ternary plots for each basin (Figure 9). In the intra-arc aquifers above 3,000 m, groups I, II and III dominate, showing the impact of near surface production of carbonic acid in the recharge zone and silicate weathering at depth in the hydrothermal system. Whereas the Río Loa shows increasing TDS downstream due to evaporative concentration or possibly halite dissolution, the Salado catchment is dominated by hydrothermal Cl that is progressively diluted downstream until the two rivers equilibrate solutes near their confluence. Group III waters appear to be a discrete geothermal source of Cl, located along the Paniri- Toconce-Leon volcanic lineament, and mixed with significant HCO₃-rich recharge water. Shallow groundwater in the Turi Basin is shown to be a mix of groups II and III, together with varying mixes of group I, whereas deep groundwater in the Turi Basin is dominated by group II alone.

In the Calama Basin the principal processes operating are degassing and precipitation of calcite, evaporative concentration and localized contamination of SO₄ and NO₃ in the vicinity of Calama. In the northern sector of the Calama Basin, shallow and deep groundwater as well as shallow groundwater in the southern sector, show strong similarities with surface waters indicating interchange between the river and aquifers. However, the deep southern sector groundwater (which is confined), shows greater affinity with group II than with the rivers. Since the southern lower aquifer groundwater is quite distinct, any hydraulic connection with the upper aquifers or river should be visible as a change in flow rate or chemistry of the latter. However, there are no significant changes either of flow rate or hydrochemistry in the Calama Basin, or in the rivers passage through the Calama gap, so that any outlet from the southern lower aquifer is potentially deep, and/or slow giving time for flows and chemistry to re-equilibrate before discharging.

Passing downstream, evaporative concentration continues across the forearc, and finally, below Quillagua, a sharp increase in TDS may indicate brine leakage from the Salar de Llamara or dissolution of marine aerosol deposits.

Several possibilities exist for the source of water in the intra-arc aquifers: a) direct precipitation infiltration and recharge, b) slab/mantle dehydration, c) deep hydrous crustal melt exsolution and associated metamorphic fluid release, d) near-surface melt decompression and outgassing, e) thermally induced upflow of deeply penetrating meteoric water through buried evaporites or carbonates, and f) far-field recharge. Evaluation of the plausibility of these sources may be approached by considering the variable and invariable components of discharge to rivers from intra-arc aquifers which can be observed in the hydrographs as well as in hydrochemical and isotopic data.

5.2 Variable recharge

The variable component of baseflow changes from year to year and can be correlated with annual rainfall (see Supplementary Material S2), implying that it represents groundwater discharge following direct precipitation recharge. Variable precipitation recharge (group I water) extends over much of the Western Cordillera in suitable conditions (high-intensity, long-duration precipitation coupled with porous, high-infiltration capacity soils) and its intra-arc aquifers, except where overprinted by group II and III waters. (Figure 8).

Baseflow recession analysis (section 3.2) suggests that direct precipitation recharge to intra-arc aquifers above 3,000 m is around 16–23 mm/a (Table 1). This value compares favorably with two alternative and independent methods of estimating direct recharge. The chloride mass balance technique (CMB), estimates of 1–20 mm/a above 3,000 m but is invalid in areas where volcanic Cl is present (Houston, 2006c) A soil moisture budgeting algorithm (SMB), estimates 28 mm/a above 4,000 m (Houston, 2009), and 39 mm/a in the Silala catchment. These recharge amounts give rise to flows of ∼0.09 m³/s in the upper Ríos Loa-Salado today, equivalent to only 9% of the total baseflow (Table 1).

Far-field recharge will be modulated by effective precipitation in its catchment area, thus resulting in variable annual recharge, albeit buffered by distance. Far-field recharge has been suggested to feed the Río Silala above Inacaliri (locality 15) (Peach and Taylor, 2023). This is supported by the low TDS (<250 mg/L) and high HCO3:Cl (>1, as mmol/L) (Aravena et al., 2023). Water balance models (Muñoz et al., 2023; Peach and Taylor, 2023) have been used to suggest a ∼250 km² extension of the groundwater catchment beyond the surface water catchment area. During the period of record (2001–2015) the measured flow in the Río Silala was 0.17 m³/s, never dropping below 0.15 m³/s suggesting the possibility of significant time-invariant flow. Elsewhere, cross-boundary flow has been widely reported (Anderson et al., 2002; Montgomery et al., 2003; Rissmann et al., 2015; Corenthal, 2016; Boutt et al., 2021). However, neither far-field recharge nor winter precipitation is likely to be able to supply intra-arc aquifer baseflow throughout the year as a result of the overwhelming concentration of precipitation in summer (DJF).

5.3 Invariant recharge

The absence of any significant annual or decadal variation in the time-invariant component of baseflow suggests a long-lived source. Above 3,000 m, two groups of time-invariant recharge water can be distinguished; groups II and III (Figure 8B).

Group II waters represent a mixture of meteoric and thermal water. The deeply circulating meteoric component supplies modern ¹⁴C with low Cl, whilst the thermal component is ¹⁴C_DIC free but Cl rich. The thermal component is assumed to arise from crustal sources (Aguilera et al., 2012; Mancini et al., 2019) and/or recycled brine at depth (Risacher et al., 2003).

Group III waters are associated with intra-arc aquifer leakage without significant surface expression. Although these waters appear to be of geothermal origin, as evidenced by ¹⁴C_DIC <10 pmC, their lack of surface expression is hypothesized to be the result of argillic capping, hence potentially restricting degassing and leading to retained (high) HCO₃.

5.4 Potential sources of invariant recharge

Despite the identification of two endpoint water types related to the volcanic arc, their sources remain obscure. However, by combining hydrogeological information with geothermal studies it is possible to speculate on the ultimate origin and evolution of these water types. The recycling of water in convergent margins via the descending slab and associated hydrous melts in the crust, leading to the release of aqueous fluid and vapour near surface suggests that these processes may play a role in generating recharge to intra-arc aquifers (Sparks et al., 2008; Fialko and Pearse, 2012; Ward et al., 2014; Pritchard et al., 2018; Gottsmann et al., 2022). The ascent of H₂O and OH entrained in magma as a fluid, vapour or bound leads to CO₂, Cl, and SO₄ amongst other species, being mobilized (Zellmer et al., 2015).

Ingebritsen and Manning (2002) calculate that average worldwide subduction dehydration could release between 10⁻⁴-10⁻⁵ m³/s/km² aqueous fluids, which, applied to the area of the upper Ríos Loa and Salado would amount to 0.04–0.27 m³/s, seemingly sufficient to contribute substantially to aquifer recharge. Furthermore, Schurr et al. (2003) and Wang et al. (2019) show that rates vary, with some areas, notably the Western Cordillera having more focused fluxes. A thick crustal passage for volatiles might be indicated by the high As concentrations in the Salado-Toconce-San Pedro subcatchments (Tapia et al., 2021), as well as He isotope ratios, interpreted by Hoke et al. (1994) and Barry et al. (2022) to indicate ∼20% mantle and ∼80% crustal volatiles in the discharge from El Tatio and the Western Cordillera. On the other hand, Risacher et al. (2011) argue that the high Cl:Br ratio of thermal waters at El Tatio (average ∼2,500 as mmol/L), precluded any subducted seawater component. Yet the same samples have a mean Na:Cl ratio of 0.86, that of seawater.

Low seismic velocities and resistivities at depths of 4–7 and 15–25 km beneath the APVC (Altiplano-Puna Volcanic Complex) have been interpreted as super-hydrous melts (≤10 wt% H₂O) of the APMB (Chmielowski et al., 1999; Zandt et al., 2003; Ward et al., 2014; Laumonier et al., 2017; Pritchard et al., 2018). Laumonier et al. (2017) estimate the total aqueous mass in the APMB to be 1.4x10¹⁶ kg, which, assuming a fluid density of ±2000 kg/m³ gives a volume of ∼7,000 km³. In similar conditions Hack and Thompson (2011) demonstrate that fluid ascent rates between 10⁻⁴–10⁵ m/s depend on the type of flow (porous vs. fissure) and are controlled by pressure, temperature and melt composition. With such reservoir volumes and ascent rates it would easily be feasible to supply continuous recharge to the intra-arc aquifers.

Geochemical models for the hydrothermal system at El Tatio suggest a deep, hot reservoir whose composition is the result of lower crustal magma fractionation leading to andesitic upper crustal fluids and rock-water interactions, coupled with mixing of shallow meteoric waters recharged at higher elevations (Giggenbach, 1978; Cortecci et al., 2005; Munoz-Saez et al., 2018).

The Paniri-Chao-Leon-Toconce volcanic chain (Figure 1) trends NW-SE displacing the general N-S trend of the volcanic arc. It partially separates the Pre-Andean basins from the intra-arc basins. Forming the Loa-Salado surface watershed to the northeast of the intra-arc aquifers, and parallel to the Paniri-Toconce trend, is the Cordon de Inacaliri-Volcan Linzor trend. These volcanic trends are thought to be controlled by the deep crustal faults associated with the Calama-Olcapato-El Toro and Lipez-Coranzuli lineaments facilitating the ascent of mantle and upper crustal APMB magmas (Godoy and Kojima, 2013; Godoy et al., 2018; Veloso et al., 2020; Hudson et al., 2022). Type III waters are intimately associated with the Paniri-Toconce trend, so it is hypothesized that this pathway, modulated by an extensional crustal regime, has allowed the ascent of aqueous fluids from various levels in the lithosphere.

Low resistivity zones have been located beneath the APVC and notably beneath El Tatio and the San Pedro-Paniri-Toconce trend (Letelier et al., 2021; Montecinos-Cuadros, et al., 2021). These low-resistivity zones occur at depths of <60 m and 200–400 m hinting at possible reservoirs for groups II and III waters. Several hydrothermal, and fumarolic vents associated with argillic alteration are located within and adjacent to the watershed of the Salado (Figure 8) (Scandiffio and Alvarez, 1992; Mancini et al., 2019; Morata et al., 2020; Taussi et al., 2021; Godoy et al., 2022) offering potential pathways for magmatic water to ascend, and mix with surface and groundwater.

The Pastos Grandes Caldera, part of the APMB, and only ∼10 km to the NE of the San Pedro Basin is provisionally discounted as an important source for waters in Río Loa catchment as it is associated with carbonate rich fluids that precipitate significant calcareous deposits (Bougeault et al., 2020; Muller et al., 2020), that are not found in the latter under present day conditions. It may however, be the source of small volumes of high-HCO_3, low-Cl waters found in perched and near surface aquifers in the Silala subcatchment (Aravena et al., 2023).

High surficial heat flows (>100 mW/m²–Ibarra and Prezzi, 2019) coupled with a geothermal gradient up to 250 ^°C/km at Sol de Mañana (Sullcani, 2015) provides a driving mechanism (but not the only) for fluid upflow through the upper crust. Upflow through sediments underling the APVC potentially alter water chemistry. Risacher and Fritz, (1991), Risacher et al. (2003) and Risacher et al. (2011) consider the Cl content of thermal waters to be buried halite and gypsum, but these are recycled salts, and do not explain the ultimate source of Cl for such evaporites. Here, it is suggested that the ultimate source is found deeper in the crust as a result of magmatic-metamorphic exsolution from the APMB as proposed by Hovland et al. (2018a, 2018b) and found by Rissmann et al. (2015) for the Monturaqui-Negrillar-Tilopozo aquifer that feeds the Salar de Atacama to the south of the Calama Basin.

If voluminous limestones underly the APVC cover, they could represent a deep aquifer and a source for carbonates in groundwater and baseflow. Within the upper reaches of the Salado-Toconce, the Cretaceous Lomas Negra Formation outcrops, containing four limestone beds, cumulatively <3 m thick in a 500 m thick redbed sequence (Marinovic and Lahsen, 1984). Conceivably correlating with this Formation (Charrier and Reutter, 1994), small inliers of Late Cretaceous-early Paleogene Yacoraite Formation occur >120 km to the SE of El Tatio in Argentina (Rubiolo et al., 2003), and the equivalent Molina Formation some 85 km NE of El Tatio in Bolivia (Pacheco and Ramirez, 1997). Although these may represent a continuous extension of lacustrine/marine conditions underlying the current volcanic arc, both these formations show facies interpreted as lacustrine shorelines (Rohais et al., 2019; Deschamps et al., 2020; Mutti et al., 2023; Villafañe et al., 2023), which might negate any extension to the Loa catchment. The Na-normalized plot of HCO₃:Ca (Figure 5B) does not indicate any significant carbonate in Salado-Toconce baseflow, other than soil/plant CO₂ generation, nor do the HCO₃:Cl and Ca:Cl ratios suggest the presence of deep carbonates (Figures 5viii, 5x). Calcite saturation is attained below the confluence of the Toconce with the Salado at ∼3,000 m, but there is no suggestion that this is the result of decarbonation of underlying limestones, so sedimentary limestone is provisionally discounted as a carbonate source.

5.5 Potential structural controls on hydrology

Apart from causing abrupt changes in lithology and related hydraulic properties the impact of faults and folds on groundwater flow and related surface water is poorly known. The inherent anisotropy of fault permeability is commonly associated with a low-k core surrounded by a high-k damaged zone resulting in reduced cross-flow and enhanced lateral flow (Bense et al., 2013). In this situation there is likely to be a gradient decrease upstream of the fault and increase downstream. Flow nets for the northern aquifers in the Calama Basin (Jordan et al., 2015; Herrera et al., 2021) show this configuration across the West Fissure. Similarly, upstream gradient decrease and downstream increase occurs across the Aquina fault system (Houston, 2006c) that lies between the Lipez-Coranzuli and Calama-Olcapato-El Toro fault systems with the same orientation. The cause of separation between the northern and southern aquifers in the Calama Basin, identified in the lower aquifer hydrochemistry, is unknown, but may be controlled by the monocline or by the combined control of several crustal faults and lineaments (Figure 2).

In addition to horizontal groundwater flow distortions, surface water ponding, and localized base-level revision (incision), there is every reason to suppose that the major faulting systems that control the location of hydrothermal sources (Veloso et al., 2020) also facilitate the vertical ascent of groups II and III type waters. It must also be assumed that open fissures will enhance flow, so that the state of stress will impact groundwater flow. The upper crustal stress regime in northern Chile is variable in space and time. East-west extension dominates current conditions in the Coast Range and forearc but rotates to NE-SW from the Pre-Cordillera eastwards (Lanza et al., 2013; Norini et al., 2013; Tibaldi et al., 2017). However, it is possible that in the short term, changes in hydraulic properties due to earthquake deformation may play an important role in controlling surface and groundwater flow (Ingebritsen and Manga, 2019). Using GPS data Klotz et al. (2001) suggest deformation is dominated by the earthquake cycle which causes compression, followed by rebound in the direction of plate convergence. Furthermore, Scheihing (2023) demonstrates a substantial increase in transient groundwater flow following a 7.8 M_W earthquake in the forearc, confirming the relevance of seismic activity to surface and groundwater flow.

6 Conclusion

Using an integrated approach to analyze surface water flows, chemistry, and isotopes in the Río Loa catchment during 1999–2000, several storage and flux elements can be recognized. The annual cycle of flood, interflow and baseflow is critical to understanding both flows and chemistry, especially since such elements follow a non-linear cycle over the hydrological year. Baseflow in the Ríos Loa and Salado is overwhelmingly sourced from intra-arc aquifers at >3,000 m, and there is sufficient active storage in these aquifers to maintain perennial baseflow during the long dry season from April to October. These aquifers receive recharge from precipitation infiltration, hydrothermal sources, and it is argued, diffuse upflow from magmatic volatile exsolution and slab dehydration.

Within the Pre-Andean basins there is considerable interaction between the rivers and phreatic aquifers. The complex geology means that aquifer units tend to be somewhat compartmentalized and may be perched, unconfined, semi-confined, or confined and wholly, partially or not-at-all in contact with each other or with surface water and phreatic and perched aquifers as a result of arc volcanism and tectonic activity. Their interactivity with other ground and surface waters is often related to the annual cycle with precipitation and water level/stage acting as controlling factors. Despite river-aquifer exchange there is little evidence for significant additional recharge to the Pre-Andean and forearc aquifers, other than via downstream transfers from upstream sources.

The intra-arc aquifers location, geometry and physico-chemical properties are all inadequately established. Current research tends to be focused on Pre-Andean and forearc aquifers, whereas there is a need for research to be re-focused on these high-altitude intra-arc sources.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author contributions

JH: Conception, Implementation and analysis. Writing draft and final.

Funding

The author declares that funding for the flow gauging, chemical and isotopic sampling and analysis program in 1999 and 2000 was provided by Nazca S.A. The funder was not involved in the study design, collection, analysis, interpretation of data, the writing of this article, or the decision to submit it for publication.

Acknowledgments

Thanks to Linda Godfrey for insightful discussions and comments on a previous version. Also to John Harman for reading and critically commenting on an early draft. The DGA are thanked for their help in obtaining the precipitation samples.

Conflict of interest

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2023.1310088/full#supplementary-material

References