Geochemical connectivity among mine tailings, soils, sediments, and waters in a historical mining basin integrating multivariate statistics and exploratory random forest modeling

1 Introduction

River basins are dynamic geochemical systems in which rocks, sediments, soils, and waters are interconnected through weathering, erosion, transport, and deposition. These interactions regulate the mobilization and redistribution of metal (loid)s across the environment (Gaillardet et al., 2013). Water and solid matrices provide complementary geochemical information because they record different environmental processes: water chemistry reflects mobilization, flow paths, redox conditions, and seasonal variability, whereas soils and sediments preserve longer-term source signatures, depositional mixing, and historical redistribution (Domingo et al., 2023; Mukherjee et al., 2024). Source apportionment is particularly challenging in regions where natural geological enrichment and anthropogenic activities produce overlapping geochemical signatures, such as historical mining districts. In ore deposits, mineralized host rocks, mine wastes, enriched sediments, and soils frequently share elemental assemblages inherited from a common metallogenic origin. Furthermore, mine tailings and other historical residues remain exposed to weathering, allowing metal (loid)s to be continuously redistributed through surface runoff, fluvial transport, sediment remobilization, and water-rock interactions. Consequently, elevated elemental concentrations alone are insufficient to discriminate between geogenic and anthropogenic inputs, and may be misclassified as contamination, particularly in hydrothermal or volcanic terrains (Fornasaro et al., 2022). Moreover, source apportionment provides the geochemical context needed to interpret ecological and human health risks associated with metal (loid) exposure, because the dominant contributor to metal (loid) concentrations is not necessarily the dominant contributor to environmental risk. Instead, the relative importance of each source depends on the toxicity, mobility, bioavailability, and exposure pathways of the metal (loid)s involved (Lei et al., 2022; Kumari et al., 2025; Wang et al., 2025).

Developing integrated geochemical and statistical approaches to identify dominant geochemical sources and evaluate their relative contributions has therefore become a major research priority in environmental geochemistry (Kumari et al., 2024; Aránguiz-Acuña et al., 2025; Sims et al., 2025; Ferrari et al., 2026). Environmental studies increasingly combine geochemical characterization with multivariate statistical and machine learning approaches. Methods such as Principal Component Analysis (PCA) and Hierarchical Cluster Analysis (HCA) are widely used to identify elemental associations and common geochemical controls. In contrast, machine learning algorithms such as Random Forest (RF) offer complementary opportunities to explore source-related geochemical patterns in complex, multi-matrix datasets. RF is particularly advantageous because it can accommodate nonlinear relationships, highly correlated variables, and heterogeneous data structures while ranking the variables that most strongly influence classification performance (Breiman, 2001; Cao et al., 2023; Kumar et al., 2026). Despite these advantages, RF also has recognized limitations. Variable importance measures based on Mean Decrease in Impurity (MDI) may be biased when predictors are highly correlated or differ in their variability (Strobl et al., 2007). Moreover, RF is primarily a predictive classification algorithm rather than a mechanistic source-apportionment model, so variable importance should be interpreted cautiously, particularly in environmental systems where multiple processes influence geochemical signatures. These considerations highlight the importance of using RF to complement conventional multivariate analyses, together with appropriate model validation procedures, and geochemical interpretation (He et al., 2022; Ferrari et al., 2026).

The Guanajuato River Basin (GRB) represents a complex geochemical system where naturally mineralized lithologies coexist with centuries of intense mining and metallurgical activities within the Guanajuato Mining District (GMD). The district hosts low-sulfidation epithermal systems associated with the Veta Madre system and related structures that have been exploited for nearly five centuries (Buchanan, 1979). Historical records indicate that Guanajuato was a premier silver-producing region during the eighteenth century, generating substantial volumes of mining waste and profoundly altering the basin’s fluvial architecture (Antúnez-Echagaray, 1964; Ferry, 2019; Rueda-Garzon et al., 2024). This shared metallogenic origin creates substantial overlap between natural and mining-related geochemical signatures because the legacy mine wastes were derived directly from the same hydrothermal mineralization that naturally enriches the basin’s bedrock and surficial media. Previous studies in the GRB have substantially improved the understanding of environmental components, including mineralized bedrock, mine tailings, sediments, and water geochemistry (Carrillo-Chávez et al., 2003; Díaz-Delgado et al., 2009; Ramos-Arroyo et al., 2013; Chaparro et al., 2020; Lecomte et al., 2022; Loria et al., 2022). To address these challenges, the present study provides an integrated basin-scale assessment encompassing bedrock, mine tailings, soils, stream sediments, overbank sediments (fluvial terraces), surface waters, and groundwater. This approach combines geochemical characterization, multivariate statistical analyses, exploratory RF analyses, and ecological and human health risk assessment to provide a comprehensive evaluation of present-day metal (loid) distribution across the basin. The specific objectives were to: (i) characterize the hydrochemical and geochemical variability across the principal environmental matrices of the basin; (ii) evaluate elemental associations among environmental matrices using multivariate statistical analyses; (iii) complement the multivariate geochemical interpretation through exploratory Random Forest modelling to identify leading predictors associated with geogenic and mining-related signatures; and (iv) assess the environmental and human health risks associated with the current distribution of metal (loid)s. We hypothesize that environmental matrices in the GRB can be differentiated by their geochemical affinity to geogenic and mining-related signatures, with intermediate affinities reflecting mixed-source contributions, and that an integrated multi-matrix approach can identify these affinities and their implications for environmental contamination and risk.

2 Methodology

2.1 Study area and sample collection

