Diverging Water Ages Inferred From Hydrodynamics, Hydrochemical and Isotopic Tracers in a Tropical Andean Volcano-Sedimentary Confined Aquifer System

Hydrogeology in the Andes cordillera reflects its complex geological history. In most cases, groundwater flows through fractures and faults that compartmentalize the volcanic material, and through the primary porosity of the volcano-sedimentary material. The volcanic mineral context and geothermal environment mark the groundwater chemistry, especially in the high concentrations of specific trace elements. This study focuses on the complex system of the Tumbaco – Cumbayá – Los Chillos aquifer, in the vicinity of the Ilaló volcano near Quito (Ecuador). Hydrodynamic, geochemical and isotopic tools were used to assess the chemical characteristics of water and its origin, identify the recharge areas, and estimate the transit time of water using simple methods and scarce data. Results revealed two distinct aquifers, one in the volcanic cone located in the center of the study area, and the other in the volcano-sedimentary series of the Tumbaco – Cumbayá – Los Chillos valley. The volcanic aquifer is characterized by a high mineralization, a recharge zone between 2400 m asl and 3100 m asl, and radiocarbon concentrations lower than 20 pmc. The volcano-sedimentary aquifer seems to behave as a partly disconnected system, between the north and the south of Ilaló volcano, and also with a great heterogeneity, maybe due to the presence of lenses of volcanic ash. It has an intermediate mineralization, a mean recharge zone between 2,300 and 2,700 m asl, and 14C activities between 45.4 and 87.4 pmc in apparent contradiction with the hydrodynamic mass balance.


INTRODUCTION
Volcanic and volcano-sedimentary aquifers often result from a complex geological history inducing a heterogeneous distribution of primary and reworked volcanic materials. Before and after their deposition, topography and tectonics, as well as erosion, influence the nature and geometry of the deposits (Martí et al., 2017;Oyarzún et al., 2017). These different stages impact the groundwater circulation, in various types of rocks with textures ranging from mainly porous to only fractured (Oyarzún et al., 2017;Moreno et al., 2018;Vincent et al., 2019). These geological stages may also modify conditions and fluxes of recharge and discharge (Wannous et al., 2016) either by compartmentalizing the system (Moreno et al., 2018) and forming water lenses (Wannous et al., 2016) or by interconnecting aquifers (Gonzales et al., 2019). The physical and chemical groundwater characteristics are influenced by the dissolution kinetics of minerals of the rock matrices, which can be enhanced by high temperature generated by regional volcanism or local geothermal flow (Wrage et al., 2017), and by interaction with deep volcanic gases (Federico et al., 2017;Mukherjee et al., 2019).
Among the dissolved elements present in volcanic waters, some may constitute a threat to public health such as uranium, selenium, vanadium or arsenic (Bundshuh et al., 2010;Cinti et al., 2015;Wrage et al., 2017;Godfrey et al., 2019;Mukherjee et al., 2019). Several studies in South America examined the arsenic geochemical behavior and its potential sources, especially in Chile, Argentina, Bolivia, Peru and Ecuador (López et al., 2012;Tapia et al., 2019). Cumbal et al. (2010) documented arsenic concentrations in numerous springs and surface waters and investigated the evolution of As concentrations in rivers and sediments in Ecuador. They found a large variability in As concentrations (2 to 969 µg l −1 ) in geothermal spring waters within 150 km from Quito.
Facing the natural complexity of geological structures and groundwater circulation, hydrogeologists require accurate and dense information. In contrast, field reality in the Andes is frequently scarcity and questioned reliability of data (Comte et al., 2012;Gómez et al., 2019). This implies to develop a modest approach, that gathers all previous information in a large range of scientific disciplines, reanalyse it, supplement it with new data and propose explanations that are probably not definitive, like Zhu (2000), Phillips and Castro (2003) and Awaleh et al. (2017) did in volcanic areas of other continents.
In Ecuador, where volcanism and tectonism are still very active, the few studies focusing on hydrogeology combine most often different approaches. In the Galapagos islands, Auken et al. (2009) and Pryet et al. (2012) used geology and geophysics for a better characterization of the volcanic aquifer, as did Zarroca et al. (2015) in the South tropical Andes. Complementing the field information with a numerical modeling is an efficient way to explore the various hypotheses open [e.g. Guzmán et al. (2016), Carrión-Mero et al. (2021), Chidichimo et al. (2018)].
Our study focuses on a small regional aquifer system of the North Tropical Andes near Quito, the capital of Ecuador. The metropolitan city of Quito, faced with surface water contamination, industrial development and doubling of population in the last decades, needs to develop its water supply system, and looks into groundwater to meet its needs. The Tumbaco-Cumbayá-Los Chillos aquifer is a potential resource in this prospect. Because of the abundance of surface waters, few hydrogeological studies were carried out in the area. As in many volcanic parts of the world, data is scarce and/or of poor quality, which significantly limits the knowledge of the geometry, hydraulic connections and storage capacity of this aquifer system, as well as groundwater quality, especially regarding arsenic levels (Cumbal et al., 2010).
This study aims at understanding the hydrodynamic behavior of this aquifer system, the connection between the different layers and the processes controlling the groundwater mineralization. Some specific issues such as the origin of arsenic, areas of recharge and the groundwater transit times will be detailed, as far as allowed by the present state of knowledge.

