Discriminant model and hydrogeochemical processes for characterizing preferential flow paths in four interconnected volcanic aquifers in Costa Rica

The volcanic Barva and Colima multiaquifer system is crucial for the drinking water supply of about 1.7 million people, about 30% of the population in Costa Rica. The complexity of this system and the simultaneous occurrence of anthropogenic and natural processes complicate the understanding of hydrochemical and hydrogeological dynamics, both essential for groundwater protection. This research aimed to develop a model for the discrimination of groundwaters according to the main withdrawn aquifer and flow path, assess interaquifer connections, and evaluate the main hydrochemical processes governing water quality. Samples (571) from 38 sampling sites, collected quarterly from 2016 to 2020, were analyzed for nitrate, major ions, and silica. Principal component analysis and discriminant analysis exhibited and validated sample grouping according to the primary aquifer system captured, i.e., Upper Barva, Lower Barva, Upper Colima, and Lower Colima, and the occurrence of two flow paths within the Lower Barva and Upper Colima aquifers. Hydrochemical and statistical analyses showed resilience to seasonal chemical variation in deeper aquifers and also three groundwater mixing processes. Lower Barva groundwater enriches in bicarbonate, magnesium, calcium, sodium, and silica, mainly due to weathering of mafic andesitic-basaltic lavas. The Upper and Lower Colima showed higher silica and major ion content, except for calcium, indicating longer residence times, dissolution from felsic andesitic minerals, and calcium adsorption by normal ion exchange. Weathering of aluminosilicates is the primary process governing groundwater quality in the four studied aquifers, whereas cation exchange, interaquifer leakage, and anthropogenic processes might modify the groundwater chemistry.


Introduction
In many regions of the world, volcanic aquifers have become valuable resources in providing water for human consumption and economic activities (Bertrand et al. 2010;El Maghraby 2015;Hosono et al. 2020;Koh et al. 2018;Morán-Ramírez et al. 2016). Approximately 65% of the population in the Central Valley (Valle Central), Costa Rica, equivalent to 1.7 million people, depends on groundwater resources from the Barva-Colima volcanic aquifer system (Reynolds- Vargas et al. 2006). This multilayered system includes the Upper Barva (UB), Lower Barva (LB), Upper Colima (UC), and Lower Colima (LC). Despite the importance of groundwater for public supply in Costa Rica (Madrigal-Solís et al. 2020a), there is scarce information on the impact of the hydrogeochemical processes on water quality. Preliminary studies on groundwater hydrogeochemical characterization have been developed in the northwestern region of Valle Central (Central Valley;Madrigal-Solís et al. 2017), the southern part of Valle Central (Vargas and Mora 1999), the eastern-central region of the country (Madrigal-Solís et al. 2020b), and the Central Pacific coast (Sánchez-Gutiérrez et al. 2020).
According to the predominant minerals, water-rock interactions directly influence the geochemical evolution of groundwater, while other natural processes might modify its chemistry. The saturation of dissolved constituents and their further precipitation into secondary clay minerals, calcite or dolomite, may lower the concentration of calcium, magnesium, sodium, and other components in the liquid phase (Kim et al. 2020). The adsorption of calcium and magnesium on the clay surfaces releases additional sodium and potassium through the normal ion exchange process, while the reverse ion exchange may cause the opposite (Custodio and Llamas 1983). In the presence of additional CO 2 inputs to the aquifer, e.g., from degassing due to volcanic activity, high partial pressures of CO 2 may accelerate the dissolution of calcium bicarbonate, preventing calcite precipitation and leading to higher bicarbonate and calcium contents (Locsey et al. 2012). The leaching of chemical substances from anthropogenic sources also modifies groundwater chemistry. In the Barva and Colima aquifers, Madrigal-Solís et al.
(2017) described a general increase in chemical constituents as groundwater transits through volcanic formations due to the natural dissolution of minerals and anthropogenic activities, while Reynolds- Vargas et al. (2006) found a temporal increase of nitrate in several springs and wells extracting water from the Barva aquifer due to the extensive use of septic tanks and nitrogen fertilizers in coffee plantations.
Moreover, in the Barva-Colima system, the identification and characterization of each water level are challenging, as the four aquifers present variable depths and thicknesses, groundwater mixing processes (BGS/SENARA 1988;Madrigal-Solís et al. 2017), and similar rock chemical compositions (Losilla et al. 2001). Furthermore, the simultaneous occurrence of anthropogenic and natural geochemical processes in Barva and Colima aquifers hinders the relative importance of each process and complicates the identification of the captured aquifers in selected sampling sites. In this situation, determining the hydrogeological unit from which the water samples are withdrawn represents a significant problem when investigating the hydrochemical processes affecting the water quality of each aquifer. Hydrochemical and statistical analyses enable, as a first step, the discrimination of the hydrogeological unit, the identification of distinctive flow paths within each unit, and subsequent examination of the processes affecting the water quality in each aquifer. In this regard, linear and quadratic discriminant analysis (DA) is a multivariate statistical practical procedure broadly used to verify a preconceived classification of groundwater samples according to the corresponding hydrogeological unit. The DA generates a model by defining the most critical variables that differentiate these units, providing the linear or quadratic discriminating functions that distinguish each aquifer, and classifying each sample according to these functions. DA has been successfully applied in volcanic systems, combined with other multivariate methods (Barzegar et al. 2019;Panagopoulos et al. 2016) and machine learning procedures (Ha et al. 2021;Sajedi-Hosseini et al. 2018). Moreover, principal component analyses (PCA; Gan et al. 2018;Panagopoulos et al. 2016;Pazand and Javanshir 2016), correlation analyses (Morán- Ramírez et al. 2016;Kim et al. 2020), factor analysis (Gan et al. 2018); hierarchical cluster analysis (Agbasi and Egbueri 2022; Egbueri 2020; Gan et al. 2018); hydrochemical diagrams (Ebrahimi et al. 2022;Fenta et al. 2020;Hosono et al. 2020;Shishaye et al. 2020;Wisitthammasri et al. 2020), base exchange indices (Haji et al. 2021;Mechal et al. 2017), multiple linear regression (Egbueri and Agbasi 2022), and the assessment of saturation indexes of aqueous mineral phases (Esteller et al. 2017;Hosono et al. 2020;Morán-Ramírez et al. 2016) are widely employed to facilitate the interpretation of geochemical processes with large hydrochemical data series in volcanic aquifers.
Although most research on groundwater in Costa Rica focuses on the Central Valley, previous studies provide a general characterization of the hydrogeology and water quality in the UB, LB, UC, and LC aquifers (e.g., Arias-Salguero et al. 2006;BGS/SENARA 1988;Foster 1993;Gómez-Cruz 1987;Losilla et al. 2001;Protti 1986;Ramírez and Alfaro 2002;TAHAL 1990). Likewise, Reynolds-Vargas and Fraile-Merino (2009) and Sánchez-Murillo et al. (2022) estimated the approximate ranges of recharge elevations and apparent groundwater age for several sites of the Barva and Colima aquifers, while Madrigal-Solís et al. (2017) provided a preliminary hydrogeochemical characterization. Nevertheless, these investigations do not supply in-depth information on the hydrogeochemical processes or interaquifer connections governing groundwater quality, e.g., hydrolysis of aluminosilicates, ion exchange, mineral precipitation, and mixing of groundwaters. Thus, the present research proposes a comprehensive understanding of the controls on water quality along the flow paths in Barva and Colima aquifers, which is necessary to contribute to the sustainable use of groundwater resources.
Given that groundwater management requires sufficient knowledge of the hydrogeochemical dynamics of the aquifers captured for drinking water supply, this research aimed to (1) develop a descriptive hydrogeological model for the classification of groundwaters in terms of the corresponding aquifer and preferential flow path and (2) assess the main hydrogeochemical processes, e.g., water-rock interactions, groundwater mixing processes, ion exchange, and equilibrium processes governing the quality of water in the Barva and Colima multilayer aquifer system. The results will provide a better understanding of the fundamental hydrogeochemical behavior and the controls on water quality along the flow paths of the most significant aquifers in Costa Rica for implementing improved management policies. Finally, recognizing representative sampling sites of the typical conditions of each Barva and Colima aquifer along the flow paths is essential for interpreting chemical and isotopic data in the present and future hydrogeological studies, e.g., contamination processes, identification of zones of recharge, and transit times, especially in the presence of complex multiaquifer systems with similar geological characteristics.