The GRB is located in central Mexico and forms part of the Lerma-Chapala hydrological system. The basin drains the Sierra de Guanajuato toward the Silao-Romita aquifer and has been influenced by mining and metallurgical activities associated with the GMD (Figure 1). The geology of GRB is dominated by four principal lithostratigraphic assemblages: the Sierra de Guanajuato mafic-ultramafic complex, the Mesozoic volcano-sedimentary sequence, the Eocene Guanajuato Conglomerate, and the Cenozoic volcanic Group. These units record the tectono-magmatic evolution of the Guerrero Terrane and subsequent Cenozoic volcanism associated with the development of the GMD (Aranda-Gómez and McDowell, 1998; Nieto-Samaniego et al., 2005). Mineralization occurs predominantly in low-sulfidation epithermal Ag-Au vein systems controlled by the Veta Madre fault system and related NW-SE-trending structures. Hydrothermal activity generated mineral assemblages enriched in Ag, Au, Pb, Zn, Cu, As, Sb, and Hg, representing an important potential geogenic source of metal (loid)s in the GRB (Buchanan, 1979; Camprubí and Albinson, 2006).

A total of 71 samples representing different environmental matrices were collected throughout the GRB, including bedrock samples (RS, n = 12), mine tailings (RJ, n = 13), soils (RP, n = 5), stream sediments (RL, n = 15), overbank sediments (RT, n = 6), surface water (RD, n = 13), and groundwater (RG, n = 7). Bedrock samples were obtained from representative lithological units to characterize potential geogenic sources. The mine tailings samples were collected from historical deposits originally generated by amalgamation and cyanidation and subsequently reworked by bulk flotation processes (Ramos-Arroyo and Siebe-Grabach, 2006). Soils were sampled from profiles developed on the principal geological units of the GRB. Stream sediment samples were collected from depositional zones along the active channel of the main drainage. Overbank sediment samples were obtained from a ∼2 m-high floodplain terrace previously dated to pre-industrial depositional periods ranging from 2450 ± 40 yr BP to 9410 ± 40 yr BP (Miranda-Avilés et al., 2009). All solid samples were obtained using a stainless-steel shovel, stored in polyethylene bags, transported to the laboratory, air-dried, homogenized, and sieved to <2 mm before analysis.

Water samples were collected during the dry season (March 2025). Surface water samples (RD, n = 13) were obtained from three upstream reservoirs (Soledad PS, Esperanza PE, and Mata PM reservoirs) and one downstream reservoir (La Purísima PP). Groundwater samples (RG, n = 7) were collected from production wells in the lower basin. Electrical conductivity (EC), pH, temperature, and alkalinity were measured in situ. Total alkalinity was determined by acid titration using H₂SO₄ following the Hach method 8203, while pH, temperature, and EC were measured using the Hanna Instruments multiparameter meter. Water samples were filtered through a 0.45 μm MF-Millipore mixed cellulose ester (MCE) membrane filter. Aliquots for cation analysis were preserved with 2% ultrapure HNO₃ (v/v), whereas non-acidified samples were retained for anion determination. Water samples were stored in pre-cleaned high-density polyethylene (HDPE) bottles, transported at approximately 4 °C, and refrigerated until analysis.

2.2 Chemical analysis

All sample preparation and chemical analyses were conducted at the Environmental Geochemistry Laboratory of the Instituto de Geociencias, UNAM (Juriquilla, Mexico). Dissolved major anions in non-acidified water samples were determined via ion chromatography (IC, Dionex IC 1100 equipped with an AS14A-4 mm column). Dissolved metal (loid)s in acidified water samples and extractable elemental concentrations in solid samples were quantified by inductively coupled plasma optical emission spectroscopy (ICP–OES, Thermo iCAP 6500 Duo View).

Prior to analysis, solid samples underwent microwave-assisted acid extraction using a CEM MarsXpress system following a modified EPA 3051 protocol (4H1a1 method; Carrillo-Chávez et al., 2003). Quality assurance and control (QA/QC) procedures included processing procedural blanks, duplicate analyses, and certified quality control standards (QCS-26) in accordance with EPA Method 200.7 (Campisano et al., 2017). Elemental recoveries ranged from 93% to 109%, and analytical precision, expressed as coefficients of variation, was ≤5.42%.

2.3 Environmental quality assessment

A combination of pollution indices and human health risk metrics was used to evaluate environmental quality and potential impacts associated with metal (loid) occurrence. For the environmental quality indices, and for descriptive statistics, concentrations below the analytical limit of quantification (LOQ) were replaced by LOQ/2. Descriptive statistics, including mean, median, standard deviation, minimum, maximum, and Fisher’s skewness coefficient, were calculated to summarize the geochemical dataset and characterize elemental variability.

2.3.1 Water quality assessment

The Heavy Metal Pollution Index (HPI) proposed by Mohan et al. (1996) was used to assess the effects of dissolved metal (loid)s on water quality. HPI was calculated using Equation 1:

H P I = ∑ i = 1 n W i Q i ∑ i = 1 n W i (1)

Where the unit weight (W i) was defined as the inverse of the corresponding drinking water guideline value (WHO, 2022), and the sub-index (Q i) was calculated following Equation 2:

Q i = M i − I i S i − I i × 100 (2)

Where M i is the measured concentration of the i t h parameter, I i is the ideal value, assumed to be zero, and S i is the corresponding drinking-water guideline value (WHO, 2022). Higher HPI values indicate greater deterioration of water quality. According to Prasad and Bose (2001), HPI values >100 represent the critical pollution threshold.

The Heavy Metal Evaluation Index (HEI) was used to assess the overall degree of metal (loid)s contamination relative to drinking-water standards (Edet and Offiong, 2002). HEI was calculated following Equation 3:

H E I = ∑ i = 1 n M i S i (3)

Increasing HEI values indicate greater contamination intensity. Because no universally accepted HEI classification exists, contamination levels were classified using the relative approach proposed by Bhuiyan et al. (2010), whereby samples were classified as low (< mean HEI), moderate (≥ mean HEI and ≤2 × mean HEI), and high (>2 × mean HEI).

2.3.2 Exposure assessment

2.3.2.1 Non-carcinogenic health risk assessment

Potential human exposure to dissolved metal (loid)s was estimated using Average Daily Dose (ADD) calculations, considering ingestion and dermal exposure pathways (Equations 4 and 5), following methodologies established by the US EPA (2004), EPA (2011).

A D D i n g = C w a t e r × I R × E F × E D B w × A T (4)