Geography of the Study Area
The study area is located roughly 10 km to the south east of Quito, in the center of the Interandean Valley, near the equatorial line ( Figure 1A), and includes the Tumbaco-Cumbayá and Los Chillos valleys ( Figure 1B). It covers 700 km 2 (about 50 × 19 km, oriented SW -NE) and has the following limits: to the north, the junction between the San Pedro and the Chiche river canyons; to the south, the Pasochoa volcano; to the east, the "Cordillera Real"; and to the west, the San Pedro river canyon. Urban structures and agricultural activities are important, especially in the northern part. The most important city in this zone is Tumbaco and the total population of the study area was about 322,000 inhabitants in 2010.
The study area topography is essentially mountainous, with altitudes ranging from 2,200 to 4,200 m asl (average of 2,400 m asl). In its center, the Ilaló volcano (3,169 m asl) separates the Los Chillos valley in the south from the Tumbaco -Cumbayá valley in the north. The Pasochoa volcano (4,200 m asl) at the southern limit is the highest point of the study area ( Figure 1B). The most important water courses are the rivers Pita, San Pedro, Machángara and Chiche ( Figure 1B). Erosion by the rivers San Pedro and Chiche has cut deep canyons, from 200 m deep in the study area to 600 m further to the north.
Average monthly precipitation for the period 1963-2019 was calculated using data from six stations operated by INAMHI (Instituto Nacional de Meteorología e Hidrología) in the study area and in the Quito metropolitan area. The Inter Tropical Convergence Zone latitudinal shift causes two annual wellmarked wet seasons: the first maximum is observed from March to May (50% of the annual precipitation) and the second one from October to December (30%). Precipitation is higher in the south (1,700 mm.a −1 ) and lower in the north (780 mm.a −1 ). Average monthly temperature from 17 INAMHI stations shows locally a relative stability all year long. It is generally colder in the south (15.5 • C, average altitude of 2,500 m asl), and warmer in the north (16.7 • C, average altitude 2,300 m asl), with a temperature gradient of −0.65 • C/100 m. The mean monthly relative humidity varies between 75% (September) and 86% (February). The annual potential evapotranspiration using the Penman-Monteith method (Allen et al., 1998) is about 800 mm.a −1 (Manciati, 2014).

Geological Context
The geological map of reference for the study area was established by the Dirección General de Geología y Minas (DGGM, 1980). The sedimentary deposits distribution in this map has been challenged by the results of K-Ar dating of the volcanic Ilaló formation made by Barberi et al. (1988). The other formations were never dated radiometrically, but only based on the stratigraphy. The geological formations ( Figure 1C) are subsequently presented from bottom to top.
The Volcánicos Ilaló formation (Lower Pleistocene) is observed near the Ilaló volcanic cone and dated about 1.62 ± 0.16 Ma (Barberi et al., 1988). This confined aquifer is mainly made of breccias and lava. Its lateral extension and thickness are unknown. Its hydrothermal waters are exploited for recreational activities and, formerly, for human consumption. Several thermal springs are present around the Ilaló volcano.
The Guayllabamba volcano-sedimentary formation (Lower -Middle Pleistocene) is observed at depths of 170 to 200 m in the northern part of the study area and outcrops at the bottom of the Chiche River canyon. It extends to the north beyond the limits of the study area. Its local lithology is made of lava flow and pyroclastic sediments, about 200 m thick. This formation is considered as an aquifer, exploited by the Guayllabamba people 20 km to the north, but there is no well-exploiting this layer in the study area.
The Chiche volcano-sedimentary formation (Middle -Upper Pleistocene) extends over the entire study area, and further to the north out of the study area. In the north it is entirely cut up by the Chiche and San Pedro rivers. Its lithology is extremely heterogeneous: it is made of conglomerates, thick sand, some compacted volcanic ash layers and tuffs well-stratified. In the study area its thickness varies from 230 m to the south of the Ilaló volcano to 70 m in the north. This formation is the main focus of the present research. Its water was distributed to the Tumbaco-Cumbayá-Los Chillos population for domestic consumption but is no longer exploited since 2006 because of high arsenic concentrations (mean content between 50 and 60 µg.l −1 ). The hydraulic continuity of the Chiche formation north and south of the Ilaló volcano is questioned due to the presence of the volcanic cone which greatly separates the north and south parts of the sedimentary volcano series. The vertical connectivity within the Chiche formation is also an important characteristic to assess, as the volcanic ash layers could behave as vertical hydraulic barriers. However, the lack of finely resolved stratigraphic logs complicates the resolution of this question. The hypothesis of a multilayer system will be discussed in the later sections.
The Cangahua formation (Upper Pleistocene to Holocene) covers a large part of the study area, and is present over the whole Interandean Valley. The typical lithology consists of altered tuffs, often intercalated with volcanic ash and pumice with a silty texture. It is thicker in valleys and depressions (about 50 m) than in upper zones (20 -30 m). Considered by Quantin and Zebrowsky (1996) as an impervious top layer of the Chiche formation, it could prevent direct diffuse recharge to the Chiche aquifer.
In this tectonically active zone, preferential directions of faults are NE-SW and NW-SE. Major faults are (i) to the west side, the Quito fault and the Machángara River fault; (ii) the San Pedro River fault joining the Chiche River fault (Figures 1B,C). All these structures are linked to different vertical movements upward and downward that command geomorphology and hydrogeology (Rios-Sanchez et al., 2012). The Ilaló and Pasochoa volcanoes are part of an aligned volcanic complex. Their slopes are cut by deep radial canyons, and their calderas collapsed westwards. Their imperfect cones are considered as typical of old volcanic structures (Cornejo, 1983).