Study site
The study area is located in the central and northern zones of the Valle Central of Costa Rica (9°57′31″-10°8′1″ N latitude and 84°3′44″-84°15′17″ W longitude). It includes the Barva and Colima multiaquifer system, i.e., Upper Barva (Los Angeles and Bambinos subaquifers), Lower Barva, Upper Colima, and Lower Colima aquifers ( Figs. 1 and 2). The area exhibits an average slope of 14%. Forest and pastures dominate in the uplands, accounting for 24 and 15%, respectively, of the total area overlying the shallower Barva aquifer. In the middle zones of the subbasins, from 1,000 to 1,400 m asl, the main land uses include coffee and other agricultural activities, accounting for 26 and 12% of the total area, respectively. About 23% of the area was covered by urban and industrial activities, mostly from 800 to 1,000 m asl (Madrigal-Solís et al. 2019).
The rainy season extends from May to November-mid-December and the dry season from mid-December to April (Maldonado et al. 2013;Taylor and Alfaro 2005). Scattered rainfall arrives from the Caribbean slope from November to April, mainly in the upper zone. Precipitation is highly altitude-dependent. From a data set comprising 54 years of observations, mean values in the upper area range between 90 mm (SD = 22) and 370 mm (SD = 81) during the dry and wet seasons, respectively, while mean values in the lower zones vary from 18 mm (SD = 20) to 257 mm (SD = 88; IMN 2014). Temperature is relatively stable throughout the year and highly dependent on altitude (Taylor and Alfaro 2005); in the upper zone, close to 10 °C, while in the lowlands, between 20 and 24 °C (IMN 2014).

Geology
The geological and hydrogeological settings in the study area are described in detail by Williams (1952), Echandi (1981), BGS/SENARA (1988), and Denyer and Arias (1991). Based on the stratigraphic terminology by Williams (1952), the area comprises, from bottom to top (Denyer and Arias 1991;Protti 1986;Fig. 2): (1) Colima Formation, andesitic lava flows covering most of the Valle Central, including La Libertad Member (two lava flows, covering a minor sector at the southeastern part of the study area and, thus, not included in this investigation or Fig. 2), Lower Colima Member (several lava flows), Puente Mulas Member (layer of tuffs and ignimbrites), and Upper Colima Member (sequence of lavas and pyroclastites); (2) Tiribí Formation, pyroclastic deposits of tuffs and ash-flows covering most of the Valle Central; and (3) Barva Formation, subdivided into six lithological units: Bermúdez Member (several fractured lava flows), Porrosatí and Carbonal members (layers of pyroclasts, weathered lithic and ash tuffs, and lahars), Los Bambinos Member (lava flows), Los Ángeles Member (lava flows), and Crater Member (thick ash and pyroclasts). In addition, at the southwestern margin of the study area, Debris avalanche El Coyol, a superficial deposit of polymictic volcaniclastic material associated with a debris flow event, is overlying Tiribí Formation (Méndez and Hidalgo 2004).
The origin of the formations is related to asynchronous eruptive events, i.e., of different ages, explaining the modest geochemical differences in their chemical composition. In the Colima Formation, the lavas of the Lower Colima Member are trachyandesites and porphyritic andesites, containing phenocrystals of plagioclase, e.g., anorthite, Puente Mulas Member is composed of tuffs and ignimbrites, whereas the lavas from the Upper Colima Member are aphyric andesites with only 4% phenocrysts (Kussmaul 1988)- Table S1 of the electronic supplementary material (ESM). The Tiribí Formation corresponds to a pack of chemically variable tuffs and ignimbrites, ranging from andesites transitional to shoshonitic dacites (Kussmaul 1988;Kussmaul and Sprechmann 1982), basaltic-andesites, trachytes, and trachydacites (Hannah et al. 2002), with trachyandesites being the dominant composition (Pérez et al. 2006), bearing 1-10% of plagioclases (Table S1 of the ESM). In the Barva Formation, the lavas of the Bermudez Member exhibit highly fractured porphyritic basaltic andesites (Arredondo Li and Soto 2007;Echandi 1981) with plagioclase and anorthite phenocrysts (Table S1 of the ESM). The lavas of the Los Bambinos Member are porphyritic andesites, showing a higher proportion of plagioclase phenocrysts than Los Angeles (Arredondo Li and Soto 2007; Table S1 of the ESM), while Los Angeles Member lavas are basaltic andesitic (Rojas et al. 2017), of porphyritic aphanitic texture (Arredondo Li and Soto 2007).
In terms of major and trace elements, Tiribí tuff rocks are similar to Colima Formation, indicating a common source, while the younger lavas from Barva Formation are more mafic (Pérez et al. 2006). The lavas from Colima Formation and Tiribí differ from the Barva Formation since (1) in Colima and Tiribí, the K 2 O/Na 2 O ratio is approximately one and higher than in Barva Formation, (2) Colima and Tiribí contain more K 2 O and Na 2 O than Barva Formation at the same silica concentration, although all three formations are considered K-rich and, (3) Barva Formation contains more CaO and MgO than Colima and Tiribí, being Barva Formation, a typically calc-alkaline association (Kussmaul 1988).

Hydrogeology
The fractured lavas and high permeability breccias in Valle Central form the Barva and Colima multiaquifer system (Fig. 2;BGS/SENARA 1988). The Barva multiaquifer system (Figs. 1 and 2), hereafter referred to as Barva aquifer, comprises an area of approximately 190 km 2 . It comprises: (1) Upper Barva (UB) aquifer, which includes the Los Bambinos and Los Angeles sub-aquifers embedded in the Los Bambinos and Los Angeles Members, and (2) Lower Barva (LB) aquifer within the Bermúdez Member (Ramírez and Alfaro 2002). The transmissivity of the aquifer ranges between 100 and 500 m 2 /day, with a storage coefficient of 0.1 (Gómez-Cruz 1987;Ramírez and Alfaro 2002); hydraulic conductivity varies between 1 and 10 m/day, the saturated thickness between 50 and 100 m (Ramírez and Alfaro 2002), and the depth to water level between 25 and 60 m (Arias-Salguero et al. 2006). The recharge of UB occurs by direct infiltration and influent rivers (Ramírez and Alfaro 2002). In LB, the recharge occurs by infiltration from influent rivers, percolation from UB through the Carbonal and Porrosatí aquitard, and direct infiltration in the remaining areas (BGS/ SENARA 1988;Gómez-Cruz 1987). Preferential recharge occurs in the upper parts of the aquifer (Reynolds- Vargas and Fraile-Merino 2009;Sánchez-Murillo et al. 2017;Sánchez-Murillo et al. 2022). Discharge occurs mainly through springs, as the base flow of rivers, and through the fractured tuffs and ignimbrites of the Tiribí aquitard, which are sufficiently permeable to allow infiltration from the Barva aquifer into the Colima underlying aquifers (BGS/SENARA 1988;Foster et al. 1985;Reynolds-Vargas and Fraile-Merino 2009). The Colima multiaquifer system is situated within the Colima Formation and separated from the Barva aquifer by the Tiribí aquitard (Fig. 2). It comprises the aquifers (1) Upper Colima (UC), of approximately 170 km 2 , and (2) Lower Colima (LC), of 230 km 2 (Arias-Salguero et al. 2006). Tuff and ignimbrites of Puente Mulas aquitard separate both aquifers (Fig. 2). Their boundaries have not been established due to similarities in rock compositions with Barva aquifer and their varying thickness and depths (BGS/SENARA 1988;Losilla et al. 2001). The UC exhibits transmissivities between 500 and 9,500 m 2 /day (Losilla et al. 2001). It is confined when overlain by the Tiribí Formation, while functioning like an unconfined aquifer in the remaining extension. In the western portion of the Valle Central, it recharges from the Barva aquifer, through the fractured tuffs and ignimbrites of the Tiribí Formation (Arias-Salguero et al. 2006;BGS/SENARA 1988) and by direct infiltration in areas where it is not confined (TAHAL 1990). The LC aquifer is confined by the Puente Mulas aquitard, whereby part of its recharge occurs by percolation through this unit and part by direct infiltration (TAHAL 1990). The LC exhibits a storage coefficient of 1.4 × 10 -4 and a transmissivity ranging from 500 to 7,500 m 2 /day (Losilla et al. 2001). The discharge from the Colima aquifers occurs through springs and along the Virilla River.