A D D d e r m = C w a t e r × S A × K p × E T × E F × E D × 10 − 3 B W × A T (5)

Where Cwater is the concentration in water (mg/L), IR is the ingestion rate (L/day), EF is the exposure frequency (days/year), ED is the exposure duration (years), BW is the body weight (kg), AT is the averaging exposure time (days), SA is the exposed skin surface area (cm²), Kp is the dermal permeability coefficient (cm/h), ET is the exposure time (h/day). Exposure parameters were selected in accordance with the US EPA (2011) recommendations for chronic exposure scenarios (Supplementary Table S1).

Potential non-carcinogenic effects were evaluated using Hazard Quotient (HQ), Equation 6:where RfD is the reference dose (mg/kg*day), HQ values <1 indicate negligible risk, whereas HQ > 1 suggests potential for adverse health effects (US EPA, 2004). To assess cumulative non-carcinogenic effects from simultaneous exposure to multiple metal (loid)s, the Hazard Index (HI) was calculated according to Equation 7:HI values > 1 indicate potential combined non-carcinogenic health risk.

H Q = A D D R f D (6)

H I = ∑ H Q i (7)

2.3.2.2 Carcinogenic health risk assessment

Lifetime carcinogenic risk (CR) was estimated for As because it was detected in water samples and an oral carcinogenic slope factor was available for the evaluated exposure pathway, using Equation 8:where SF is the cancer slope factor (1.5; US EPA, 2002), carcinogenic risks were classified as negligible (<1.0E-06), acceptable (1.0E-06 to 1.0E-04), and potentially unacceptable (>1.0E-04) according to USEPA guidelines.

C R = A D D × S F (8)

2.3.3 Contamination and ecological risk assessment

The degree of metal (loid)s enrichment in mine tailings (RJ) was evaluated using the Contamination Factor (CF), with preindustrial overbank sediments (RT) representing the local preindustrial geochemical background. The individual Ecological Risk Factor (Er) and the Potential Ecological Risk Index (RI) were calculated according to Equations 9–11 proposed by Håkanson (1980): where T r i is the toxic response factor for element i, T r i adopted in this study, was 1 for Zn, 2 for Cr, 5 for Cu, Pb, and Ni, 10 for As, 30 for Cd, and 40 for Hg (Håkanson, 1980). C i is the measured concentration, and C b is the corresponding background concentration. According to Håkanson (1980) Er < 40 indicates low ecological risk, 40–80 moderate risk, 80–160 considerable risk, 160–320 high risk, and >320 very high ecological risk. For cumulative ecological risk, RI < 150 indicates low ecological risk; 150–300, moderate ecological risk; 300–600, considerable ecological risk; and >600, very high ecological risk.

R I = ∑ E r i (9)

E r i = T r i × C F i (10)

C f i = C i C b (11)

2.4 Multivariate statistical analysis and machine learning approaches

2.4.1 Quality control and data processing

Before multivariate statistical analyses, the raw dataset was reviewed for missing and censored observations. No true missing values were present in the dataset; all non-quantifiable analytical results were explicitly reported as below the analytical limit of quantification (<LOQ) and were therefore considered censored observations. Variables with >40% overall censored values were excluded from the dataset to minimize potential distortions in the multivariate structure (Helsel, 2012). For the retained variables, observations below the LOQ were treated using median imputation according to the corresponding environmental matrix. Geochemical concentrations were subsequently log10-transformed and standardized using Z-score to reduce skewness and ensure comparability among variables prior to multivariate statistical analyses (Reimann et al., 2008).

To evaluate whether the multivariate results were sensitive to the treatment of < LOQ observations, an additional sensitivity analysis was performed using three alternative censored data treatments: the median imputation approach used in the main analysis, LOQ/2 substitution, and multiple imputation by chained equations (MICE) as a multivariate iterative imputation approach. For each treatment, the resulting dataset was log10-transformed and Z-score standardized, and PCA and Spearman correlation analyses were repeated using the same analytical procedures. Robustness was evaluated by comparing the variance explained by the first two principal components, the main PCA loading patterns, and the overall Spearman correlation structure among treatments. Similarity between Spearman correlation matrices was quantified using Pearson correlation coefficients between the corresponding pairwise Spearman coefficients.

2.4.2 Multivariate statistical analysis

Elemental associations were evaluated using Spearman rank correlation coefficients (Spearman, 1904). The correlation analysis included 18 variables for the water dataset, corresponding to 153 unique pairwise tests, and 19 variables for the solid-matrix dataset, corresponding to 171 unique pairwise tests. To control multiple comparisons, p-values were adjusted using the Benjamini–Hochberg false discovery rate (FDR) procedure (Benjamini and Hochberg, 1995), with statistical significance defined as q < 0.05. Principal Component Analysis (PCA) was applied separately to a solid matrix and water dataset to identify dominant geochemical controls and elemental associations. Additionally, Hierarchical Cluster Analysis (HCA) was used to evaluate similarities among samples and environmental matrices using Euclidean distance and Ward’s minimum variance method (Ward, 1963; Jolliffe and Cadima, 2016).

2.4.3 Exploratory random forest analysis

A supervised Random Forest (RF) classifier was implemented in Python 3.10 using the scikit-learn library (version 1.2.2; Pedregosa et al., 2011) as an exploratory tool to evaluate geochemical affinity between environmental samples and geogenic and mining-related reference materials (Breiman, 2001). The model was trained using bedrock samples (RS, n = 12) to represent the geogenic reference class and mine tailings samples (RJ, n = 13) to represent the mining-related reference class (Supplementary Table S2). The RF model was implemented using 100 decision trees (n_estimators = 100) and a fixed random seed (random_state = 42) to ensure reproducibility. Complementary validation procedures, including the out-of-bag (OOB) error estimate and repeated stratified 10-fold cross-validation (five repetitions), were used to assess internal classification stability and model performance.