MATERIALS AND METHODS
The pre-existing geological information comes mainly from unpublished technical reports of the DGGM, and academic studies about geology and hydrogeology of surrounding aquifers. This information was compared with lithological and technical sections of wells.
In the study area, the Ilaló aquifer and the Chiche aquifer are reached by 147 registered wells, of which 107 are private and not accessible. The other wells were used for groundwater observation. Information about the construction profiles (technical logs) of 24 wells is available. Among them, two are in the volcano piedmont, they screened in the Ilaló formation. In the other wells, the screens are open to all the groundwater horizons without distinction. In 19 wells, the groundwater level was measured monthly: 12 wells in Tumbaco -Cumbayá (September 2009-October 2012, total of 372 measurements), and seven wells in Los Chillos (February 2011-October 2012, 109 measurements). Among them, three have been equipped with automatic level recorders since September 2010. A differential GPS land survey allowed to determine their elevation with an accuracy of 1 cm. This information was used to draw a piezometric map (Figures 2A,B), extrapolating contour lines by a triangulation method with the nearest neighbors. All Ilaló wells are under pressure, and some are artesian but no manometric measurements are available. The results from 11 well-pumping tests performed at the time of the well-construction, and lasting over 24 h, were used to estimate transmissivity and permeability. All hydrodynamic data were evaluated with basic statistical parameters (maximum, minimum, average values and standard deviation calculations) allowing to detect possible measurement errors, in which rare cases the doubtful data were removed.
A large part of the chemical dataset used in the present study was provided by the Water Utility Company of Quito (Empresa Pública Metropolitana de Agua Potable y Saneamiento -EPMAPS) which has been surveying 24 wells for major ions and trace elements every 6 months since 2003. Among the chemical analyses (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011) available from the EPMAPS laboratory, some were discarded because of inadequate ionic balance (greater than 20%), or when the major ions analyses were incomplete. A selective subset of 100 analyses was considered in the present study. Physical parameters such as electrical conductivity (EC, error ± 0.5%), TDS (± 0.2 mg.l −1 ), pH (± 0.01), temperature (± 0.1 • C) and occasionally redox potential Eh (± 0.3 mV) were measured in the field. Samples dedicated to chemical analysis were filtered at 0.45 µm. Alkalinity and chloride were measured in the laboratory by titration (± 5%). Sulfate, nitrate and dissolved silica were measured by spectrophotometry (±5%). Phosphate and fluoride were determined by colorimetry (± 7%). Cations Ca 2+ and Mg 2+ were measured by complexometry (±5% and ±10%, respectively), and Na + and K + cations by flame atomic absorption (± 0.03 mg.l −1 ). Dissolved oxygen was determined in the laboratory with the Winkler method. Minor and trace elements including Al, Cd, Co, Cu, Cr, Fe, Li, Mn, Ni, Zn were measured using flame atomic absorption spectrophotometry and arsenic (As) was obtained by atomic absorption spectrophotometry using the hydride generation method (error ±0.03 mg.l −1 ).
Three additional geochemistry campaigns were performed by the Institut de Recherche pour le Développement (IRD) in April 2010, July 2010 and October 2012. Among the 11 samples collected, 5 originated from the Chiche North aquifer, 2 from wells presumed to mix waters from the Chiche and Ilaló aquifers, and 4 (1 well and 3 springs) from the Ilaló aquifer.
They were analyzed at the HSM (HydroSciences Montpellier) water chemistry laboratory, University of Montpellier, France. A few samples were duplicated for an inter-comparison exercise with the EPMAPS laboratory. For these samples, besides major and trace element analysis, the speciation of arsenic (reduced As(III) and oxidized As(V) species) and iron (total dissolved Fe and Fe(II)) was also conducted to gain information on the processes controlling As release into groundwater. In the field, these samples were filtered at 0.22 µm for cations, anions, total dissolved arsenic, arsenic speciation, total dissolved iron, iron II and trace elements. Samples were preserved with ad hoc reagents to avoid chemical or microbial evolution during transport and storage (Gallagher et al., 2004). Samples for cations and trace elements were added one or three drops of suprapur nitric acid in 10 mL and 60 mL plastic bottle, respectively. For arsenic speciation, samples were taken in a 10 mL plastic bottle preserved with 100 µL of 8,7 M acetic acid and 100 µL of an EDTA 50 g/L solution. For Fe(II), the sample was taken in a 10 mL plastic bottle and preserved using a 500 µL buffer of ammonium acetate and 1 mL of 0.5% in mass of phenanthroline chloride. All samples were preserved in a +4 • C cooler and covered with parafilm. Samples for Fe II were also protected with aluminum foil to avoid photochemical oxidation. Dissolved oxygen was measured in the field with the colorimetric method (error ± 0.02 mg.l −1 ). For this series of samples, ionic balance was lower than 7%. Alkalinity was measured in the field by titration (± 0.1 mg/l). Major ions were measured by ion chromatography (± 0.02 mg.l −1 ). Total Fe and As contents were determined by ICP-MS (Fe detection limit = 0.3 µg.l −1 , As detection limit = 0.03 µg.l −1 , uncertainty ± 1%). Fe II contents were obtained by spectrophotometric method (wave length λ = 510 nm, detection limit = 88 µg.l −1 , uncertainty ±5%). Arsenic speciation was determined using high performance liquid chromatography (HPLC)-ICP-MS (detection limit = 0.2 µg.l −1 for As(III) and 0.4 µg.l −1 for As(V), uncertainty ±5%).
Between January 2010 and September 2012, four rain gauges (Rumipamba INHAMI station on the east flank of Pasochoa volcano, La Tola INAMHI station in the piedmont of Ilaló volcano, Cumbayá and Las Peñas stations, see location in Figure 1B) were run at elevations between 2,360 and 2,970 m asl, for a monthly isotopic survey of rainwater in Tumbaco -Cumbayá and Los Chillos valleys, amounting to a total of 112 samples. Data from Quito INAMHI station were also available for this study since 2001 (see Figure 1B for location). Aquifers were sampled through 40 wells: 23 Chiche North and South, 13 Ilaló, and 4 wells under pressure located within 2 km from the volcano, where the borehole water may result from a mixing of North Chiche and Ilalo aquifers, and hereafter referred to as mixed wells. The sampling was conducted during wet and dry seasons for 18 O and 2 H measurements, representing a total of 72 samples collected in the months of April and August of years 2010 and 2011. Water isotopes were analyzed at the LAMA laboratory of HSM, on an Isoprime mass spectrometer using the CO 2 equilibration and Dual Inlet methodology for δ 18 O measurements (error ± 0.06 ‰), and water reduction over chromium at 1050 • C using a Eurovector PyrOH analyzer for continuous-flow δ 2 H measurements (error ± 0.8‰). Water stable isotope compositions are expressed on the VSMOW scale. Samples for radioactive isotopes were taken in April 2010, July 2011 and October 2012 in 12 wells (7 North Chiche, 2 mixed and 3 Ilaló): 8 samples for 3 H (error ± 0.5 TU) and 12 samples for 14 C and 13 C (errors ± 0.5 pmc, and ± 0.05‰ vs. VPDB, respectively), all analyzed at the University of Avignon, France. All samples for stable and radioactive isotopes were taken following the protocols described in Clark and Fritz (1997). Our results were compared with those of Changkuon et al. (1989) obtained in Los Chillos and Tumbaco -Cumbayá valleys (chemical and isotopic data in 41 wells).