Water sampling and analyses
A total of 38 sites, including 10 springs and 28 wells with depths ranging from 30 to 350 m, were sampled (Table 1). From 2016 to 2020, 571 samples were collected approximately every 3 months during March-April (representative of the dry season), June (early wet season), September (midwet season), and November-December (late wet season). The physical parameters were determined in situ with a multielectrode HI98311 model (Hanna Instruments, USA) for electrical conductivity (EC in μS/cm), total dissolved solids (TDS), and temperature (T). A multielectrode HI98121 was used for pH and redox potential (ORP) measurements, and a HI 9147 sensor for dissolved oxygen (DO). Samples were collected in duplicate in 500-ml Nalgene bottles and kept refrigerated at 4 °C until analysis the following day. Bicarbonate (HCO 3 − ) was estimated by volumetric analysis according to method 2310B (American Public Health Association 2012), whereas calcium (Ca 2+ ), magnesium (Mg 2+ ), potassium (K + ), sodium (Na + ), chloride (Cl -), nitrate (as NO 3 -), sulfate (SO 4 2-), and fluoride (F -) were analyzed by ion chromatography (Dionex Thermo Scientific ICS 5000), according to the method 4110 (Rice 2012). In addition, silica (SiO 2 ) was analyzed in 53 samples collected during June 2016, January 2021, and April 2021 by spectrophotometry according to the 4500-SiO 2 C method (Rice 2012). Only samples with a charge balance error of 5% or less were used.

Statistical analyses
Hydrochemical and physical variables from groundwater samples were assessed through the statistical software SPSS (version 25) to conduct: (1) descriptive analyses, (2) normality and homogeneity of variances tests using the Shapiro-Wilk and Levene's tests, (3) comparisons between wet and dry seasons within each sampling site to assess whether seasonal or annual data should be used in subsequent statistical and hydrochemical analyses, (4) comparisons between wet and dry seasons within each aquifer and flow path and, (5) comparisons between the four studied aquifers and flow paths. A t-test was applied to compare two data sets when data complied with the assumptions of the normal distribution; otherwise, the Mann-Whitney test was used.
The lithological descriptions were obtained from the National Service for Groundwater, Irrigation and Drainage, SENARA (Costa Rica), and from the Empresa de Servicios Públicos de Heredia, ESPH S.A. (Costa Rica) databases. For the selection of the sampling sites, preliminary identification of the primary aquifer was performed based on the approximated depths of the boreholes and the available hydrogeological description. The principal component analysis (PCA) and discriminant analysis (DA) were applied to assess and validate this prior classification, differentiate the physical and chemical variables in each aquifer, and evaluate the presence of different flow path systems by using the statistical software XLSTAT (Addinsoft 2021). The physical and chemical mean values were log-transformed (Log 10 ) and standardized (z-scores, where z i = (x i -mean)/SD) to homogenize the effect of variables measured on different scales and obtain a log-normal distribution (Reimann et al. 2008).
The DA was employed to validate the prior classification of samples, according to the aquifer withdrawn, using the available hydrochemical data from all samples. The DA algorithm eliminated variables showing multicollinearity from the analysis. It tested the discriminant ability of the selected independent-predictor variables by examining the differences in means among groups using Wilk's Lambda test (Rao's approximation). A Box test (with a Chi-square asymptotic approximation) was employed to select the DA model. When equality among the covariance matrices is present, the linear model should be selected; otherwise, the quadratic model should be chosen. The DA estimated the group centroids; a sample close to a centroid was classified as belonging to that group. Discriminant functions classified the samples into each group. Furthermore, a validity test was performed, in which the discrimination functions were estimated for all except for 20% of the samples, i.e., validation samples. Then, the discriminant functions were used to predict the classification of the randomly selected validation samples. The validity of the model is verified through the percentage of succeeded classifications.

Hydrogeochemical diagrams
Hydrochemical data of all individual samples were used to produce Piper diagrams (Piper 1953), Pearson correlation coefficients, and bivariate diagrams to assess the influence of water-rock interactions, mixing, and ion-exchange processes on groundwater quality. The mean values were included in the diagrams to facilitate the graphical representation of each aquifer. In the analysis of ion-exchange processes, two base exchange indices were used: chloro-alkaline imbalance indices 1 and 2 (CAI-1 and CAI-2), according to Eqs. (1) and (2) (Schoeller 1967): These indices usually obtain a positive or negative value close to 0. The exchange of bases is defined by the existence of an increasing or decreasing trend, not by an absolute value (Custodio and Llamas 1983). Negative values of CAI-1 and CAI-2 indicate normal cation exchange, i.e., calcium and magnesium are adsorbed onto clay minerals while sodium and potassium are released into groundwater. Positive values of CAI-1 and CAI-2 suggest reverse ion exchange, i.e., sodium and potassium retention on clay minerals and release of calcium and magnesium into the groundwater (Haji et al. 2021). In addition, a higher absolute value of CAI-1 indicates a higher degree of cation exchange (Pan et al. 2019). Moreover, the relation of [Na + + K + -Cl -] vs. [Ca 2+ + Mg 2+ -SO 4 2--HCO 3 -] was employed to corroborate the relevance of ion exchange. A linear relation with a slope of -1.0 indicates ion exchange is an important control of water chemistry (Fisher and Mullican III 1997;Fenta et al. 2020;Hosono et al. 2020). Finally, the software package EQ3/6 for hydrogeochemical modeling (Wolery 1992) was used to estimate the saturation indexes and partial pressure of CO 2 (g) in 53 samples collected during June 2016, January 2021, and April 2021.

Results
The hydrochemical parameters of the 571 groundwater samples collected in 38 springs and wells in the study area are shown in Table S2 of the ESM, with the minimum, maximum, mean, median, standard deviation (SD), and coefficient of variation (CV) estimations. The Shapiro-Wilk test showed evidence of normal distribution in almost all constituents within each sampling site (Table S2 of the ESM).