Variable importance was evaluated using Mean Decrease in Impurity (MDI), based on the Gini criterion, complemented by permutation importance (Supplementary Table S3). The trained model was subsequently applied to classify RP, RL, and RT samples to evaluate their relative affinity for each reference class. Model outputs were expressed as RF affinity scores. These scores were interpreted as an indicator of relative geochemical affinity rather than as definitive quantitative source contributions.

Sample-level uncertainty in RF affinity scores was evaluated using 1000 bootstrap iterations. In each iteration, the original training dataset comprising RS and RJ samples was resampled with replacement, and an RF classifier was fitted to the resulting bootstrap sample. The fitted model was then used to estimate the probability of affinity to the tailings reference class for each sample, generating a distribution of 1000 probability estimates per sample. Median probability scores and the 2.5th and 97.fifth percentiles of these distributions were used to define empirical 95% uncertainty intervals. The bootstrap resampling seed varied among iterations (random_state = i), whereas the RF initialization seed was held fixed (random_state = 42). Thus, variability among repeated predictions reflected changes in the composition of the bootstrap training samples. Samples with uncertainty intervals entirely above or below the 0.5 probability threshold were assigned to the Tailings or Rock reference class, respectively, whereas samples with intervals spanning 0.5 were considered uncertain/unassigned (Supplementary Table S4). Because RP, RL, and RT may integrate overlapping geogenic and mining related, their RF scores were interpreted as exploratory indicators of relative affinity within the continuum defined by the two reference classes.

3 Results

3.1 Hydrochemical and physicochemical characteristics

Due to limited hydrochemical variability observed within individual reservoirs, surface waters results are presented in four groups: PS, PE, PM, and PP. Descriptive statistics of physicochemical parameters, major ions, and selected metal (loid)s in RD and RG waters were summarized in Table 1. Water pH ranged from 6.33 to 8.22. RD waters exhibited higher mineralization than RG waters, as reflected by greater mean EC and TDS values (Table 1). RD waters exhibited higher mean Ca²⁺ (41.3 mg/L), Mg²⁺ (18.4 mg/L), and SO₄ ^2- (62.0 mg/L) concentrations than RG waters, whereas RG waters showed relatively higher Na⁺ (46.4 mg/L) and Si (27.2 mg/L) concentrations. The cation sequence was Ca²⁺ > Na⁺ > Mg²⁺ > K⁺ in RD waters and Na⁺ > Ca²⁺ > K⁺ > Mg²⁺ in RG waters. HCO₃ ⁻ was the dominant anion in both groups, with mean concentrations of 167 and 176 mg/L for RD and RG waters, respectively. The anion sequence was HCO₃ ⁻ > SO₄ ^2- > Cl⁻ > NO₃ ⁻ in RD waters, and HCO₃ ⁻ > SO₄ ^2- > NO₃ ⁻ > Cl⁻ in RG waters.

TABLE 1. Descriptive statistics of physicochemical parameters, major ions, and selected trace elements in surface water (RD) and groundwater (RG) samples from the Guanajuato River Basin.

Values reported as minimum (Min), maximum (Max), mean, median, standard deviation (SD). Drinking water guidelines and regulatory values from the World Health Organization (WHO, 2022), the United States Environmental Protection Agency (US EPA, 2009), and Mexican regulations (NOM-127-SSA1-2021; Secretaría de Salud, 2022) are provided for comparison. * Values indicate LOQ/2 substitutions for concentrations below the analytical limit of quantification (LOQ) used in the calculation of descriptive statistics.

Most metal (loid) concentrations were below the LOQ in both RD and RG waters. Among the detected analytes, As, Fe, Mn, Ni, and Zn exhibited the highest detection frequencies. As was quantified in 30.8% of RD samples and 57.1% of RG samples; corresponding values were 84.6% and 28.6% for Fe, 100% and 28.6% for Mn, 61.5% and 28.6% for Ni, and 53.8% and 100% for Zn, respectively. Overall, Fe, Mn, and Ni were more frequently detected in RD samples, whereas As and Zn showed higher detection frequencies in RG waters. Most samples complied with drinking-water guidelines established by the US EPA (2009), Mexican regulations (NOM-127-SSA1-2021; Secretaría de Salud, 2022), and the WHO (2022). However, As concentrations above the corresponding guideline value occurred in both water groups and were more pronounced in RG waters, whereas elevated Mn and Fe concentrations were primarily observed in RD waters.

The Piper diagram (Figure 2) showed clear hydrochemical differences between RD and RG waters. Upstream reservoirs (PS, PE, and PM) were characterized by Ca-Mg dominated compositions with a contribution of SO₄ ^2-, whereas RG waters showed greater Na + K and HCO₃ ⁻. PP waters, located downstream, exhibited an intermediate hydrochemical composition between the upstream reservoirs and RG waters. Although HCO₃ ⁻ was the dominant anion in all samples, RD waters showed a greater contribution of SO₄ ^2- than RG waters.

3.2 Geochemical signatures and source apportionment

3.2.1 Multivariate statistical analysis of water chemistry

To further evaluate the relationship among physicochemical parameters and dissolved constituents, Spearman correlation analysis and PCA were applied to the hydrochemical dataset. Al, As and Ni, which exhibited >40% overall censored observations, were excluded from the PCA and Spearman correlation analyses (Supplementary Table S5). The first two principal components explained 75.8% of the total variance, with PC1 and PC2 accounting for 47.8% and 28.0%, respectively (Figure 3). PC1 primarily separated most RG waters, associated with Si and F⁻, from RD samples associated with EC, TDS, Ca²⁺, Mg²⁺, SO₄ ^2-, HCO₃ ⁻, and Cl⁻. PC2 further differentiated RD samples among the four reservoir groups (PS, PE, PM, and PP). PP samples were primarily associated with HCO₃ ⁻, Cl⁻, and K⁺, whereas PE and PM showed strong association with Fe and Mn. PS samples were primarily associated with EC, TDS, Mg²⁺, and SO₄ ^2-.