Water Budget Calculations
An approximate water budget was calculated to estimate the infiltration and flow velocity for the North Chiche aquifer (Figure 3) and to infer water transit time from a hydrodynamic point of view, to be compared with radioactive isotopes dating and to assess the heterogeneity of the system. Due to the limited number of field measurements and the large variabilities observed, the hydrodynamic calculations presented in this section give orders of magnitude only. The North Chiche aquifer was selected because of the denser hydrodynamic information available, and the more accurate geology profile. The main uncertainties to complete this rough estimate concern the mean hydraulic conductivity and porosity distribution, and the possible multilayered structure of the aquifer. To present all these analyses, the Chiche aquifer will be treated as a single unit.
Infiltration was calculated from a theoretical discharge flow rate using Darcy law (equation 1), respecting the continuity law in steady state.
where e m is the saturated thickness (average value of 100 m); n is the porosity. This range of values leads to an estimate of V stock between 70 and 700 hm 3 . Using this information, we calculated the water balance in the North Chiche aquifer. Discharge flow rate was estimated from Equation 1 at 0.09 m 3 .s −1 (Q max = 0.35 m 3 .s −1 , Q min = 0.04 m 3 .s −1 ) or 3 hm 3 .a −1 (Q max = 11 hm 3 .a −1 , Q min = 1.3 hm 3 .a −1 ), not taking into account water output from pumping, because piezometric levels do not show significant variations during the year. An example of high-frequency groundwater level variations in a well (Patagua -C105, see location in Figure 1B) exploiting the North Chiche aquifer is shown in Figure 4. In spite of visible transient reactions to pumping, static groundwater levels in the North Chiche aquifer varied slightly: < 0.05 m in a month, and 0.40 m in a year (Figure 4).
If Cangahua is considered an impermeable top layer in the valleys, infiltration may take place in the Ilaló piedmont, between 2700 and 2400 m a.s.l., in an 18 km 2 area (Figure 3). If we consider a diffuse infiltration, the net infiltration has to be 160 mm.a −1 (600 mm.a −1 , 70 mm.a −1 ), i.e. between 75% and 10% of the total precipitation, to reach the calculated flow rate. These calculations using this simple method will help characterize the transit time of the Chiche aquifer.

Physico-Chemical Characterization of the Aquifers
The geochemical tools presented in this section allow to define the conditions of the chemical balance between the groundwater and its reservoir along its course. The spatial distribution of EC values (Manciati, 2014) suggests three distinct compartments. The wells intercepting the South Chiche aquifer had EC between 300 and 400 µS.cm −1 , while those in the North Chiche aquifer showed EC values between 500 and 700 µS.cm −1 (Supplementary Table 1). Wells and springs around Ilaló volcano, presumably intercepting the Ilaló aquifer, had higher EC, from 1,000 to 3,000 µS.cm −1 . These values are consistent with the general groundwater flow directions inferred from the piezometric data of the Chiche aquifer (Figures 2A,B).
Surface groundwater temperature was between 17 and 19 • C in South Chiche wells, and between 19 and 22 • C in North Chiche, higher than the local mean annual air temperature of about 16.7 • C. Wells connected to the Ilaló aquifer and those supposed to mix Chiche and Ilaló have notably higher temperatures, between 27 and 42 • C (Supplementary Table 1). In the Chiche aquifer, North and South, pH varied between 7.0 and 8.5, while it ranged between 6.2 and 8.1 in the Ilaló aquifer (Supplementary Table 1). We observed in Ilaló wells an excess of CO 2 , due to emanating gases from the Ilaló volcano passing through faults. Redox potential (Eh) was typical of a transitional state, for both Chiche and Ilaló aquifers: between 70 mV and 190 mV for Chiche and between 170 mV and 225 mV for Ilaló during field measurements conducted from February 2011 to October 2012. Supplementary Table 1 shows some of the values during wet and dry seasons for wells intercepting the two aquifers.
Saturation indexes (SI) were calculated as the logarithm of the ratio between the ion activity product and the equilibrium constant K of the corresponding mineral: log (SI) = log ion activity product K = log ion activity product −log (K) Minerals containing sodium, magnesium, manganese combined with sulfate, chloride or fluoride are under-saturated. Some minerals including calcium carbonate, iron carbonate and manganese carbonate are in equilibrium. Water is oversaturated in calcium magnesium minerals (dolomite), in silicates as quartz, talc and chalcedony, and in a ferrous mineral, goethite, which could explain precipitation of iron in wells. Trace element analyses by ICP-MS showed that Chiche aquifer has a mean arsenic concentration varying between 1 µg.l −1 and 30 µg.l −1 . The South Chiche aquifer had the lowest concentrations. The Ilaló aquifer mean arsenic concentrations varied between 50 µg.l −1 and 60 µg.l −1 , with a maximum value of 600 µg.l −1 in the La Merced well where water temperature is 37 • C. Among the 143 samples from the North Chiche and Ilaló aquifers, 90% of wells and 70% of springs have arsenic concentrations higher than the NTE INEN 1108:2011 Ecuadorian norm (INEN, 2011) of 10 µg.l −1 . Interconnection between Ilaló and North Chiche aquifers is strongly suggested by arsenic concentrations.

Water Stable Isotopes
Rainwater isotopic composition from the four stations operated in the study area showed a clear seasonal signal following that of monthly rain amounts, with wetter months leading to a more depleted isotopic signal, consistently with the amount effect whereby a more intense rainfall causes a higher depletion of the remaining atmospheric humidity (Figure 7). The four stations The weighted mean isotopic composition measured for the four stations during <3 years shows a clear correlation with the elevation of the stations, with a δ 18 O gradient of −0.16‰ every 100 m of elevation (R² = 1.00). The robustness of this correlation is somewhat limited by the short size of the data set. Rumipamba station, which showed the most depleted isotopic signal, is the highest (2,970 m asl), ∼600 m above the other rain gauges (from 2,340 to 2,480 m asl). The correlation between altitude and isotope is mostly driven by this station. The higher amount of annual rainfall at this location, the distance between this station and the northern ones, and the orography in between, prevent from assigning the isotopic difference solely to the altitude effect. Previously, the altitudinal gradient was calculated by Aranyossy et al. (1989) in the Quito area and Los Chillos valley, with a value of −0,33‰ per 100 m in δ 18 O. This value, based on the measurements of nine springs located at different elevations, and consistent with our data given the uncertainties, is used later in this study to assess the mean recharge altitudes of the aquifers.
Likewise, Changkuon et al. (1989) also measured some variability within the Ilaló and Chiche aquifers, and observed that the Ilaló aquifer was more depleted in heavy isotopes than South Chiche in the Los Chillos valley.
In the δD vs. δ 18 O plot (Figure 8), the Ilaló and South Chiche aquifer isotopic compositions are distributed close to the Global Meteoric Water Line (GMWL) as well as the local meteoric water line. Conversely, the North Chiche data points diverged notably from the GMWL, with a slope (δD/δ 18 O) of 3.5, that will be discussed in the later section.