Principal component analysis
The PCA was applied to evaluate relationships between variables in each aquifer and the potential occurrence of different flow path systems. Mean values were employed since, for most of the chemical parameters: (1) the Shapiro-Wilk test did not show evidence of nonnormality (Table S2 of the ESM), and (2) the Mann-Whitney U test found no significant differences between seasons per sampling site (Table S3 of the ESM). Variables with factor loadings explaining most of the cumulative variance in the first two principal components (PC) were selected, including major ions, nitrate, and temperature, all of which showed factor loadings above 0.6. The Bartlett Chi-square test confirmed that variables are correlated but not orthogonally (df = 36, p < 0.001), whereas the Kaiser-Meyer-Olkin measure of adequacy was 0.755. The first two principal components described 85.89% of the cumulative variance of the data set (Table 2), with eigenvalues >1.
The first function F1 showed strong factor loadings in TDS, calcium, magnesium, sodium, potassium, bicarbonate, chloride, sulfate, and temperature (from 0.715 to 0.910; Table 2; Fig. 3a), representing mineralization due to the weathering of volcanic lithology. The second function F2 exhibited a strong positive loading only in nitrate (0.849) and, to a lesser extent, in calcium (0.463), representing anthropogenic activities such as inputs from nitrogen fertilizers and sewage leaching from septic tanks/ drainage systems. Weaker negative loadings were found in sulfate, sodium and potassium (from -0.451 to -0.611), showing that the relative importance of mineralization due water-rock interactions is negatively correlated to F2, i.e., to anthropogenic activities.
The F1-F2 plot exhibited three subgroups of samples within the Barva aquifer (Figs. 3a and 4a): 1. The UB samples and LB samples in the highlands, located in the negative portion of F1 and F2 and opposed to variables related to water-rock interactions. These samples represent young groundwaters in the preferential recharge zone, and thus low mineralization and temperatures, small transit times and short flow paths. 2. Samples along a northern flow path, hereafter referred to as Lower Barva flow path 1 (LB1), located mainly in the positive part of F2 and close to nitrate. 3. Samples along a southern flow path, referred to as Lower Barva flow path 2 (LB2), located on the positive side of F1, close to nitrate, but also to calcium, chloride, magnesium, and bicarbonate variables.
The F1 and F2 plots also showed two subgroups among the Colima samples (Figs. 3a and 4b,c):  1. Samples along a north-south flow path, referred to as Upper Colima flow path 1 (UC1), were strongly related to potassium and sodium, increasing their contents as groundwater flows from RN07 and RN35, at a higher elevation, to discharge site RN43. 2. Samples located along a northeast-southwest flow path, referred to as Upper Colima flow path 2 (UC2), in the positive section of F1 and F2. These samples are related to nitrate, calcium, chloride, magnesium, and bicarbonate, and lay between LB2 and UC1. 3. Samples from LC are highly related to potassium and sodium and are closely associated with the UC1 samples. As opposed to the other flow paths, mineralization decreased as groundwater traveled from sampling site RN09, at a higher elevation, to discharge site RN15, which is also evident by the decrease in the mean TDS values, from 213 in RN09, to 123 ppm in RN15 (Table S2 of the ESM).