Spearman correlation analysis further quantified the relationships among the hydrochemical variables (Figure 4). Strong positive correlations were observed between SO₄ ^2- and Mg²⁺ (ρ = 0.96; q < 0.001), F⁻ and Si (ρ = 0.85; q < 0.001), EC and Mg²⁺ (ρ = 0.84; q < 0.001), EC and SO₄ ^2- (ρ = 0.80; q < 0.001), Ca²⁺ and Mg²⁺ (ρ = 0.79; q < 0.001), and Ca²⁺ and SO₄ ^2- (ρ = 0.77; q < 0.001). Positive correlations were also observed between Fe and Mn (ρ = 0.70; q < 0.01) and F⁻ and Na⁺ (ρ = 0.69; q < 0.01). In contrast, Zn showed strong negative correlations with EC, Ca²⁺, Mg²⁺, and SO₄ ^2- (ρ = −0.76 to −0.90; q < 0.001).

3.2.2 Solid matrices, geochemistry, and multivariate analysis

Major element concentrations varied among solid matrices, with Al, Fe, Ca, Mg, and K occurring at relatively high concentration (Table 2). RJ samples were characterized by elevated Ca, Mn, and S concentrations, whereas RP samples showed higher Al and Fe concentrations. RL exhibited the highest median Mg and Na concentrations, while RT and RS samples displayed intermediate major-element compositions.

TABLE 2. Major element concentrations in solid matrices from the Guanajuato River Basin.

Values are reported as median (minimum-maximum) concentrations for mine tailings (RJ), soils (RP), stream sediments (RL), overbank sediments (RT), and bedrock (RS). Upper Continental Crust (UCC) values from Rudnick and Gao (2013) are included as a general geochemical reference.* Values indicate LOQ/2 substitutions for concentrations below the analytical limit of quantification (LOQ) used in the calculation of descriptive statistics.

Metal (loid)s distribution also varied among matrices (Table 3). RJ samples were enriched in Cu, Cd, Sb, and Pb, whereas RT samples showed the highest median concentration of Cr, V, and Ni. RL and RP samples exhibited relatively elevated Zn and As concentrations, respectively. RS samples displayed wide concentration ranges, reflecting the heterogeneous composition of the lithological units. Pb exhibited the widest concentration range, reaching a maximum of 13402 mg/kg in RJ samples, and 160 mg/kg in RL samples.

TABLE 3. Trace element concentrations in solid matrices from the Guanajuato River Basin.

Values are reported as median (minimum-maximum) concentrations for mine tailings (RJ), soils (RP), stream sediments (RL), overbank sediments (RT), and bedrock (RS). Upper Continental Crust (UCC) values from Rudnick and Gao (2013) are included as a general geochemical reference.* Values indicate LOQ/2 substitutions for concentrations below the analytical limit of quantification (LOQ) used in the calculation of descriptive statistics. < LOQ: indicates matrices in which all samples were below the analytical limit of quantification.

Spearman correlation analysis quantified the relationships among elemental concentrations (Figure 5; Supplementary Table S6). The first association comprised Al, Fe, Cr, V, and Ni, which showed strong positive correlations (ρ = 0.67 to 0.87; q < 0.001). A second association involved As, Cd, Cu, Pb, and Zn, including strong positive correlations between Cd-Cu (ρ = 0.75; q < 0.001), As-Cd (ρ = 0.71; q < 0.001), Cu-Zn (ρ = 0.70; q < 0.001), and Pb-Cd (ρ = 0.61; q < 0.001).

PCA explained 52.7% of the total variance, with the first two principal components accounting for 35.8% and 16.9%, respectively (Figure 6). PC2 distinguished two contrasting elemental associations. Positive PC2 scores were associated with Al, Fe, V, Cr, and Ni, whereas negative PC2 scores were primarily associated with Pb, Cd, Zn, and As. Most RJ samples clustered in the positive PC1 and negative PC2 region, aligned with the S, Pb, Cd, As, and Zn vectors. RT samples clustered predominantly toward Al, Fe, V, Cr, and Ni association, whereas RP showed a more dispersed distribution. RL samples occupied intermediate positions within the ordination space and showed relatively broad dispersion. RS samples exhibited the widest distribution across the PCA ordination space.

A sensitivity analysis using alternative treatments of censored (<LOQ) observations showed that the main multivariate patterns were stable across LOQ/2 substitution, median imputation, and iterative imputation. The variance explained by PC1 was 35.7%, 35.8%, and 35.9% respectively, while PC2 explained 17.0%, 16.9%, and 16.4%. The Spearman correlation structures also showed high agreement among treatments, with Pearson correlation of r = 0.997 for LOQ/2 vs. median imputation, r = 0.981 for LOQ/2 vs. iterative imputer, and r = 0.993 for median imputation vs. iterative imputer (all p < 0.001). Consistent patterns were also observed in the PCA biplots and Spearman correlation heatmaps across treatments, indicating that the main multivariate structure was not substantially affected by the treatment of censored observations (Supplementary Table S7; Supplementary Figure S1).

The HCA dendrogram showed several groupings among the analyzed matrices (Figure 7). RT and RP samples clustered within a relatively restricted portion of the dendrogram, whereas RL and RS samples were distributed among multiple subclusters, reflecting greater geochemical variability. RJ samples were mainly concentrated in the central portion of the dendrogram but showed variable associations with RT, RP, and RL samples. RS 06, RS 10, and RS 12 were among the most dissimilar bedrock samples and were separated from the other samples at relatively high Euclidean distances.

3.2.3 Random forest-based geochemical affinity assessment

As a complementary exploratory analysis, the RF classifier was trained using RS samples as the geogenic reference class and RJ samples as the mining-related reference class. The RF classifier achieved an OOB score of 0.88 (OOB error = 0.12), and a mean repeated cross-validation accuracy of 0.880 ± 0.098. Sensitivity analysis showed generally stable model performance across the evaluated hyperparameter settings. Maximum tree depth produced no appreciable changes in OOB score or cross-validation accuracy within each evaluated number tree, whereas models with 100 and 200 trees showed the highest mean cross-validation accuracy (Supplementary Table S3).