Radioactive Isotopes
Tritium was measured in 8 samples of water from Chiche and Ilaló aquifers. Groundwater from the two aquifers did not have detectable tritium. Only one spring, located at the north and fed by the Chiche aquifer (Chirimoyas 1 -C201 - Figure 1B) showed a measurable concentration (0.6 TU), indicating recent water ( Table 1). This result was confirmed by the radiocarbon analysis for this spring with a 14 C activity of 103.3 pmc. Radiocarbon analyses of all other samples revealed relatively old groundwaters for both Chiche and Ilaló aquifers. North (Table 1) and South Chiche (Changkuon et al., 1989) aquifers had 14 C values between 45,4 and 87.4 pmc, which correspond to an apparent age between 1000 and 6000 years. Groundwater from Ilaló aquifer showed 14 C activities lower than 20 pmc, giving an apparent age between 11,000 and 46,000 years. These ages have to be corrected using the δ 13 C values which allow to discriminate the biological carbon sources from, in this case, the deep geogenic volcanic CO 2 which can disturb calculated water ages by adding 14 C-free dissolved inorganic carbon (DIC), making water appear older than it is in reality (Vogel and Ehhalt, 1963).
We used a simple mixing model proposed by Salem et al. (1980) to correct apparent ages estimated from 14 C activities (Han and Plummer, 2016), taking into account the C coming from the volcanic source, using the following equations: where t is the age calculated from the measured 14 C activity (A), and A 0 is the initial 14 C activity of the total dissolved carbon mixture at the time of isolation of the recharge water from the soil gas phase. A 0 is calculated from equation (6), where δ 13 C DIC is the dissolved inorganic carbon δ 13 C value measured in groundwater, δ 13 C vol is the δ 13 C value of volcanic CO 2 , δ 13 C bio is δ 13 C value of soil CO 2 , and δ CO2−DIC is the isotopic enrichment factor between soil CO 2 and TDIC (total dissolved inorganic carbon, in our case mainly HCO − 3 ) at the mean annual air temperature.
The north spring (Chirimoya 1 -C201), characterized with modern water, was taken as the carbon 13 biogenic reference for the area, and the initial 100% in radiocarbon during the recharge, its TDIC being only formed by exchange between soil CO 2 and the liquid phase. This δ 13 C value is − 10.62‰ vs. VPDB (Table 1). Taking into account the fractionation factor of − 8.8‰ corresponding to 16.7 • C (mean annual air temperature), the theoretical CO 2 of biological origin might bear a δ 13 C value of ∼-19,4‰ vs. VPDB. This 13 C value is consistent with the present C 3 type vegetation, which is dominant in the study area, as grass (Lollium perenne and Calamagrostis effuse), eucalyptus and legumes (beans, potatoes, and broccoli), mixed to a limited extent with CAM and C 4 type vegetation, as cacti, corn or the grass Pennisetum clandestinum (kikuyu), also observed in the study area. According to Ledru et al. (2013), vegetation cover did FIGURE 8 | δD and δ 18 O in rainwater and groundwater. Circles represent weighted mean annual rainfall isotopic composition in the study area. Squares and triangles represent groundwater isotopes from this study and from Changkuon et al. (1989). Global meteoric water line (GMWL) and regression lines for groundwater isotopes from the three aquifers are plotted. not change significantly since the first human settlements in the valley, 3000 years ago.

Aquifers Recharge
Stable isotopes of water were principally used to estimate the mean altitude of the recharge area. As expected, the highest station is more isotopically depleted than the three others, but the altitudinal range is too limited to define properly a local isotopic altitudinal gradient. Moreover, the particular climatic conditions of 2011 and 2012, with two La Niña episodes following a strong El Niño and their unusual isotopic depletion in relation with the mass effect (Figures 7, 8), may have somewhat obscured the expected altitude effect (Hurley et al., 2019). Consequently, we used the longer time series of Quito (δ 18 O interannual mean of −12.18‰) as a regional reference, together with the altitudinal gradient 0,33‰ per 100 m in oxygen18 defined by Aranyossy et al. (1989).
The mean recharge altitude for the Chiche aquifer, was thus estimated between 2,300 and 2,700 m asl, assuming climatic conditions at the time of recharge to be similar to present ones, consistently with several late Holocene environmental studies carried out at similar altitudes (Van der Hammen, 1974;Servant and Servant-Vildary, 2003;Ledru et al., 2013). This reinforces the assumption of an impervious Cangahua formation overlying the Chiche and infiltration could predominantly take place in the volcanoes piedmonts (Ilaló in the north, Pasochoa in the south), where deep gullies may cut the Cangahua formation. A similar calculation for the more depleted Ilaló aquifer indicates a recharge elevation range between 2,400 and 3,400 m asl, the upper limit being higher than the actual present-day elevation of the volcano (3,100 m asl). Current local infiltration may be limited by the volcano morphology where slopes are very steep. Infiltration may take place in the gullies localized in the western part of the Ilaló volcano, where Cangahua formation is not present. Infiltration could also have taken place before the Cangahua deposition, during the Upper Pleistocene when air temperature was colder and weather was drier (Novello et al., 2017). In this case the precipitation should have been isotopically more depleted than present, and would correspond to a direct infiltration on the slope of Ilaló volcano. This would imply very low renewal rates in the North and South Chiche and Ilaló aquifers, making, specially the Ilaló aquifer, a non-renewable resource. This hypothesis is consistent with radiocarbon results Corrected Ages = ages calculated using Salem et al., 1980 model. but does not agree with the hydrodynamic interpretation. This apparent discrepancy will be discussed below.