Discriminant analysis
Discriminant analysis (DA) was conducted to validate the prior classification of samples into aquifers and flow paths, using the same set of variables used in the PCA as independent variables. The box test indicated that within-class covariance matrices are different, and the quadratic model should be used (df = 140, p < 0.0001). The variables selected proved significant differences among group means (Wilk's Lamba test, df = 35, p < 0.0001) and, thus, their discriminant ability. The eigenvalues were estimated for the principal factors and the corresponding discrimination percentages, where the F1 described 50.1% of the discrimination power and the F2, 39.1%, leading to a cumulative discrimination power of 89.2% (Fig. 3b).
According to the correlation of the predictor variables and the factor loadings, i.e., estimated by the canonical functions, sodium, potassium, bicarbonate, magnesium, and chloride were the most significant variables with discriminant power in the F1 and nitrate and calcium in F2. The canonical correlations showed values of 0.869 and 0.835 for F1 and F2, respectively, showing similar discriminant power (Table 3). Once the canonical discriminant functions were estimated, the centroids evaluated these functions by showing a separation between the mean points of each grouping class (Fig. 3b).
The results showed similarities with those obtained from the PCA (Fig 3a): 1. In the discriminant factor space (Fig. 3b), the separation between centroids and their confidence circles distinguished the proposed grouping of samples according to the primary aquifer and flow path, e.g., LB1 from LB2 and UC1 from UC2. 2. The centroid of UB, in the positive part of F1 and F2, was significantly separated from the other centroids and opposed to the chemical variables, showing lower mineralization and shorter flow paths. 3. The centroids of LB1 and LB2, in the positive part of F1, displayed a stronger influence from the nitrate variable, while calcium showed a stronger influence on LB2. 4. The centroid of UC2, in the negative part of F1 and F2, is equidistant to LB2 and UC1, showing a plausible mixing pattern with LB2 and, thus, a more significant influence of calcium and nitrate compared to UC1.
The UC and LC centroids, in the negative part of F1, exhibited stronger influence from magnesium, chloride, sodium, and potassium, indicating longer flow paths, residence times and mineralization. According to the confusion matrix, the percentage of correctly classified samples was 92.8% from a total of the 456 training samples used to estimate the discriminant functions of the DA model (Table 4). The test of validity on the 114 remaining samples, i.e., 20% of all samples, showed that 90.4% of these samples were successfully classified (Table 4).

Hydrochemical characterization
The sampling sites and the proposed flow paths are shown in Fig. 4. The mean values were employed in the diagrams to facilitate the graphical representation of each aquifer. In general, the ion abundance of major cations and anions, in meq/L, in each aquifer and flow path was the following: According to TDS contents, ion abundance and the typical sequences for recently recharged waters, i.e., Ca 2+ > Mg 2+ > Na + > K + and, HCO 3 -> SO 4 2-> Cl - (Chebotarev 1955;Custodio and Llamas 1983), only the samples from UB were classified as slightly mineralized young groundwater. In LB and UC2, calcium was the dominant cation, while magnesium dominated in UC1 and LC. Regarding anions, bicarbonate was dominant in all flow paths. Nitrate was the second anion dominating LB and UC2, whereas UC1 and LC showed a dominance of sulfate and chloride after bicarbonate (Table 5), departing from the typical sequence for young, moderately mineralized groundwaters. The Piper diagram in Fig. S1 of the ESM shows individual groundwater samples and their median concentrations per sampling site. Groundwater from UB was classified into two facies, from mixed Ca 2+ -Mg 2+ -SO 4 2--HCO 3 to Ca 2+ -Mg 2+ -HCO 3 -. Most of the samples from the other aquifers were classified

Variation on hydrochemical parameters
Regarding physical parameters and major anions, important variations were found within and between aquifers (Table 5; Fig. 5; Tables S2 and S4 of the ESM). The mean values of ORP and DO were relatively high (217 mV and 6.27 mg/L, respectively), suggesting general oxidic conditions in all aquifers. The EC ranged between 33 to 728 μS/cm, while the pH ranged between 5.2 to 7.9, with a mean value of 6.58. In general, calcium content varied from 4.5 to 37.0 mg/L; sodium, from 2.8 to 32.4 mg/L; magnesium, from 1.1 to 18.4 mg/L, and potassium, from 1.1 to 12.7 mg/L, with CVs varying from 14% for calcium in LC to 49% for magnesium in UB (Table 5). The major anions also exhibited a wide range of values. Bicarbonate varied from 19 to 173 mg/L; nitrate, from 0.3 to 55.6 mg/L; chloride, from 0.9 to 31.6 and, for sulfate, from 0.7 to 38.4 mg/L, with CVs ranging from 14% for bicarbonate in LC to 83% for nitrate in UB (Table 5).

Comparison between aquifers and flow paths
When comparing aquifers and flow paths, the Mann-Whitney U test showed significant differences in almost all physical and chemical parameters (Table S4 of the ESM), including differences between LB1 and LB2 and between the UC1 and UC2, except for potassium (median (M) = 4.5 mg/L). Unexpectedly, more similarities were observed between LB2 and UC2 regarding nitrate (~20 mg/L), calcium (~21 mg/L). and sulfate (~6 mg/L; p < 0.001). The UB showed, in general, the lowest contents of TDS, EC, and major ions (p < 0.001), whereas LC exhibited the highest median values, except for nitrate and calcium (p < 0.001). Median chloride concentrations in LC doubled the contents found in UC and were three and nine times higher than in LB and UB, respectively (Table 5).

Seasonal variation
Seasonal differences were not significant when comparing dry versus wet season datasets within aquifers ( Fig. 5

Spatial evolution in aquifers and flow paths
Hydrogeochemical spatial variation is observed in Stiff diagrams in Fig. 6a, b, silica and magnesium concentration maps in Fig. 6c, d and biplots 7a-d. In Barva aquifers, low   major ion and silica contents were observed in the highlands and in UB and increasing gradually towards the lowlands (Fig. 6a, c). Bicarbonate, calcium, and magnesium exhibited the greatest increases (Fig. 7a). Groundwater in UC1 and UC2 also showed increasing concentrations along the flow paths, except for calcium in UC1, thus the sites at lower elevation displayed higher magnesium relative to calcium contents (Figs. 6b, d, 7b, c). In UC2, calcium maintained its dominance as the main cation (Fig. 7c), as in the Barva aquifer. Conversely, major ions and silica in LC decreased as groundwater traveled westward to lower zones, from RN09 to RN15 (Fig. 6b, d), while showing lower calcium relative to magnesium content (Fig. 7d) and, thus, similarities with UC1. Distances along the flow paths were estimated based on the first sampling site of each flow path and the corresponding phreatic and piezometric contour line (Fig. 4).

Weathering of aluminosilicate rocks
Dissolved CO 2 increases the acidity of groundwater through the production of weak carbonic acid (H 2 CO 3 ), enhancing the bicarbonate contents, the alteration of the minerals in the rocks, and the dissolution of chemical species, which, in turn, releases more bicarbonate. The diagram of [calcium + magnesium + sodium + potassium] vs.
[bicarbonate] (Fig. 7e) shows the relative importance    of bicarbonate as a balancing anion (Haji et al. 2021;Mechal et al. 2017). Samples fell below the 1:1 equiline and progressively departed from the line as groundwater is more mineralized, implying that bicarbonate is the dominant ion to balance major cations only in the recharge zone. Moreover, the chloride vs. bicarbonate/ chloride biplot (King et al. 2014) displays a lower relative abundance of bicarbonate, while chloride increases in LB, UC1, and UC2, which can be used to describe larger transit times (Fig. 7f).  Samples that fell around the stoichiometric lines for minerals dissolution denote a common origin of the involved parameters through the weathering of the dominant geologic formations and reveal additional sources of ions due to geochemical processes, other than rock dissolution in the study area. Hydrolysis of silicates and aluminosilicates yields kaolinite clays and major ions. Figure 8a shows UB and LB samples mainly falling between the 1:2 line for the anorthite dissolution (Eq. 3) and the 1:3.4 line for the augite weathering (Eq. 4) in the calcium: bicarbonate bivariable plot, whereas an excess in bicarbonate is observed with respect to calcium-rich plagioclase dissolution on a molarto-molar relation according to the dissolution reactions: (3) CaAl 2 Si 2 O 8 (anorthite) + 2CO 2 + 3H 2 O → Al 2 Si 2 O 5 (OH) 4 (kaolinite) + Ca 2+ + 2HCO 3 − The UC1 and LC samples lie slightly above the 4:1 line for pyroxene (hypersthene) dissolution (Fig. 8a), showing a major dissolution from ferromagnesian minerals, based on Eq. (5): In Fig. 8b, samples from UB and LB lie below the 1:1 line of calcium:sodium for the dissolution of calcium-sodium plagioclases (i.e., from andesite to oligoclase), indicating an (4) CaMg 0.7 Al 0.6 Si 1.7 O 6 (augite) + 3.4CO 2 + 4H 2 O → 0.3Al 2 Si 2 O 5 (OH) 4 (kaolinite ) + Ca 2+ + 0.7Mg 2+ + 3.4HCO 3 − + 1.1H 4 SiO 4 (silicic acid)  excess of calcium probably from calcium-rich plagioclase dissolution (anorthite) or anthropogenic sources, while UC2 and LC samples fell above, probably due to a major contribution from sodium-rich plagioclases, as follows (Eq. 6): In Fig. 8c, samples from UB and LB fell around the 0.7:1 line for magnesium:calcium ratio for augite dissolution, whereas Colima samples lie significantly to the right from the augite and fall closer to the general dissolution line for pyroxenes, indicating significant contribution from this mineral according to Eq. (7): In the biplot of calcium vs. sulfate (Fig. 8d), some samples from deep wells and others close to the Barva volcano in UB, UC1, and LC lay to the right of the 1:1 line for anhydrite. In contrast, LB and UC2 samples lay far below, indicating the low presence of sulfate-rich minerals to calcium. Samples from LB1, LB2, and UC2 showed high correlations between calcium and sulfate (Fig. 8d). Figure 8b-d also suggests a deficit in calcium in UC1 and LC with respect to the mineral dissolution.
Hydrogeochemical relations of Ca 2+ /(Ca 2+ + SO 4 2-) and Mg 2+ /(Mg 2+ + Ca 2+ ) with high molar ratios (> 0.5) denote the importance of calcium-rich plagioclase and ferromagnesian minerals weathering, respectively (Locsey et al. 2012   only in 54 and 92% of samples in UC1 and LC, respectively, while mean values are also higher in these flow paths. Figure 8e shows samples from UB, LB, and UC2, at higher elevations, located in the negative part of the CAI-1 and CAI-2 diagram, and samples at lower elevations in the positive part of the plot. To verify the importance of ion exchange, the diagram [Na + + K + -Cl -] vs. [Ca 2+ + Mg 2+ -SO 4 2 -HCO 3 -] showed slopes from -1.15 to -1.29 in LB2, UC1, and LC, relatively close to -1 (Fig. 8f). In contrast, the slopes in UB, LB1, and UC2 clearly deviated from -1 (-0.58, -2.7, and 0.7, respectively; Fig. 8f). Additionally, low molar ratios (<0.5) of calcium/bicarbonate might indicate cation exchange (adsorption of calcium and magnesium and release of sodium) or bicarbonate enrichment due to silicate weathering, while higher molar ratios may indicate an excess of calcium and magnesium from other sources, e.g., due to reverse ion exchange (Drever 1997;Hounslow 1995). From 87 to 100% of the samples in UB, LB1, and LB2 exhibited high molar ratios (>0.5), with LB1 showing the higher mean value (0.87). A lesser percentage, approximately 67% of samples, showed a high molar ratio in UC1. Conversely, low molar ratios (<0.5) were observed in 83% of all samples in LC.

Mineral saturation indexes
The silica in the 53 samples employed to estimate the saturation indexes (SI) ranged from 15.5 to 107.2 mg/L (Table S5 of the ESM). Groundwater is enriched with silica (Table S5 of the ESM), particularly during the dry and early wet season, resulting in most samples being oversaturated in all aquifers (SI from -0.75. to 1.2) with respect to silicate minerals, including tridymite (Fig. 9a), quartz (Fig. 9b), chalcedony, cristobalite (slightly lower than for chalcedony), and coesite, and slightly undersaturated for SiO 2 (am) ( Table S5 of the ESM). The SI in the Barva and Colima groundwaters increased close to equilibrium with respect to calcite (Fig. 9c), dolomite, aragonite, and magnesite minerals (Fig. S4a-d of the ESM) towards the lower basin, particularly in LC, while one sample in RN09 showed saturation with respect to dolomite ( Fig. S4a and Table S5 of the ESM). Partial pressures for dissolved CO 2 exhibited higher values towards the lower basin (Fig. 9d), while negatively correlated to pH (Fig. S4e of the ESM) and positively correlated to calcium (Fig. 9e). Negligible correlation coefficients between fluoride and calcium in all aquifers and flow paths (from 0.01 and 0.46) were found.

Seasonal variation
Statistical analyses of samples collected every 3 months for 5 years suggested that seasonality may not be a major control of groundwater chemistry during the study period since no significant differences were found between dry and wet seasons per site or per aquifer. However, higher site-specific CV in the shallower UB than sites from other aquifers suggest a more important seasonal variation in the uplands (Fig. S3 of the ESM). The lesser hydrochemical variability within sampling sites from LB and Colima aquifers indicates the occurrence of important resilience to temporal variations in recharge, typical of larger flow systems (Clark 2015). These conditions are especially evident in LC aquifer, being the only one showing no significant seasonal differences in any parameter.
Conversely, higher physical-chemical and hydrochemical variability in UB groundwaters reflects lower inertia in the recharge zones, where flows may be shallower and shorter and, therefore, more dependent on seasonality (Clark 2015), despite the similarities between dry and wet seasons within each sampling site. The hydrochemical temporal variation in UB agrees with the seasonal fluctuations in the water table, larger than 5 m, described by Foster et al. (1985) and Foster (1993) for Barva aquifer. However, Salas-Navarro et al. (2019) described a "damping effect" based on the minor temporal variations in daily and weekly isotopic composition in two springs above 2,200 m asl in UB, compared to more considerable variations in rainfall composition. EC remained close to the mean in Flores Spring, showing no relation to the seasonal pattern in precipitation or volumetric discharge, while the opposite was observed in Sacramento Spring at higher altitude. In deeper Colima aquifers, Foster et al. (1985) and Foster (1993) found slight or no seasonal variations in groundwater levels due to the buffering effect caused by the high storage capacity of the interbedded tuffs, which contributes to steady recharge rates. This coincides with the observations of Vega-Arce (2019), who conducted a study using groundwater level data from 2013 to 2018 in 12 wells in the study area and found no seasonality in all the seven deep wells that captured a mixture of UC and LC or exclusively LC.

Discriminant model
The discrimination between the two subgroups of samples within the LB, i.e., LB1 and LB2 flow paths, and the UC, i.e., UC1 and UC2 flow paths observed in the PCA (Fig. 3a), validated the DA model due to the separation of the group centroids and their confidence circles on the factor axes (Fig. 3b). The smaller distance between the LB2 and UC2 centroids in the discriminant factor space, even smaller than the distance between LB1 and LB2, evidenced high similarities in their chemical characteristics, which are more influenced by nitrate and calcium variables. The correct classification of most training and validity samples denotes the predictive capacity of the DA model to sort preclassified groundwater samples such as in other studies carried out in volcanic aquifers (Barzegar et al. 2019;Panagopoulos et al. 2016). Thus, the DA provided a validated model to classify old and new samples from the Barva and Colima multiaquifer system as a practical tool for the selection of sampling sites representative of the hydrogeochemical conditions of each aquifer along preferential flow paths and their proper individual characterization in terms of the natural and anthropogenic processes.

Hydrogeochemical evolution and mixing of groundwater
The PCA, Stiff, bivariable plots, and Piper diagrams in Figs. 3a, 6 and 7 and Figure S1 of the ESM showed a substantial spatial geochemical evolution in the studied aquifers. In LB2 and  a c b d f Fig. 9 Relations of a distance along the flow path vs. SI-tridymite, b SiO 2 vs. SI-quartz, c calcium vs. SI-calcite, d distance along the flow paths vs. Log pCO 2 , and e Log pCO 2 vs. calcium for groundwa-ter in the Barva and Colima multiaquifer system. Dry season refers to March survey, early wet season to June survey, and late wet season to January survey UC, the mineralization increases as groundwater travels from the recharge areas towards the lower zones of the basin, with longer flow paths in a NE-SW direction. LB1 increasing mineralization is observed if the UB samples are considered as pertaining to this flow path, representing its recharge chemical signal. In UC and LC, longer transit and residence times may explain higher contents in major ions, except for calcium. However, in LC, major constituents decrease as groundwater progressively approaches the lower zones. Although, several geochemical processes may cause some mineral species to decrease in ion content, the LC samples approach similar concentrations to the observed concentrations in LB and UC2 as the water approaches the discharge zone near (RN15), in terms of chloride, sodium, magnesium, and sulfate (Figs. 6c,d and 7b,c). This may suggest a groundwater mixing pattern from the overlying formations. In total, three potential groundwater mixing processes are observed in the PCA (Fig. 3a) and several hydrochemical scatter plots (Figs. 7e, f and 8a-d), as follows: 1. UB recharging water into LB, since samples from LB continued to follow the chemical evolution trend of UB. The PCA showed that the LB sampling sites at higher altitudes are adjacent to the samples from UB, in the negative section of F1 and F2, while the remaining samples are positioned in the positive section of F1 and F2, reflecting a spatial variation of groundwater hydrochemistry from the preferential recharge zone to lower zones in the basin. The vertical groundwater movement from UB towards underlying LB was also proposed by Protti (1986) and BGS/SENARA (1988). 2. LB2 probably percolating into UC2. The UC2 samples falling adjacent or among the LB2 samples in the F1 and F2 plot for the PCA and in most of the hydrochemical biplots suggest the mixing of groundwaters probably due to the vertical movement of water from overlying LB2. 3. LB2 and UC2 percolating into LC as groundwater approaches the discharge zone towards site RN15. Approaching the discharge zone, the hydrochemistry is gradually becoming more similar to that of UC2 and LB2 (Fig. 8a-d), evolving from Mg 2+ -HCO 3 towards Ca 2+ -Mg 2+ -HCO 3 facies. This indicates a potential flow from these overlying aquifers as the water flows in an E-SW direction. Although their permeability is low, groundwater may vertically move through the tuffs and ignimbrites of Tiribí and Puente Mulas formations (Arias-Salguero et al. 2006;BGS/SENARA 1988;Reynolds-Vargas and Fraile-Merino 2009).

Mineralization due to water-rock interactions
In the Barva and Colima aquifers, the dissolution of CO 2 might enhance the contents of bicarbonate, which is the predominant ion in the study area. Elevated CO 2 on Barva and Colima can originate from dissolution in the rain recharging the aquifer, CO 2 released in the soil during organic matter decomposition, and degassing due to active volcanism and hydrothermal activity, as observed in the Central Volcanic Cordillera (Aiuppa et al. 2014). Dissolution of CO 2 releases H 2 CO 3 and bicarbonate, according to Eqs. (8) and (9): The production of H 2 CO 3 may explain the decline in pH values with increasing CO 2 in Barva and Colima samples (Fig.  S4e of the ESM). Since H 2 CO 3 accelerates the weathering of aluminosilicate minerals and, thus, bicarbonate dissolution (Kim et al. 2020, Mechal et al. 2017, Mgbenu and Egbueri 2019, water-rock interactions may be responsible for most of the increasing contents of bicarbonate along the Barva and Upper Colima flow paths (Fig. 7a-c), and the increase of calcium with increasing CO 2 (Fig. 9e). However, the excess of bicarbonate with respect to the dissolution of anorthite (Fig. 8a), the predominant mineral in Barva and Lower Colima formations, may indicate additional sources. Although the incongruent dissolution of albite, hypersthene, and pyroxene release more bicarbonate per mol of calcium than that of anorthite, the proportion of these minerals in the Barva and Colima Formation is lower (0-9%, Table S1 of the ESM). Other possible sources related to anthropogenic activities may be further evaluated, e.g., the percolation of wastewaters, particularly due to the widespread use of septic tanks (AyA 2016).
The increase in bicarbonate and major cations along the LB and UC flowpath is mainly driven by the weathering of aluminosilicates, enhanced by the presence of CO 2 and the further production of H 2 CO 3 . The correlations found in the diagram of [calcium + magnesium + sodium + potassium] vs.
[bicarbonate] (Fig. 7e) indicate that bicarbonate is the balancing anion, suggesting a common origin related to the hydrolysis of silicate minerals, particularly in the upper basin. In volcanic aquifers in Ethiopia and Mexico, Haji et al. (2021), Mechal et al. (2017) and Morán-Ramírez et al. (2016) obtained similar correlations in the diagram of [calcium + magnesium + sodium + potassium] vs. [bicarbonate], which lead to the conclusion that silicate hydrolysis was the primary hydrochemical process responsible for the increase in bicarbonate and major cations along the flow paths. In Barva and Colima, the apparent importance of bicarbonate as a balancing anion decreases along the flow paths due to higher contents in nitrate and chloride, in the case of LB and UC2, and to a more elevated chloride and sulfate content in UC1 and LC, hindering the evidence of a common origin, e.g., silicate weathering, of bicarbonate and major cations.
The chloride vs. bicarbonate/chloride biplot (Fig. 7f) showed larger transit times in LB and UC as chloride increases and the relative importance of bicarbonate to chloride decreases. This general pattern of increasing salinity might reinforce the hypothesis of LB and locally unconfined UC recharging preferentially in the upper zone, where rainfall is more abundant. However, in LC, more significant transit and residence times are represented mainly by the chemical characteristics in RN09 since groundwater seems to experience a mixing process towards the discharge zone in RN15, probably due to percolation from overlying aquifers, causing the lowering of chloride concentration along the flow path.
The water of the shallow and perched UB aquifer exhibits the typical hydrochemical evolution of young groundwaters with slight but nonsignificant seasonal chemical variations. As it infiltrates into LB and travels through the assessed flow paths, the water exhibits hydrochemical characteristics corresponding to older groundwaters, with lesser seasonal variability, compared to UB. Higher calcium contents were found in samples in the discharge zone of the Barva aquifer compared to those in , probably due to the dissolution of the basaltic-andesitic and andesitic rocks in the Barva Formation. This observation coincides with high proportions of calcic feldspars, mostly anorthite (44-80%), reported by Rojas et al. (2017) and Arredondo- Li and Soto (2007), given that mafic or basic rocks, such as basalts, exhibit calcic rather than sodic feldspars and higher ferromagnesian mineral contents. However, additional calcium in LB may also be related to its percolation from limed agricultural fields or through wastewater infiltration due to the widespread use of septic tank drainage systems in the study area. The contribution of calcium from anthropogenic sources requires further examination.
The low concentrations of sulfate in LB samples may correspond to the reduced abundance of these ions in the silicate rocks (Custodio and Llamas 1983) and, thus, sulfide oxidation is a nonsignificant process. The nonaltitudinal distribution of samples with higher sulfate and the moderate correlation between sulfate and temperature may suggest a volcanism-related origin of the sulfate only in samples nearby the Barva crater (UB samples) and from deeper aquifers (UC2 and LC). Sulfur-rich mineral deposits may occur around fumaroles and in the proximity of deep veins in andesitic stratovolcanoes due to active volcanism (Zimbelman et al. 2005). However, these minerals are not described in the study area; moreover, the general oxidant conditions of groundwater, due to the elevated DO and ORP values, may prevent the occurrence of major redox reactions to altering the water chemistry in the studied aquifers.
The hydrogeochemical characteristics in UC1 and LC, with higher bicarbonate, chloride, magnesium, sulfate, and sodium (Figs. 5,6,Fig. S2 and Table S4 of the ESM), might broadly correspond to that commonly observed in mineralized groundwaters of longer residence time (Custodio and Llamas 1983). Furthermore, groundwater in UC1 showed the typical chemical composition of deep flow paths, with longer flow paths and lesser groundwater alteration by mixing processes, whereas UC2 flow path exhibited mixing with groundwater from the LB2. The LC showed deep and older groundwater characteristics solely at the easternmost sites (RN09, 37 and 47), while mixing with UC2 occurs as the groundwater moves westward.
In UC1 and LC, sodium, potassium, and bicarbonate enrichment might be the result of groundwater moving through acidic fractured volcanic rocks such as pumice, ignimbrite, ash, and trachyte (Mechal et al. 2017). Felsic or acid rocks contain more than 66% silica, potassium and sodiumrich feldspars and plagioclases, and low Fe-Mg contents. Since the lavas in Colima Formation, classified as andesites and trachyandesites, contain more alkaline feldspars and higher Na 2 O and K 2 O content than the Barva Formation (Kussmaul 1988), water-rock interactions may be responsible for the major contents in sodium and potassium in LC. Particularly, further enrichment in Colima might have also been obtained due to the infiltration of groundwater through the acidic ignimbrites, tuffs, and pumice-rich layers of the overlying Tiribí Formation (BGS/SENARA 1988). Tiribí lithologies are acid, based on their SiO 2 content (from 55.1 to 69.2 wt%; Hannah et al. 2002), and comprise andesites, trachyandesites, and dacites (Pérez et al. 2006) and elevated contents of Na 2 O and K 2 O (Kussmaul 1988). Higher magnesium concentrations in UC1 and LC, compared to LB, could be related to longer flow paths and residence times and consequent hydrolysis of the ferromagnesic minerals, i.e., hypersthene, augite, olivine. Lower calcium content in Colimas, even when the residence time of groundwater in the Colima aquifers is more significant than in LB (Sánchez-Murillo et al. 2022), and despite the presence of calcium-rich plagioclases (anorthite), might be the result of ion exchange processes (see section 'Ion-exchange processes') and lower vulnerability to anthropogenic sources.

Nitrate contamination
Mean values of nitrate in LB and UC2, i.e., 15.3 and 13.1 mg/L as NO 3 - (Table 5), indicate the alteration of the water quality in these aquifers. In addition, multivariate statistical analyses also showed that the LB and UC2 samples were mainly related to the nitrate and calcium variables, suggesting contamination due to a higher vulnerability to anthropogenic activities in the Barva aquifer and the interconnected UC2. Nitrate, commonly considered a minor ion, is the second major anion in LB and UC2, implying pollution from anthropogenic sources and possible groundwater infiltration from LB to UC2, as preliminarily observed by Madrigal-Solís et al. (2017), due to leaching of N-fertilizers and wastewaters from septic tank systems (Reynolds-Vargas et al. 2006). Groundwater infiltrating from LB2 to UC2 and UC2 to LC represents a pollution hazard for the deeper Colima aquifers, in terms of nitrate and other organic and inorganic pollutants. This topic is the subject of an expanded article on groundwater contamination.

Ion-exchange processes
This study shows normal and reverse ion exchange as plausible nondominant processes in the Barva and Colima aquifers. The results derived from estimated base exchange indices CAI-1 and CAI-2 (Fig. 8e) suggest normal ion exchange (calcium and magnesium adsorption by clays and potassium and sodium released into groundwater) near the preferential recharge zone progressively evolving into a reverse ion exchange process (potassium and sodium adsorbed and calcium and magnesium released) in LB and UC2. Most samples from UC1 and LC are located on the negative side, indicating normal ion exchange, which might partially explain the lower calcium concentration in these groundwaters. Negative CAI-1 and CAI-2 values were also reported in a volcanic aquifer by Haji et al. (2021), accompanied by higher sodium/calcium ratios in deeper groundwater due to cation exchange.
Furthermore, the relation of [Na + + K + -Cl -] vs. [Ca 2+ + Mg 2+ -SO 4 2-− HCO 3 -] showed that LB2, UC1, and LC were the only flow paths with slopes close to -1 (Fig. 8f), suggesting that cation exchange may influence their water quality substantially, according to the criteria by Fisher and Mullican III (1997). Therefore, particularly higher sodium contents and lower calcium in LC (Fig. 5) may result from calcium adsorption onto clay minerals and sodium released into groundwater; nevertheless, as the sodium contents in the LC are substantially higher, an additional amount of sodium is expected to derive from the weathering of sodium-rich silicate minerals. As observed in LC, Fenta et al. (2020) also found increasing sodium and decreasing calcium with depth in volcanic aquifers in Ethiopia. However, in contrast to LC, the slope observed by these authors in the relation of [Na + + K + -Cl -] vs. [Ca 2+ + Mg 2+ -SO 4 2-HCO 3 -] was far from -1, concluding that higher sodium was primarily due to the weathering silicate minerals and, to a lesser extent, to cation exchange. Conversely, in UB, LB1, UC2 flow paths, slopes deviated from -1 denote that cation exchange is not a dominant process, although it may occur to some extent.
In addition, most of the samples in LB and UC2 showed high calcium/bicarbonate molar ratios (>0.5), indicating the possible occurrence of a reverse ion exchange process or additional sources of calcium, according to Drever (1997) and Hounslow (1995). In contrast, low molar ratios in LC may indicate adsorption of calcium and magnesium due to normal ion exchange in 83% of the samples.
Although normal ion exchange might partly explain the significantly higher sodium and lower calcium contents in Colima, the samples also exhibited an excess of magnesium with respect to calcium along the dissolution line of augite and pyroxene (Fig. 8c), which could hinder the importance of the proposed ion exchange process. When assuming that normal ion exchange is involved, both calcium and magnesium would be adsorbed on the clay surfaces and, thus, low concentrations in groundwater would be expected. However, the lower holding capacity of magnesium by the negative charges on the clay's surface and its higher solubility in water, compared to that of calcium (Custodio and Llamas 1983) may cause calcium to be adsorbed to a greater extent.
Furthermore, all samples from Barva and Colima showed high Ca 2+/ (Ca 2+ + SO 4 2-) ratios (>0.5), indicating that the source of calcium could be related to plagioclase feldspars, as observed by Locsey et al. (2012) in basaltic groundwaters in Australia. However, only samples in UC1 (54%) and LC (92%) showed high Mg 2+ /(Mg 2+ + Ca 2+ ) ratios, reaffirming the larger importance of the weathering of ferromagnesian minerals in the Colima compared to Barva aquifer. Enrichment in magnesium with respect to calcium in volcanic aquifers due to the weathering of mafic minerals, such as augite, was also found by Edmunds et al. (2002) in Mexico City. Mechal et al. (2017) also found limited calcium in groundwaters interacting with acidic volcanic rocks due to the lack of plagioclase, olivine, and pyroxene, such as ignimbrite, pumice, tuff, ash rhyolite, and trachyte.

Saturation indexes
According to the SI, mineral precipitation may not explain the lower calcium contents in UC and LC, since no evidence for the precipitation of aragonite, calcite, or dolomite was found; most samples were undersaturated with respect to their solid and gas phases and only some samples approached equilibrium in Barva and Colima ( Fig. 9c and Fig. S4 of the ESM). In addition, relatively high CO 2 partial pressures (ranging from -2.97 to -1.47, Table S5 of the ESM) diminish the probability of calcium precipitation since precipitation of CaCO 3 is favored in natural waters with low dissolved CO 2 conditions (Liñán et al. 2021;Locsey et al. 2012). High CO 2 partial pressures release H 2 CO 3 through Eq. (8), shifting Eq. (10) to the right, increasing the dissolution of calcium-rich minerals. In contrast, low CO 2 partial pressures move Eq. (10) to the left, favoring calcite precipitation, as follows: Furthermore, low correlations between fluoride and calcium in all aquifers also revealed the low possibilities of fluorite precipitation and subsequent calcium precipitation. Thus, (10) CaCO 3 + H 2 CO 3 ↔ Ca +2 + 2HCO 3 − lower calcium concentrations in LC with respect to the dissolution of pyroxenes and plagioclases and compared to LB may be more related to normal ion exchange process and to possible lesser calcium contents in LC andesitic lithologies.

Conclusions
The present study provided the basis for understanding the natural and anthropogenic processes governing groundwater quality in the UB, LB, UC, and LC aquifers. This was achieved by (1) developing a discriminant model to classify samples according to the primary captured aquifer and flow path and (2) assessing the hydrochemical processes in each flow path.
The proposed discriminant model identifies the principal aquifer from which old and new samples are collected, enabling the selection of adequate sampling sites for the assessment of water quality and further in-depth investigations in four aquifers of a fractured multilayered volcanic system. The model minimizes the interpretation problems when using mixed water samples from various aquifers or undetermined aquifers, improving the present and future investigations into natural geochemical processes, pollution processes, recharge zones, residence times, length of flow paths, and subsequently, specific management plans for UB, LB, UC, and LC. The combination of PCA and DA in developing a discriminant model may apply to other multiaquifer systems of diverse lithological origins.
This study revealed that weathering of aluminosilicate minerals is the primary hydrogeochemical process governing groundwater quality through silicate hydrolysis in the Barva and Colima aquifers, whereas natural and anthropogenic processes distinctly modify the water chemistry in each aquifer and flow path. Moreover, hydrochemical and statistical analyses supported three vertical mixing processes: UB partially recharges the LB, the LB2 percolates towards the underlying UC2, and UC2 infiltrates into LC as groundwater travels westward.
In the Barva aquifers, the hydrochemical evolution can be described as follows (Fig. 10): (1) younger CO 2 -enriched groundwater in UB in the recharge zones releases H 2 CO 3 and bicarbonate, (2) as groundwater partially recharges the LB and travels towards lower zones, the acid accelerates the weathering of aluminosilicate minerals in andesitic-basaltic rocks, producing further groundwater enrichment in terms of bicarbonate, calcium, silica, and other major ions, (3) nitrate increase as groundwater travels underneath agricultural fields and urban areas lacking sewage systems due to the leaching of fertilizer and wastewaters.
Although the study has not covered the complete flow paths in the UC and LC aquifers, the evolution can be described as follows (Fig. 10): (1) in the recharge zones (when unconfined), groundwater is bicarbonate and H 2 CO 3 enriched due to dissolved CO 2 inputs from rainfall, organic matter decomposition, and dissolution of CO 2 from volcanic activity, (2) water-rock interactions enrich groundwater, especially in bicarbonate, silica, magnesium and, to a lesser extent, in calcium, (3) a large proportion of the calcium is adsorbed onto clays, (4) sodium and bicarbonate progressively increase due to the weathering of sodium-rich feldspars and the percolation through felsic volcanic rocks such as andesitic lavas, tuff, ignimbrite, ash, and pumice in Barva and Tiribí formations, while sodium is also released due to normal cation exchange, particularly in LC, and (e) potassium contents remain low even in the presence of potassiumrich plagioclases, since a portion may be irreversibly bound to clay minerals. The hydrochemical characterization along the UC1 flow path and that observed at the highest elevation site in LC is most consistent with that of typical, more mineralized, deep-circulating groundwaters with minimal seasonal chemical variations.

Recommendations and perspectives for future research
The enhanced knowledge of the hydrodynamic and hydrogeochemical behavior of individual aquifers and flow paths in a complex multilayer volcanic scenario may improve the delineation of management approaches and policies for protecting groundwater quality and quantity in the northwestern sector of the Central Valley. Shallow and moderately young groundwaters from Barva aquifers, exhibiting direct infiltration, a more significant seasonal chemical variability, and higher nitrate contents, present higher vulnerability to the leaching of contaminants of anthropogenic origin compared to the Colima aquifers. Given that the contamination of the Barva aquifers also affects the quality of the Colima aquifers, management strategies should prioritize the protection of Barva. Nevertheless, for the adequate management of each groundwater reservoir, future investigations should also focus on the hydrogeochemical processes of anthropogenic origin controlling water quality in both Barva and Colima aquifers individually.
Groundwater management should include general and specific protection strategies for each hydrogeological unit, considering that water quality in both Barva flow paths is distinctly affected by environmental pressures, the interaquifer connections, e.g., between LB2 and UC2, and between UC2 and LC in the southwestern section, and the subsequent percolation of contaminants from Barva to Colimas aquifers. In addition, to further improve the understanding of seasonal variability of the UB and LB, it is critical to study chemical and isotopic water quality for periods longer than 5 years, given the recurrence of climate-altering patterns such as El Niño-Southern Oscillation (ENSO). Finally, it is also recommended to include additional UC and LC monitoring sites, mainly to the north and east, to evaluate more extensive flow paths and recharge zones, which are poorly known for the Colima aquifers.