Variable importance based on the MDI identified S as a leading predictor for discrimination between the two reference classes, followed by Cd and As (Figure 8). Permutation importance further identified S as the strongest individual predictor (Supplementary Table S3). At the matrix level, median RF affinity scores for RT and RP samples were shifted toward the geogenic reference class (61% and 63%, respectively), whereas RL samples occupied an intermediate position, with median scores of 49% geogenic and 51% mining related affinity (Table 4). RL samples also showed the broadest range of RF affinity scores, indicating greater variability in their relative position between the two reference classes. Sample level median FR affinity scores and their empirical 95% uncertainty intervals derived from 1000 bootstrap model fits showed variable uncertainty among individual samples, particularly within the intermediate environmental matrices (Supplementary Table S4).

TABLE 4. Summary of Random Forest affinity scores (%) describing relative affinity to geogenic (RS) and mining related (RJ) reference classes across solid matrices of the Guanajuato River Basin.

Values represent the median RF, affinity scores and observed sample range within each environmental matrix. Scores indicate relative affinity to the geogenic (RS) and mining-related (RJ) reference classes and should not be interpreted as prediction uncertainty. Sample level median scores and 95% uncertainty intervals derived from repeated model fits are provided in Supplementary Table S4.

Spatially, RL samples located near the historical mining district generally showed lower median affinity scores and coincided with higher values of the standardized S + Cd + As index (Figure 9). Elevated index values also occurred at some downstream locations, indicating that this geochemical signature was not restricted to the immediate vicinity of historical tailings deposits.

3.3 Environmental assessment and human risk

3.3.1 Water quality assessment

HPI values ranged from 56.9 to 279 in surface waters and from 56.5 to 350 in groundwater samples (Table 5). According to the criterion proposed by Prasad and Bose (2001), surface water samples RD 04, RD 05, RD 06, RD 07, RD 08, RD 09, RD 10, RD 11, and RD 13, as well as groundwater samples RG 02, RG 03, RG 06, and RG 07, exceed the critical HPI value of 100. The highest HPI values were recorded in RG 06 (350), RG 07 (293), and RD 04 (279), indicating the greatest degree of metal (loid) related water quality impairment among the samples analyzed. The highest HEI values were also observed in RD 04 (23.3), followed by RD 07 (10.0), RD 05 (9.28), and RD 06 (8.08), while in groundwater, the highest values corresponded to RG 06 (4.51) and RG 07 (4.14).

TABLE 5. Water quality and human health risk indices for surface water (RD) and groundwater (RG) samples from the Guanajuato River Basin.

HPI: Heavy Metal Pollution Index, HEI: Heavy Metal Evaluation Index, HI: Hazard Index, CR (As): Carcinogenic Risk for As. HPI and HEI were calculated using dissolved As, Mn, Ni, and Zn concentrations. HI and CR were calculated following the US EPA (2011) methodology. *As concentrations below the analytical limit of quantification (LOQ=0.0143 mg/L), LOQ/2 was used to estimate CR values. These values were not interpreted as confirmed carcinogenic risk exceedances. Bold values indicate risk estimates based on quantified concentrations (≥ LOQ), which were included in the risk assessment.

3.3.2 Ecological risk assessment

Mine tailings represented the principal ecological risk reservoir within the GRB (Table 6). RI values ranged from 24.2 (RJ 03 days) to 299 (RJ 05 b). According to Håkanson (1980), RJ-1, RJ-2, RJ-3, and RJ-6 were classified as low ecological risk deposits, whereas RJ 05 and RJ 07 approached the moderate risk category. Individual ecological risk factors identified As and Cd were important contributors, whereas Cu dominated RI of RJ 05 and RJ 07. Cr, Ni, and Zn generally showed comparatively low individual ecological risk factors.

TABLE 6. Contamination factor, ecological risk factor, and potential ecological risk index for historical tailings deposits (RJ) from the Guanajuato River Basin.

CF: Contamination factor, Er: Ecological risk factor, RI: Potential ecological risk index. Indices were calculated according to (Håkanson, 1980).

3.3.3 Human health risk assessment

Non-carcinogenic risk assessment indicated that most RD and RG samples remained within acceptable limits (Table 5). HI values ranged from 0.312 to 0.687 in RD waters and from 0.311 to 1.93 in RG waters. Samples RG03, RG06, and RG07 exceeded the reference threshold of 1, indicating the potential non-carcinogenic health effects under the evaluated long-term exposure scenario.

Carcinogenic risk assessment focused on As (Table 5). As was quantified in 40% of the water samples, whereas 60% of the measurements were below the analytical LOQ. Therefore, the CR results require careful interpretation. For samples with As concentrations below the LOQ, CR values were estimated using the midpoint (LOQ/2) substitution approach and were considered screening estimates only (Supplementary Table S8). Accordingly, CR was interpreted only for samples with quantified As concentrations. For samples with quantified As concentrations, CR values ranged from 2.54E-04 to 7.84E-04 (Table 5). The highest values were observed in groundwater samples RG 06 and RG 07, followed by RG 03 and RG 02, while elevated CR values were also identified in downstream surface water samples (RD 09-RD 11 and RD 13). According to US EPA (2004) guidelines, these values exceeded the upper acceptable carcinogenic risk level (1.0E-04), indicating potential carcinogenic concern under the evaluated long-term exposure scenario.

4 Discussion

4.1 Geogenic and mining related controls on trace element distribution

Maximum concentrations of As, Cd, Cu, Pb, and Zn were observed in individual RS samples, consistent with the known low sulfidation epithermal mineralization of the GRB. Buchanan (1979) described the GMD as a fossil geothermal system characterized by hydrothermal assemblages containing sphalerite (ZnS), galena (PbS), chalcopyrite (CuFeS₂), pyrite (FeS₂), Cd-bearing sulfides, and mineral phases associated with As and Hg, indicating that these elements are intrinsic components of the geological framework of the GRB. The occurrence of similar metal (loid)s in preindustrial RT samples further demonstrates that these elements were present prior to intensive mining, indicating that present day distribution of trace elements within the GRB cannot be attributed exclusively to mining activities but also reflects the natural geochemical signature of the mineralized bedrock (Miranda-Avilés et al., 2009; Ramos-Arroyo et al., 2013; Rueda-Garzon et al., 2022; Ponce-Hernández et al., 2025).