Hydrodynamics and Connectivity Between Aquifers
The water table contour map implies the general direction of streamlines from south to north (Figures 2A,B). The question whether Ilaló volcano represents a hydraulic barrier between South and North Chiche is difficult to debate. Geologically, the two parts of the Chiche aquifer are in the same formation, which may explain the similarities in their chemical properties, but the volcano constituted a physical barrier during the sedimentary deposition post lower Pleistocene, at least partially. There are no existing wells on the sides of Ilaló where to test this connection by field work. It is highly probable that Ilaló limits at least substantially the flow between Chiche South and North. Levels of groundwater and rivers were compared (Manciati, 2014). In the northern part, where deep river canyons cut entirely the Cangahua and Chiche formations, the North Chiche water table was always higher than rivers by at least 30 m even in the dry season. In the southern part of the Chiche aquifer, the relationship between groundwater level and the rivers inside the basin was not so clear. However, the piezometric contour maps confirm that River Chiche in the north (Figure 2A), and River San Pedro in the south (Figure 2B) are major drainage axes. Comparison between groundwater levels in the Ilaló formation and rivers is not meaningful, considering local and regional gradients and the distance between wells and rivers.
Groundwater is generally between 20 and 50 m below the ground surface and the steeper hydraulic gradients (2% to 5%) are close to valleys and springs. Wells reaching the Ilaló formation, the Chiche formation in the south (Los Chillos valley), and wells mixing Chiche and Ilaló are under pressure; some are even artesian (Figures 2A,B).
Using the information of water balance in the North Chiche aquifer, discharge flow rate was estimated as 0.09 m 3 .s −1 (Q max = 0.35 m 3 .s −1 , Q min = 0.04 m 3 .s −1 ) or 3 hm 3 .a −1 (Q max = 11 hm 3 .a −1 , Q min = 1.3 hm 3 .a −1 ). According to this calculation, the annual flow rate is much higher than the flow extracted from the current pumping (<7% of average annual flow rate). The distance between the recharge area and the outlet area for North Chiche being about 6400 m (Figure 3), transit time of groundwater may be estimated, using the range of acceptable values, of the order of a few decades to centuries depending on the mean porosity considered (between 2% and 20%). This estimate agrees well with the comparison of the aquifer volume (V stock , between 70 and 700 hm 3 ) and the calculated annual input flow (3 hm 3 .a −1 ). The lower limit of 5 years appears obviously to be much too short. Effective porosity is probably much lower than the upper bound of 20% we used in this preliminary calculation. As explained in the geological settings section, North Chiche aquifer is rather heterogeneous and the compacted volcanic ash layers could form lenses and may reduce both porosity and vertical permeability. These results will be compared and confronted with the radioactive isotopes measurements hereafter.

Mineralization and Geothermalism
Wells intercepting Ilaló and North Chiche aquifers, or a mix of both, show water temperatures higher than average air temperature, suggesting a potential geothermal effect presumably related to the activity of Ilaló volcano (Supplementary Table 1). This influence can be observed up to 5 km northward from Ilaló volcano. In the Ilaló aquifer the presence of sodium in water may be explained by an average content of 53% of andesine in geological formations, in mineral matrix and in the phenocrysts (andesine is at 60% a sodium plagioclase, albite NaAlSi 3 O 8 ). At a higher water temperature, the dissolution process is enhanced, possibly increasing the sodium content in groundwater. The excess of sodium in water found in the north Chiche aquifer may then not come from the local mineralogy, but from a deep flux through the Ilaló volcano faults. The prevalence of Mg in the south Chiche aquifer could result from interactions of water with iron-magnesium bearing rocks like basalts and andesites with olivine, pyroxene (augite -(Ca,Mg,Fe) 2 (Si,Al) 2 O 6 , hypersthene -(Mg,Fe)SiO 3 ), hornblende which are present in the geological formation.
In the Korjinski alumino-silicate activity diagram at 25 • C (Figure 9A), data from South Chiche, North Chiche and Ilaló aquifers have been plotted. Data points from the three compartments show a passageway between kaolinite and pyrophillite. Theoretically, if the equilibrium is reached in the aquifer system, then pyrophillite is not stable, so the final product of aluminosilicate activity is kaolinite (Figure 9B).
Arsenic is usually found mainly in the form of arsenopyrite (FeAsS), a mineral formed in high temperature hydrothermal systems (≥250 • C) associated with pyrite (FeS 2 ), pyrrhotite (Fe (1−x) S (x = 0-0.2)) or cassiterite (SnO 2 ). Arsenic may also be present in hydrothermal systems as low temperature hydrothermal minerals (≤150-200 • C) orpiment (As 2 S 3 ) and realgar (AsS) (Pokrovski et al., 1996). Several studies in the Andes showed arsenic to be contained in volcanic glass, often present in geological formations (Peters, 2008). This could be the case for Chiche and Ilaló aquifers, volcanic glass being part of the matrix of these formations (Smedley et al., 2005;Coomar et al., 2019). Positive correlations were found with B and Fe in volcanic regions around the world López et al., 2012) suggesting a common origin of these elements. In the North Chiche aquifer a trend has been found between arsenic and boron (Figure 10A), although boron was not measured systematically in the aquifer system so there is insufficient data to support the significance of the correlation. A principal components analysis was performed on the chemical dataset, and is presented in Figure 10B. The first three components represent 91% of the variance. In the Ilaló aquifer, a possible correlation was found between arsenic and iron (R 2 = 50%, not shown). This correlation suggests the association of As with iron sulfides or iron oxides (López et al., 2012). This last hypothesis is consistent with mineralogical analyses of the Ilaló formation, where magnetite (Fe 3 O 4 ) was observed in the matrix and phenocrysts, with hematite (Fe 2 O 3 ) as a secondary mineral. These secondary minerals may indicate alteration processes occurring at high temperatures. The isotopic data of Ilaló groundwater did not reveal a geothermal exchange in 18 O, contrary to the North Chiche groundwater, which stable isotope data points diverge from the GMWL, with a slope (δD/δ 18 O) of 3.5 (Figure 8). This slope is inconsistent with an evaporation effect, but it may trace either direct geothermal water-rock exchange processes at relatively low temperatures, or a mixing line representing aquifer water of meteoric origin contaminated to some extent with geothermal water. Indeed, a specific 18 O enrichment may result from isotopic exchange with isotopically enriched oxygen-bearing minerals when water comes in contact with rocks at high temperature (Tassi et al., 2016;Ma et al., 2017;Cabassi et al., 2019). Further research is needed to understand the circulation of the Ilaló and North Chiche groundwaters in possible interaction with deep hot water reservoirs connected with several aquifer layers by faults.

Groundwater Age Estimates
Current mean tritium ( 3 H) value in Quito rainwater is 3 TU, which corresponds to the natural atmospheric level in the tropical zone (IAEA-GNIP, 2012). The absence of tritium in the water samples (except at Chirimoyas) suggests that the groundwater is older than the beginning of atmospheric nuclear bomb tests in 1950. Using a piston flow model would result with the infiltration having taken place anytime before 1950, whereas using a mixing model would produce a much older groundwater age (several hundreds of years).
TDIC δ 13 C values of both the Chiche and Ilaló aquifers are generally enriched compared to the biogenic reference established at Chirimoyas 1 spring. Considering the geological context, this enrichment in 13 C may be due to a deep volcanic CO 2 input. In the literature, δ 13 C values reported for deep volcanic CO 2 gas cover a large range, depending on its origin. Deep CO 2 can originate from (i) organic matter trapped in the sediments with a δ 13 C composition between −30 and −20‰ vs. VPDB (Evans et al., 1981), (ii) metamorphism of carbonate marine rocks between −3.7 and +2‰ vs. VPDB in δ 13 C (Barnes et al., 1984), and (iii) Earth mantle and crust with δ 13 C between −10.9 and −4.0‰ vs. VPDB, with a mean value around −6‰ (Kyser, 1986). The 13 C enrichment observed in Ilaló and Chiche groundwaters excludes the first and the third hypotheses because of the exceedingly negative values. Deep CO 2 could have its origin in the metamorphism of marine carbonates or recycled gas from magma produced in the subduction zone possibly connected to Quito faults system, which is a geological formation older than Ilaló volcano. As stated by Saffer and Tobin (2011), in subduction zones such as Ecuador, the persistent seismic movements and the high pressures, combined with the porosity and hydraulic characteristics of faults, can favor interactions between oceanic flow and other processes inside the continent. According to Alvarado et al. (2014), using geomorphological studies, GPS data and seismic methods, and Reyes et al. (2020), who worked in the same area using a magnetotelluric technique, this fault is a very old structure that is possibly connected to the oceanic terrains accreted to the continent. Schipper et al. (2017) measured the δ 13 C of CO 2 emanating from three volcanoes in the Andes of South Peru and Chile. They found δ 13 C values from − 3.7‰ to + 2.3‰ vs. VPDB. Their results suggest that deep CO 2 coming from these volcanoes have a connection with old oceanic structures. We used the range of values from Schipper et al. (2017) as our geogenic 13 C end member.
The initial interaction between water and CO 2 at high pressure and temperature in the volcanic system probably produces a strong acidification of water and the quantitative conversion of the dissolved gas to H 2 CO 3 , with negligible isotope fractionation between CO 2 and H 2 CO 3 . The two origins of TDIC in the system are therefore: a biogenic end member defined by the point with the highest 14 C activity (Chirimoyas 1 -C201, A 14 C = 103.3 pmc) where δ 13 C value is −10.62‰ vs. VPDB, and the geogenic pole with δ 13 C values between −3.7 and +2.3‰ vs. VPDB (Figure 11). Water ages were therefore calculated with the geogenic end-member spanning the −3.7‰ to +2.3‰ vs. VPDB range (Table 1), giving uncertainties of the order of ± 10% of the age computed.
Three wells intercepting the North Chiche aquifer line up along the south-to-north flow line showed in Figure 2A, namely FIGURE 9 | Korjinski diagram at 25 • C or alumino-silicate activity diagram (A) before equilibrium and, (B) after equilibrium when pyrophillite is not stable. Frontiers in Water | www.frontiersin.org FIGURE 11 | Measured 14 C activity (A 14 C) as a function of measured TDIC δ 13 C for North Chiche wells (blue points), Ilaló wells (red points), and mixed wells (green points). The vertical red line represents the biogenic TDIC δ 13 C signature. The purple arrow represents the acceptable range for deep geogenic TDIC δ 13 C values. The blue lines represent the range of mixing lines between the soil biogenic and geogenic end-members. A 0 values for all stations would stand on the biogenic/geogenic mixing line.
the stations C103 (El Carrizal), C101 (Las Acacias), and C104 (La Esperanza -MICEI). The station C105 (Patagua) stands approximately two km to the west from these stations. The water ages calculated for this group of stations show a clear northward increase, from a few hundred years for C103, to 2793 ± 51 years for C101 and 1728 ± 932 years for C105, and up to 3611 ± 1,036 years for the north most point C104. These estimates lead to a horizontal flow velocity of the order of 2 m.a −1 , one order of magnitude less than the 15 to 1,500 m.a −1 groundwater flow velocity estimated from hydraulic properties.
In the two other wells intercepting North Chiche (Diego Andrade -C109 and EEQ -C204) groundwater was dated in the same range, between 3,666 ± 388 and 1,696 ± 70 years. Ilaló groundwater was dated between 10,700 ± 1,000 and 45,700 ± 800 years, in spring I203 (Cununyacu) and the well I102 (El Tingo), respectively, and the mixed wells had intermediate apparent ages (7,300 ± 300 and 12,800 ± 1,400 years). The modern water from C201 (Chirimoyas 1) spring may be linked with a fractured zone in the northern part of the study area allowing direct recharge from the surface.
Groundwater hydrodynamic calculations suggest a mean age of a few decades for the North Chiche aquifer and the total transit time between recharge and the furthest discharge area was estimated at a maximum of 430 years. On the other hand, radiocarbon dating indicates ages up to ∼3,600 years. According to the 14 C method, the net infiltration has to be much lower than those calculated from hydrodynamics (7 mm.a −1 , minimum infiltration = 1 mm.a −1 , maximum infiltration = 20 mm.a −1 ), representing at most a few percent of the annual precipitation. On the other hand, if the hydrodynamic calculations are relevant, it implies that the TDIC contamination from CO 2 coming from the Ilaló volcano is much higher than expected, contaminating the Chiche groundwater, and that it bears a δ 13 C isotopic signature in a lower range (−4 to −8‰ PDB) than presumed. Other independent tracers would be necessary to date the groundwater. This significant difference in groundwater age may be attributed to the scarcity of data used to construct the geochemical and hydrodynamic hypotheses, the lack of knowledge on the extent of heterogeneity of the aquifer, and chiefly spatial variations of permeability, porosity and hydraulic conductivity coefficient distributions as well as few radiocarbon dating data.
To reconcile these large discrepancies, it can be proposed that less permeable layers may exist inside the Chiche formation. In the northern part of the aquifer, springs outcrop at different elevations, which suggests important vertical variations of permeability in the Chiche formation and the existence of distinct aquiferous layers supporting the hypothesis of a partly disconnected system. There is also few pumping test data to accurately estimate the distribution of permeability values. However, the pumping test of two neighboring wells located in North Chiche, along a flow line, gave consistent permeability values (i.e. K El Carrizal = 1.2 x 10 −5 m/s and K Las Acacias = 1.9 x 10 −5 m/s). In this particular case, a decrease of permeability toward the north is not well-corroborated.
If we consider a multilayer aquifer, we can calculate an average permeability taking into account the difference showed by the hydrodynamic and the isotopic tracers approaches. The estimated horizontal flow velocity is 2 m.a −1 using the apparent ages results, leading to an average permeability of 1.9 x 10 −6 m.s −1 , suggesting that the heterogeneity in the Chiche formation could be greater than expected. Using the same information, we can also calculate a recharge rate using Darcy law. In this context, Q = 0.03 m 3 .s −1 , 25% lower than the minimum rate we have obtained using the available hydrodynamic data.
Using the scarce information existing in this area and the simple approaches we used allowed to determine a range of hydrodynamic values characterizing the system. The Chiche formation appears to constitute a heterogeneous, partly disconnected aquifer, providing locally a good water productivity, albeit of questionable quality, possibly linked to deep circulations rather than horizontal flow. Future studies may be conducted focusing on the spatial characteristics of hydraulic discontinuities and the structural properties of the aquifer.