RT samples provide a valuable reference for evaluating geochemical conditions before large-scale mining and metallurgical activity. The As, Cd, Cu, and Zn concentrations in RT samples indicate that weathering of naturally mineralized lithologies has supplied trace elements to the fluvial system before the onset of intensive mining. However, the contrasting geochemical compositions of RJ and RL samples suggest that historical mining has likely modified their distribution. Compared with RT samples, RJ exhibited higher concentrations of ore-related metal (loid)s, whereas Cr, Ni, and V remained primarily associated with natural lithological sources. Together with the consistent geochemical associations identified by PCA and the partial overlap among solid matrices revealed by HCA (Figures 6, 7), these results suggest that historical mining primarily enhanced the abundance and redistribution of hydrothermal related elements already present in the basin rather than generating distinct elemental associations.

The exploratory RF analysis provides additional evidence that RL samples exhibit mixed geochemical affinity rather than a clear association with either geogenic or mining-related reference classes. The nearly equal median RF affinity scores for geogenic (49%) and mining-related (51%) reference classes, together with their wide ranges among samples (18%–81% and 19%–82%, respectively), indicate substantial variability among individual RL samples and mixed geochemical signatures rather than definitive source attribution (Table 4). This mixed affinity is consistent with the intermediate geochemical characteristics of RL relative to RS and RJ, suggesting that stream sediments integrate both natural and mining-derived materials. Such patterns are consistent with the long-term influence of weathering, erosion, and sediment transport processes operating throughout the GRB. Historical records indicate that large volumes of finely ground wastes were deposited in tailings deposits or discharged directly into drainage channels during ore processing (Antúnez-Echagaray, 1964; Ramos-Arroyo et al., 2004). Because these materials are more susceptible to weathering and erosion than the original bedrock, they serve as an efficient mechanism for the downstream redistribution of metal (loid)s and contribute to the persistence of mining-related signatures in modern sediments (Lottermoser, 2010; García-Arreola et al., 2018, 2025).

Together, these observations indicate that the present-day geochemical configuration of the GRB reflects the combined influence of geological inheritance and historical mining activities. The overlap between mineralized bedrocks and historical mine waste makes source apportionment inherently challenging when based solely on elemental concentrations. Therefore, the RF affinity sources should be interpreted as indicators of relative geochemical affinity rather than definitive source discrimination.

4.2 Geochemical connectivity among environmental matrices

Despite the elevated concentrations of As, Cd, Cu, Pb, and Zn observed in solid matrices, most dissolved metal (loid)s remained below the LOQ in both RD and RG waters. This contrast suggests that the processes controlling metal (loid) accumulation in solid matrices differ from those governing their dissolved concentrations in water. PCA results showed that water samples were primarily differentiated by major ion chemistry rather than by trace metal (loid) concentrations (Figure 3). RD and RG waters exhibited distinct hydrochemical signatures despite occurring within the same basin. Similar hydrochemical patterns have been reported previously for the GRB and have been attributed to local water-rock interactions, residence time, and mixing processes rather than to direct dissolution of metal (loid) rich solid phases (Carrillo-Chávez et al., 2003; Lecomte et al., 2022). Accordingly, the present results characterize the dominant geochemical conditions observed during the low-flow period.

However, low dissolved concentrations do not imply the absence of geochemical connectivity among environmental matrices. Previous studies have shown that the carbonate-rich composition of Guanajuato mine tailings buffers acidity and limits acid mine drainage generation, favoring metal (loid) retention through interactions with Fe-Mn oxyhydroxides and carbonate phases (Carrillo-Chávez et al., 2003; García-Arreola et al., 2018). Under these conditions, metal (loid)s may remain preferentially associated with solid matrices and be transported primarily in particulate rather than dissolved form. As represents an important exception because its mobility may increase under circumneutral to alkaline conditions through processes such as desorption from mineral surfaces and other water-mineral interactions (García-Arreola et al., 2025). Consequently, As remains one of the few elements occurring at measurable dissolved concentrations while also showing enrichment in solid matrices, highlighting its relevance for evaluating geochemical connectivity between solid and aqueous matrices within the basin.

4.3 Environmental and human health implications

The ecological risk assessment revealed marked spatial heterogeneity among historical mine tailings deposits, indicating that mining residues do not represent a uniform environmental concern across the GRB. RJ 05 and RJ 07 deposits exhibited the highest RI values and were classified as moderate ecological risk, whereas the remaining deposits were classified as low ecological risk. Elevated RI values of RJ 05 and RJ 07 were primarily associated with As, Cd, and Cu contributing. These differences suggest that the environmental significance of individual deposits may depend on their specific geochemical composition, mineralogical characteristics, weathering history, and physical stability (Miranda-Avilés et al., 2025).

The standardized S + Cd + As index in RL samples showed higher values in the central-northern sector of the basin, corresponding to the historical mining district and coinciding with the location of several mine tailings deposits (Figure 9). Elevated index values also occurred downstream of the mining district, indicating that the geochemical signature is not confined to the primary mining area. Previous investigations have reported moderate to very high ecological risk in RL and RP samples, particularly associated with enrichment in As, Cd, Pb, and Zn (Rueda-Garzon et al., 2022; Ponce-Hernández et al., 2025). Together, these observations suggest that historical mining has contributed to localized metal (loid) enrichment and downstream redistribution of metal (loid)s through sediment transport and depositional processes (Banerjee et al., 2023).