CONCLUSIONS
Our triple approach, using hydrodynamics, chemistry and isotopic tracers, allowed a better characterization of the Tumbaco -Cumbayá -Los Chillos aquifer system. Hydrodynamic results showed limited groundwater-level variations in both South and North Chiche, linked mainly to the confined state of the aquifers which prevents direct recharge, and general flow directions from south to north. The South Chiche aquifer has Mg 2+ -HCO − 3 water type whereas the North Chiche aquifer has Na + -Mg 2+ -HCO − 3 water type, similarly to the Ilaló aquifer. Stable water isotopes identified different elevations for the recharge area of the Chiche and Ilaló aquifers. Chemical data revealed different mineralization processes for the two aquifers with a geothermal effect for Ilaló and North Chiche aquifers. Chemical characteristics of Ilaló groundwater, in particular high As contents, were also met in North Chiche, pointing to a highly probable connection with the Ilaló aquifer. In this study arsenic constitutes a marker of volcanic imprint.
Radioactive isotopes showed low renewal rates with older ages for Ilaló aquifer compared to Chiche. The northernmost area (Chirimoyas spring -C201) showed a modern age suggesting fractured zones feeding this part of North Chiche aquifer by direct rainfall infiltration. This area may not be significantly connected with the main North Chiche aquifer, which could behave like a multilayer system, although no factual data confirm this characterization at present. The radiocarbon ages in North Chiche (ages between modern and ∼3,600 BP) are markedly different from groundwater transit time estimated using a simple hydrodynamic approach (between 5 and 430 years). This discrepancy in circulation pattern of groundwater may reveal the heterogeneity of the aquifer and a reduced global porosity, possibly due to ash layers present in the formation, that could form lenses of compacted material. The presence of deep faults is assumed and vertical and horizontal connections are unknown. More research is needed to clarify the overall structure of the aquifer.
Our results contributed to the better understanding of a complex continental volcanic aquifer, but further research on the Chiche and Ilaló aquifers is necessary, especially to better understand the geological heterogeneity, the hydrodynamic parameters variation in the area, and using other chemical tracers and gases analysis, to understand the input from the volcanic system and the arsenic dynamics. These new results could be extrapolated to other similar zones in the Ecuador or Andes cordillera. These findings highlight that knowledge progresses through the confrontation and the contradictions of all the available information. In these heterogeneous systems with high spatial-temporal variability, the field measurements are fundamental. It is worth noting that short data series can lead to erroneous conclusions. The large differences observed between the hydrodynamic and the geochemical approach in the study area could be an excellent stimulation to better investigate an aquifer system that has probably revealed just a little part of its complexity and heterogeneity.

DATA AVAILABILITY STATEMENT
The original contributions generated for the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
CM, JT, and CL conceived and designed the study. CM, JT, NP, and CL did the field work. CM, NP, and CC performed the laboratory work. CM acquired the samples, developed the calculations, and analyzed the results. JT, CL, and NP provided input on method and statistical analysis. JT and CC provided input on the methods to take samples for chemical analyses. All authors provided critical feedback, helped shape the research, analysis, contributed to the article, and approved the submitted version.

FUNDING
IRD supported this research with a financial aid for technical infrastructure and chemical and isotopic analyses. SENESCYT (Francia 2008-1 program)/IRD (files 702512B and 744360J) helped for student scholarship personal support. This research was supported by the biodiversity and natural resources network VLIR NETWORK Ecuador which also funded the publication fees of this article.