Regarding water quality, HPI values exceeded the critical threshold in several surface water and groundwater samples. Among surface waters, particularly high HPI and HEI values coincided with elevated Mn concentrations, especially at RD 04 and RD 07, whereas high index values in groundwater samples RG 06, and RG 07 coincided with As concentrations exceeding drinking-water guideline values (Table 1 and 5). Similar hydrochemical conditions have been reported previously for the GMD, where circumneutral to alkaline conditions and carbonate-rich materials may limit the mobility of many potentially toxic metal (loid)s and their occurrence in dissolved form (Carrillo-Chávez et al., 2003; García-Arreola et al., 2018; Hernández-Mendiola et al., 2022; Rueda-Garzon et al., 2025).

The human health assessment revealed greater concern for groundwater than for surface water. All RD samples exhibited HI values below the accepted threshold of 1, indicating limited non-carcinogenic risk, whereas RG 03, RG 06, and RG 07 exceeded this threshold, indicating potential non-carcinogenic health effects under evaluated long-term exposure scenario. The higher risk associated with these RG samples is consistent with the elevated As concentrations observed in the lower basin and agrees with previous reports describing the occurrence of As in groundwater systems (Carrillo-Chávez et al., 2003). For samples with quantified As, carcinogenic risk estimates exceeded the upper acceptable risk level (1.0E-04), with the highest values observed in RG 06 and RG 07. Although these estimates indicate potential carcinogenic concern under the evaluated exposure scenario, their interpretation should consider the high proportion of As measurements below the LOQ and the screening level nature of the risk assessment.

The contrast between ecological and human health risk assessments indicates that environmental concerns in the GRB are spatially heterogeneous, with ecological risk concentrated in specific mine tailings deposits and human health concerns primarily associated with As in groundwater. From a management perspective, the integrated assessment provides a screening framework for prioritizing monitoring and further investigation. Based on the present findings, RG 05 and RG 07 represent priority mine tailings deposits for monitoring because they exhibited the highest ecological risk, whereas groundwater sites RG 06 and RG 07 require particular attention because of their elevated As-related health risk. Downstream RL sediments may also provide useful indicators for tracking the transport and redistribution of mining-related geochemical signatures within the basin. Continued monitoring across contrasting hydrological conditions would help evaluate temporal changes and potential impacts on water resources used for human consumption.

4.4 Study limitations and implications

The interpretation of the geochemical patterns identified in the GRB should consider the characteristics of the sampling design and analytical approach. Sample sizes varied among environmental matrices, particularly for groundwater, and the RS and RJ samples used to train the RF classifier represent a relatively limited training dataset. This introduces additional uncertainty in the estimated affinities of intermediate matrices. In addition, because the RF classifier was trained using lithological samples and mine tailings as reference classes, the observed discrimination may partly reflect differences in matrix characteristics, weathering history, and ore processing effects in addition to source-related geochemical signatures. Therefore, the RF probabilities are interpreted as indicators of relative geochemical affinity rather than definitive source attribution. Complementary mineralogical and isotope-based investigations could further constrain source relationships in mineralized catchments such as the GRB (Wiederhold, 2015; Desaulty and Petelet-Giraud, 2020; Ding et al., 2024).

Water sampling was restricted to the dry season; therefore, the hydrochemical patterns identified in this study represent low-flow conditions and may not capture seasonal variations in metal (loid) mobility and transport. This temporal constraint limits extrapolation of the observed geochemical relationships to contrasting hydrological conditions, highlighting the need for seasonal monitoring. In addition, Al, As, and Ni were excluded from the water PCA and Spearman correlation analyses because they exceeded the established censoring threshold. This limitation is particularly relevant for As because of its environmental and human health significance in the GRB. For retained variables, the treatment of censored observations by median imputation may introduce some uncertainty into multivariate relationships, particularly for variables with censoring levels approaching the established threshold. However, the sensitivity analysis using alternative censored data treatments showed highly consistent PCA and Spearman correlation structures, indicating that the main multivariate patterns were robust to the treatment of censored observations. For As, CR values based on quantified concentrations were distinguished from estimates derived from concentrations below the LOQ, which depended on the LOQ/2 substitution assumption and were therefore considered screening estimates.

The environmental and human health risk indices should be interpreted as screening-level assessments under the evaluated conditions. Site-specific factors such as contaminant bioavailability, local exposure conditions, and seasonal variability may influence actual ecological and human health risks (Lei et al., 2022; US EPA, 2004; 2011). Finally, the dataset and analysis code are not currently available through a public repository, which limits independent reproduction of the statistical and machine learning analyses.

5 Conclusion

The Guanajuato River Basin’s current geochemistry reflects both hydrothermal mineralization and historical mining. Groundwater and surface waters generally exhibited low dissolved metal (loid) concentrations, although localized As and Mn enrichment affected water quality, with As representing the principal human health concern. Historical mining has amplified and redistributed preexisting hydrothermal geochemical signatures rather than generating distinct elemental associations. Stream sediments integrate geogenic and mining related geochemical signals across the basin. Multivariate analyses identified two principal geochemical associations reflecting geogenic (Al, Fe, Cr, V, and Ni) and mining related signatures (As, Cd, Cu, Pb, and Zn), while the exploratory Random Forest analysis provided complementary support for this interpretation. S, Cd, and As were the leading predictors of geochemical affinity between the reference classes. Preindustrial overbank sediments indicate that several metal (loid)s were present before intensive mining, whereas mine tailings and some modern sediments show enhanced mining related geochemical signatures. Environmental risks were spatially heterogeneous, with ecological risk concentrated in specific historical mine tailings deposits (RJ 05 and RJ 07) and human health concerns primarily associated with As in groundwater. These findings provide a screening level assessment for identifying priority areas requiring continued monitoring, particularly groundwater and tailings deposits with elevated ecological risk. Nevertheless, these findings should be considered within the spatial and temporal scope of the present study, and future investigations incorporating broader datasets could further refine the geochemical relationships identified here.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

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

Supplementary material

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

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Footnotes

Edited by: Carla Patinha, University of Aveiro, Portugal

Reviewed by: Kwaku Boakye, Colorado State University, United States

Dickson Abdul-Wahab, University of Ghana, Ghana

Anjali Kumari, Indian Institute of Technology Dhanbad, India