Source apportionment and health risk prioritization of groundwater contamination in a coal-mining area: An integrated SOM-PMF-Monte Carlo framework.

2.1. Study area

2.2. Geology and hydrogeology

The stratigraphy belongs to the Greater Khingan–Yanshan division, specifically the Boketu–Erenhot subdistrict. Exposed strata are relatively simple, mainly including Mesozoic Jurassic formations (Hongqi (J₁h), Tuchengzi (J₂t), Manketouebo (J₃mk), Manitou(J₃mn), Baiyingaolao (J₃b)), Cretaceous Huolin River Formation (K₁d), and the Cenozoic Quaternary System (Q). Jurassic strata are widely distributed and comprise clastic sedimentary rocks (such as feldspathic quartz siltstone), volcanic rocks, and pyroclastic rocks (e.g., andesite, rhyolite, and tuff). The Cretaceous Huolin River Formation (K₁d) consists of coal-bearing clastic rocks, with lithologies including sandstone, mudstone, and conglomerate. The Quaternary and Neogene systems are mainly composed of loose sand-gravel layers and basalt.

The groundwater system is divided into three fourth-level hydrogeological units: Sharhure, Ulgai River, and Dundnuur. Based on occurrence and hydraulic properties, aquifers are classified into five systems: (1) the Quaternary unconsolidated porous and bedrock weathered-fractured porous phreatic aquifer system (Q + J); (2) the Quaternary unconsolidated porous and clastic rock porous-fractured phreatic aquifer system (Q + K); (3) the clastic rock porous-fractured phreatic aquifer system (K); (4) the bedrock fractured phreatic aquifer system (J); and (5) the clastic rock porous-fractured confined aquifer system (K). The two Quaternary unconsolidated porous-fractured aquifer systems are the most widespread and have the best water yield (max single-well yield 1000–3000 m³/d). They are shallow (<30 um) and serve as the main water-extraction targets.

Groundwater recharge mainly originates from atmospheric precipitation infiltration. In the surrounding low to medium mountain areas with exposed bedrock, precipitation rapidly infiltrates through rock fractures to form bedrock fissure water, which discharges into valleys and basins. Natural groundwater flow varies among units: Sharhure flows southwest–northeast; Dundnuur converges south–north (or southeast–northwest). In the unconsolidated and clastic rock aquifers within river valleys and basins, groundwater also receives lateral runoff recharge from the aforementioned mountainous areas and is influenced by seasonal floods and surface water infiltration. However, intensive human activities (groundwater extraction, open-pit mine dewatering) have formed a regional groundwater cone of depression in the Sharhure unit. This has significantly altered the natural flow field, causing groundwater to converge toward extraction and dewatering centers. Artificial extraction and mine dewatering have become the main groundwater discharge pathways (Wang, 2015).

The study area has 14 coal mining districts with a total annual capacity of about 35 million tons. Most are open-pit mines, with coal hosted in the fourth member of the Lower Cretaceous Huolin River Formation. The coal is high-quality aged lignite, low-sulfur (<0.5%) and low-phosphorus (<0.02%), with average calorific value>3100 kcal/kg.

3.1. Sampling and analysis

A total of 38 groundwater samples were collected from the Huolingol mining area between May 12 and 29, 2025, with 25 samples located in Huolingol City and the other 13 situated in the northern part of the adjacent Jarud Banner (Fig. 1c and d). All samples were phreatic water obtained from agricultural irrigation wells, domestic drinking water wells, and coal mine monitoring wells, with well depths ranging from 20 um to 120 um. The water table in unconfined aquifers ranges from about 2 to 25 um (Wang, 2015), ensuring all wells fully penetrate the saturated zone. Sampling points were evenly distributed, ensuring good regional representativeness. Before sampling, each well was purged for approximately 5 min to eliminate the influence of stagnant water in the well. Samples were filtered through 0.45 um membranes and stored in pre-cleaned polyethylene bottles rinsed 2-3 times with local groundwater. Triplicate samples were collected at each site. Sample preservation, transport, and pretreatment strictly followed national specifications. All samples were ultimately sent to Inner Mongolia Xin’antai Testing and Evaluation Technology Co., Ltd. for physicochemical analysis.

Analyzed parameters included pH, major cations (Mg²⁺, Ca²⁺, Na⁺, K⁺), major anions (HCO₃^–, Cl^–, SO₄^2-), permanganate index (COD_Mn), ammonia nitrogen (NH₄⁺-N), nitrate (NO₃^–), total dissolved solids (TDS), total hardness (TH), and fluoride (F^–). pH was measured on-site using a portable pH meter. COD_Mn was determined by the acidic potassium permanganate titration method. NH₄⁺-N was analyzed using ultraviolet-visible spectrophotometry (instrument model UV-9600, serial number 03005). TDS was measured gravimetrically. TH was determined via EDTA-2Na titration. Cations were analyzed by atomic absorption spectrophotometry (instrument model WFX-120B, serial number 03006). HCO₃^– was measured by acid-base indicator titration. The remaining anions (Cl^–, SO₄^2-, NO₃^–, F^–) were analyzed using ion chromatography (instrument model CIC-100, serial number 03117). Sulfate sulfur and oxygen isotopes (?³⁴S-SO₄^2- and ?¹⁸O-SO₄^2-) were analyzed using a stable isotope ratio mass spectrometer (DeltaV Plus) after sulfate precipitation as BaSO₄. Analytical precision was ±0.2‰ for ³⁴S and ±0.3‰ for  ¹⁸O. Ion charge balance errors were below 5% for all samples, indicating reliable data. Heavy metals (As, Cd, Pb, Zn, Fe, Mn, Hg) were mostly below detection limits or far below drinking water standards. They were excluded due to negligible health risks.

3.2. Self-organizing map (SOM)

Groundwater samples were classified using a hybrid approach combining Self-Organizing Map (SOM) and K-means clustering. SOM projects high-dimensional hydrochemical data onto a two-dimensional neuron grid, preserving topological relationships while enabling nonlinear dimensionality reduction (Qu et al., 2022; Appukuttan et al., 2025). The optimal map size was determined based on the quantization error (QE) and topographic error (TE). K-means clustering was then applied to the SOM neuron weights, with the optimal number of clusters determined by the silhouette coefficient and Davies-Bouldin Index (Mao et al., 2021). The analysis was implemented in R.

3.3. Spatial interpolation of hydrochemical parameters

Spatial distribution maps of 14 hydrochemical parameters were generated using inverse distance weighting (IDW, power = 2, four-sector neighborhood) in ArcGIS 10.8 Geostatistical Analyst. Leave-one-out cross-validation (LOOCV) was conducted for each parameter to evaluate performance. Cross-validation results (mean error ME, root mean square error RMSE) are summarized in Table S2.

3.4. Entropy weight-TOPSIS model

Groundwater quality was assessed using the entropy weight-technique for order of preference by similarity to ideal solution (EW-TOPSIS) method, which combines objective entropy weighting with TOPSIS ranking (Chen, 2021). To anchor the evaluation to regulatory benchmarks, virtual standard water samples derived from the Class I-IV thresholds of the Chinese Groundwater Quality Standard (GB/T 14,848-2017) were incorporated into the decision matrix (detailed thresholds are listed in Table S3, and the corresponding virtual sample parameters are provided in Table S4). Entropy weights were calculated as detailed in the Supplementary Material (Table S5). This approach shifts the evaluation from purely relative comparison to standard-based classification.

3.5. Positive matrix factorization (PMF) model

Source apportionment was performed using Positive Matrix Factorization (PMF, EPA PMF 5.0; Norris et al., 2014), a receptor model that decomposes a measured data matrix into factor contributions and profiles under non-negativity constraints (Paatero and Tapper, 1994). The model minimizes an objective function Q based on the measured concentrations and their associated uncertainties, which were calculated from the minimum detection limit and relative standard deviation (detailed equations provided in Supplementary Material). Input data included 12 chemical parameters (TDS, TH, COD_Mn, NH₄⁺-N, NO₃^–, Cl^–, HCO₃^–, F^–, Na⁺, K⁺, Ca²⁺, Mg²⁺) from all 38 sampling sites. Models were run with 4–6 factors and multiple random seeds. The optimal solution was selected based on residual analysis, Q value, and geochemical interpretability.

3.6. Monte Carlo-based human health risk prioritization

Non-carcinogenic health risks from F^–, NO₃^–, and SO₄^2- were assessed following the USEPA exposure model (USEPA, 2004), considering ingestion and dermal contact pathways for adults and children. The inclusion of SO₄^2- as a health risk driver is justified by its potential health effects. High SO₄^2- in drinking water may cause gastrointestinal effects, especially in children (Chien et al., 1968; USEPA, 2004). A reference dose (RfD) of 8.2 mg/kg/day was adopted from Su et al. (2026).

Monte Carlo simulation (10,000 iterations) was performed using Crystal Ball® software to address parameter uncertainties. Input parameters were assigned probabilistic distributions based on published literature (Chen et al., 2023; Su et al., 2026; Yan et al., 2025): ingestion rate (IR) and body weight (BW) followed normal distributions (children: IR 1.20±0.12 L/d, BW 26.8 ± 2.68 kg; adults: IR 2.50 ± 0.25 L/d, BW 70.0 ± 7.00 kg); exposure frequency (365 d/y), duration (6 y for children, 30 y for adults), and averaging time (2190 d for children, 10,950 d for adults) were point estimates. All exposure parameters were assumed independent, which is conservative and widely adopted in similar assessments (Chen et al., 2023; Su et al., 2026). Parameters and RfD are summarized in Table S13. In each iteration, hazard quotients (HQs) were calculated for each contaminant and summed to obtain a probabilistic distribution of the hazard index (HI) (Eqs. S20–S23). Sensitivity analysis identified the most influential parameters.

Deterministic HI values were also calculated using mean exposure parameters and the measured contaminant concentrations at each sampling point. To link PMF sources to health risks, factor contributions to HI were estimated from the PMF contribution matrix and source profiles (Tables S11 and S12), following the method described by (Chen et al., 2023) (Eqs. S24–S26). Spatial distribution of deterministic HI values was generated using IDW interpolation (power=2, four-sector neighborhood; Section 3.3), with LOOCV using adult HI data.

4.1. Hydrochemical classification

1. Download: Download high-res image (944KB)
2. Download: Download full-size image

Fig. 2. Self-organizing maps of 14 variables and the corresponding U-Matrix plot.

Spatially, the three clusters exhibited distinct hydrogeological characteristics (Fig. 1c): Cluster 1 primarily occurred in clastic rock porous-fractured phreatic aquifers, distributed around the depression cone formed by dewatering from open-pit coal mines, representing a typical slow-flow or stagnant runoff zone. Cluster 2 was concentrated in unconsolidated porous and clastic rock porous-fractured phreatic aquifers, mainly located in river confluence zones, corresponding to a typical groundwater stagnation zone. Cluster 3 predominantly occurred in unconsolidated porous and bedrock weathered-fractured porous phreatic aquifers, distributed in the peripheral areas of the sampling network and adjacent to bedrock fractured aquifers, reflecting the hydrogeochemical characteristics of recharge or near-recharge zones.

4.2. Contaminant characteristics and spatial distribution

1. Download: Download high-res image (721KB)
2. Download: Download full-size image

Fig. 4. (a)-(j) Box plots of hydrochemical variables for three groundwater clusters and their permissible limits. Note: The permissible limits refer to the Class III standards of the “Groundwater Quality Standards” (GB/T 14,848-2017).

1. Download: Download high-res image (1MB)
2. Download: Download full-size image

Fig. 5. (a)-(n) Spatial distribution of the hydrochemical variables of groundwater in the study area.

Cluster 1 (central-northeastern slow-flow zones) had significantly higher NH₄⁺-N, F^–, HCO₃^–, Na⁺, and K⁺ than the other clusters (Fig. 4). The elevated Na⁺, K⁺, and HCO₃^– reflect intense silicate weathering and cation exchange under sluggish flow, creating an alkaline, Na-rich environment that promotes F^– desorption and inhibits its precipitation. Consequently, this cluster showed the highest F^– concentrations among all samples, with 44% exceeding the Class III standard (1.0 mg/L), indicating a geogenic fluoride contamination issue exacerbated by natural hydrodynamic conditions. Spatially, high F^– was concentrated in the central-north (Fig. 5h), consistent with Cluster 1 distribution. Widespread NH₄⁺-N exceedances (67% above the Class III standard of 0.50 mg/L) suggest anthropogenic inputs (domestic wastewater, agricultural non-point sources).

Cluster 2 (central-eastern stagnation zones) had the highest concentrations of TDS, TH, COD_Mn, NO₃^–, SO₄^2-, Ca²⁺ and Mg²⁺ (Fig. 4). Mean TDS (1009.25±246.71 mg/L) and TH (611.75±135.32 mg/L) both exceeded the Class III standards, classifying this groundwater as highly mineralized and hard. SO₄^2- emerged as a key contaminant in this cluster, with all samples showing elevated sulfate, probably related to anthropogenic acid inputs from coal mine drainage and industrial wastewater. Correspondingly, the concentrations of TDS, TH, and SO₄^2- in the central-eastern part of the study area were significantly higher than those in other regions (Fig. 5b, c, i). Nitrate contamination was also prevalent: 50% of samples exceeded the Class III standard for NO₃^–, with high-value areas aligning closely with villages characterized by intensive agricultural cultivation (Fig. 5f). Stagnant conditions slow pollutant migration and promote accumulation of NO₃^– and SO₄^2-, making this area a mixed contamination hotspot.

Cluster 3 (eastern, western, southern peripheral recharge zones) had the lowest concentrations across all parameters (Fig. 4). All samples met the Class III standards for all measured indicators, reflecting minimal anthropogenic influence. This area is far from mining, industrial, and agricultural activities, with rapid groundwater renewal that limits contaminant accumulation. Cluster 3 represents the natural background and serves as a reference for contamination assessment.

In summary, groundwater contamination in the study area exhibits clear spatial heterogeneity: geogenic fluoride contamination dominates in Cluster 1 (slow-flow zone), while anthropogenically sourced sulfate and nitrate contamination prevails in Cluster 2 (stagnation zone adjacent to industrial and agricultural activities). Cluster 3 retains natural background quality. These patterns provide a spatial basis for subsequent source apportionment and health risk assessment.

4.3. Groundwater quality evaluation

Based on the EW-TOPSIS model calculation, the five-level classification thresholds for groundwater quality were determined using the relative closeness (Ci) of four virtual groundwater samples (Table S6). The specific classification is as follows: Ci >0.8347 corresponds to “Class I (Excellent)”; 0.7424  Ci<0.8347 to “Class II (Good)”; 0.6028 < Ci<0.7424 to “Class III (Moderate)”; 0.1934 < Ci<0.6028 to “Class IV (Poor)”; and Ci<0.1934 to “Class V (Very Poor)”. The quality categories and Ci values of groundwater samples for each cluster are presented in Fig. S1 and Table S6. The evaluation results indicate that 84.2% of the groundwater samples in the study area were rated as “Moderate” or higher, reflecting generally good overall water quality that was suitable or could be made suitable for drinking after appropriate treatment. The remaining 15.8% of the samples were classified as “Poor” (Class IV), all of which fell within Cluster 1 and Cluster 2. In contrast, all samples in Cluster 3 attained “Moderate” or higher quality categories, demonstrating significantly better water quality compared to Cluster 1 and Cluster 2.

5.1. Hydrochemical processes determining major ions

1. Download: Download high-res image (1MB)
2. Download: Download full-size image

Fig. 6. Hydrochemical processes controlling groundwater chemistry. (a-b) Gibbs diagrams, (c-d) Gaillardet diagrams, and (e-h) ionic ratio correlation diagrams: (e) Cl^– vs. Na⁺, (f) SO₄^2- vs. Ca²⁺, (g) HCO₃^–vs. Ca²⁺, (h) ?¹⁸O-SO₄^2-vs. ³⁴S-SO₄^2-.

To further identify the sources and controlling processes of groundwater chemical components, a series of ionic ratio diagrams were employed (Fig. 6e–h).

The Na⁺/Cl^– ratio diagram (Fig. 6e) shows that most samples have Na⁺/Cl^– ratios well above 1 (1.5–22.7), ruling out halite dissolution as the primary sodium source. Cluster 1 (around mine dewatering cone) had higher ratios, showing strong Na⁺ enrichment. To further distinguish whether Na⁺ enrichment is attributed to silicate weathering or cation exchange, the chloro-alkaline index (CAI-I) was quantitatively evaluated (Fig. S2). The results showed that Cluster 1 had strongly negative CAI-I, confirming intense forward Na⁺-Ca²⁺ cation exchange (Fan et al., 2024). This caused Ca²⁺ and Mg²⁺ depletion and Na⁺ enrichment, creating favorable conditions for F^– enrichment (high pH, Na-rich environment promotes F^– desorption). This was related to the flat hydraulic gradient and sluggish flow around the dewatering cone (Zhang, 2016). Clusters 2 and 3 had weaker cation exchange, with excess Na⁺ mainly from silicate weathering (albite in pyroclastic rocks).

The Ca²⁺/SO₄^2- ratio diagram (Fig. 6f) shows that most samples have Ca²⁺/SO₄^2- > 1 and deviate from the 1:1 line typical of gypsum dissolution, indicating gypsum dissolution was not dominant. PHREEQC saturation indices confirmed all samples undersaturated with gypsum (SI < 0) (Table S7), further ruling out gypsum dissolution as a significant SO₄^2- source. This finding aligned with previous studies reporting that the strata in this region were almost devoid of gypsum (Wang, 2015; Zhang, 2016). Most samples exhibited relative excess of Ca²⁺, which may originated from silicate weathering and/or the dissolution of calcium-silicate minerals driven by external acidic inputs (Zhang et al., 2023). A few samples in Cluster 1 showed elevated SO₄^2- content. Given the limited distribution of sulfide minerals in the study area (Wang, 2015; Zhang, 2016), these SO₄^2- were likely derived from anthropogenic activities, such as coal mine drainage or industrial wastewater discharge (Yang et al., 2024).

The Ca²⁺/HCO₃^– ratio diagram (Fig. 6g) further revealed distinct patterns of water-rock interactions. Most Cluster 3 and some Cluster 1 samples aligned with the 1:1 line. With no carbonate strata (Wang, 2015; Zhang, 2016), this reflects silicate weathering (e.g., anorthite) (Raju et al., 2009). Some Cluster 1 and 3 samples plotted above the line (Ca²⁺/HCO₃^–<1), indicating Ca²⁺ adsorption via cation exchange. All Cluster 2 samples plotted below the line (Ca²⁺/HCO₃^–> 1). Combined with the weak influence of gypsum dissolution (Fig. 6f), this pattern points to external acidic inputs (e.g., anthropogenic sulfuric acid) dissolving calcium-bearing silicate minerals. This process releases substantial Ca²⁺ and Mg²⁺ while introducing SO₄^2- without a corresponding increase in HCO₃^– (Li et al., 2008; Nordstrom et al., 2015), causing high hardness and SO₄^2- in Cluster 2.

Dual isotopes (?³⁴S-SO₄^2- and ?¹⁸O-SO₄^2-) were used to trace sulfate sources in groundwater across the three clusters (Fig. 6h). The isotopic ranges for different sulfate sources (coal mine drainage, atmospheric deposition, fertilizer, sewage, and evaporite dissolution) were adapted from literature (Cui et al., 2025; Su et al., 2026). Two Cluster 2 samples fell in the coal mine drainage field, indicating sulfate from oxidation of sulfur-bearing compounds in coal seams mobilized by mining activities; the other two fell in the sewage field. Cluster 3 samples clustered near atmospheric deposition, with some overlapping the fertilizer and coal mine drainage fields. Cluster 1 samples spanned mine drainage to higher isotopic values, indicating mixing of mining sulfate with agricultural/domestic sources.

Nitrate sources were further distinguished using the Cl^– vs. NO₃^–/Cl^– diagram (Fig. S3). Cluster 1 had low NO₃^–/Cl^– and low Cl^–, showing no significant nitrate pollution. In Cluster 2, a “low Cl^–, high NO₃^–/Cl^–” signature indicated agricultural fertilizer pollution, while a “high Cl^–, moderate NO₃^–/Cl^–” pattern suggested sewage/manure contamination (Su et al., 2026). Cluster 3 had low ratios, reflecting natural background with minor anthropogenic influence.

Integrating the above graphical analyses, the dominant hydrogeochemical processes governing the three groundwater clusters in the study area could be clearly delineated: Cluster 1 was primarily dominated by intense forward cation exchange, with silicate weathering providing the material basis. This process led to severe depletion of Ca²⁺ and Mg²⁺, significant enrichment of Na⁺, and the formation of typical Na-HCO₃ type water; Cluster 2 was mainly driven by anthropogenic pollution, exhibiting characteristics typical of human-altered systems. External acidic inputs promoted strong dissolution of silicate minerals, releasing large amounts of cations such as Ca²⁺, while introducing substantial SO₄^2- and Cl^–, collectively forming a highly mineralized mixed-pollution water body; Cluster 3 represented the regional natural background hydrochemical signature, controlled primarily by silicate weathering of Jurassic volcaniclastic rocks. It exhibited low mineralization and reflected the favorable hydrogeochemical environment of recharge or near-recharge zones. These process insights support source identification and risk assessment of priority pollutants (F^– in Cluster 1; SO₄^2-/NO₃^– in Cluster 2).

5.2. Source identification and apportionment

PMF analysis identified six sources as the optimal solution. The model fit well: residuals mostly within -3 to 3, Q(robust)/Q(true)=0.80. R² for 12 parameters ranged from 0.35 to 0.99 (mean 0.85) (Table S8), showing good agreement between predicted and measured values. Robustness was verified by bootstrap (100 runs, R >0.6) and DISP analyses. All six factors were stably reproduced, with mapping rates 68–100% (Table S9). DISP showed no factor swaps at dQmax=4, 5, 15, indicating no rotational ambiguity (Table S10). These diagnostics support the reliability of the six-factor solution, although Factor 4 should be interpreted cautiously due to a lower mapping rate. Factor profiles and contributions are presented in Table S11 and Fig. 7.

1. Download: Download high-res image (715KB)
2. Download: Download full-size image

Fig. 7. The PMF analytic results (a) Contribution rate of the six factors to each chemical parameters. (b) Proportion of the six factors in source apportionment (Factor 1: cation exchange-silicate weathering; Factor 2: agricultural fertilizers; Factor 3: human-enhanced silicate dissolution; Factor 4: organic matter degradation-carbonate; Factor 5: mixed coal mine drainage and industrial/domestic wastewater; Factor 6: domestic sewage and livestock manure).

Factor 1 made prominent contributions to Na⁺ (33.26%), K⁺ (43.19%), and F^– (49.68%). This pattern is consistent with intense forward cation exchange (Na⁺-Ca²⁺), leading to significant Na ⁺enrichment and corresponding depletion of Ca²⁺ and Mg²⁺. Weathering of fluorine-bearing silicate minerals (e.g., biotite, amphibole) in Jurassic volcaniclastic rocks, enhanced by prolonged water-rock interaction and cation exchange, was the main F^– source. The consistent undersaturation of all samples with respect to fluorite (CaF₂; Table S7) indicated that F^– concentrations were not limited by fluorite solubility, further supporting this interpretation. Accordingly, Factor 1 was identified as a cation exchange-silicate weathering source contributing to both major ion evolution and fluoride enrichment.

Factor 2 showed an exceptionally high contribution to NO₃^– (89.86%), far exceeding other indicators. NO₃^– enrichment in groundwater is typically linked to agricultural fertilizers, domestic sewage and manure discharge, or light industrial wastewater (Xu et al., 2025). However, organic nitrogen sources would also elevate COD_Mn (Hou et al., 2022). In Factor 2, the contribution to COD_Mn was only 7.36%, ruling out significant organic inputs. Thus, the NO₃^– in Factor 2 primarily originated from inorganic agricultural fertilizers. This factor was identified as an agricultural fertilizer source.

Factor 3 showed high contributions to Mg²⁺ (64.77%), Ca²⁺ (60.68%), and SO₄^2- (20.33%). This factor represents anthropogenically enhanced silicate weathering rather than a natural geochemical process. Exogenous H₂SO₄ from coal mine drainage (oxidation of sulfur compounds) accelerates dissolution of Ca-Mg silicates, releasing Ca²⁺, Mg²⁺, SO₄^2- (Li et al., 2008; Nordstrom et al., 2015). This interpretation is supported by SO₄^2- isotopes (Fig. 6h): samples with high Factor 3 contributions (all belonging to Cluster 2, Table S12) fell in the mine drainage field. This confirms that the acidity driving silicate dissolution originates from coal mining activities rather than natural processes. Thus, Factor 3 was therefore identified as a human-enhanced silicate mineral dissolution source.

Factor 4 was dominated by HCO₃^– (73.99%), with moderate contributions to COD_Mn (20.59%), Cl^– (28.99%), and F^– (27.04%), but low contributions to major cations such as Na⁺, Ca^2+, and Mg²⁺. This pattern does not represent typical silicate or carbonate dissolution. Instead, microbial degradation of organic matter in the vadose zone or shallow aquifer generated CO₂, driving carbonate equilibrium toward HCO₃^– formation (Cangemi et al., 2021; Yamanaka, 2012). The concurrent presence of COD_Mn, Cl^–, and F^–further supported inputs from mild domestic sewage or soil organic matter (Zhou et al., 2026). Factor 4 was thus identified as an organic matter degradation-carbonate source.

Factor 5 showed significant contributions to SO₄^2- (77.23%), Cl^– (54.58%), and Na⁺ (45.52%). Sulfate dual isotope data (Fig. 6h) indicated that SO₄^2- mainly from coal mine drainage or wastewater, while Na⁺ and Cl^– co-enrichment also reflected mixed domestic/industrial wastewater inputs (Mu et al., 2024; Yan et al., 2025). High Factor 5 samples were distributed in mining, industrial, agricultural areas (Table S12). Factor 5 was thus identified as a mixed source of coal mine drainage and industrial/domestic wastewater.

Factor 6 showed the highest contribution to NH₄⁺-N (78.15%), followed by K⁺ (29.95%) and Na⁺ (14.99%). These ions are common in domestic sewage and livestock manure: NH₄⁺-N derives from livestock farming and nitrogen-based fertilizers, while K⁺ and Na⁺ come from domestic wastewater and livestock excreta (Jiang et al., 2024; Li et al., 2020). Inadequate sewage and manure management allowed pollutant infiltration. Consequently, Factor 6 was determined as a mixed source of domestic sewage and livestock manure.

The six factor contributions were normalized (Fig. 7b), quantifying each source’s impact on groundwater contamination. The results showed that the most significant influence was attributed to the mixed pollution source of coal mine drainage and industrial/domestic wastewater (Factor 5, 22.42%), followed by silicate mineral dissolution (Factor 3, 19.31%), cation exchange-silicate weathering (Factor 1, 16.40%), organic matter degradation-carbonate (Factor 4, 15.53%), agricultural fertilizers (Factor 2, 13.98%), and domestic sewage and livestock manure (Factor 6, 12.36%). This ranking shows that anthropogenic-driven composite contamination dominated groundwater quality evolution. Mining and industrial activities exerted the strongest pressure. Human-enhanced silicate dissolution ranked second, confirming strong human alteration of natural water-rock interactions. In summary, the six PMF-resolved factors fall into three categories: (i) natural background (Factor 1, 16.40%); (ii) human-enhanced geochemical reactions (Factors 3 and 4, 34.84%); and (iii) direct pollutant inputs (Factors 2, 5, and 6, 48.76%). This classification reveals that direct human inputs and enhanced reactions together accounted for the vast majority of contamination, while natural background played a minor role.

5.3. Human health risk prioritization

5.3.1. Probabilistic risk characterization

Monte Carlo simulation method (10,000 random iterations) was used to assess non-carcinogenic health risks from F^–, SO₄^2-, and NO₃^– via ingestion and dermal contact for adults and children. The results showed that the non-carcinogenic risks associated with these three pollutants differed significantly across population groups (Fig. 8a–d). Regarding the hazard quotient (HQ) for individual pollutants, the probabilities of HQ exceeding the safety threshold (>1) in children were 4.70% for F^–, 3.60% for SO₄^2-, and 1.87% for NO₃^–, whereas the corresponding exceedance probabilities in adults were 1.92%, 2.07% and 1.13%, respectively. These findings indicated that F^– and SO₄^2- were priority pollutants in local groundwater, and children faced higher health risks from all three pollutants compared to adults. The simulated hazard index (HI) showed that children had a 33.93% probability of exceeding the safe threshold (HI>1), 1.75 times higher than adults (19.36%). Therefore, children should be a key protection target in regional groundwater health risk management.

1. Download: Download high-res image (1023KB)
2. Download: Download full-size image

Fig. 8. Human health risk prioritization and source-specific risk contribution. Cumulative distribution functions of HQ for (a) F^–, (b) NO₃^–, (c) SO₄^2- and (d) HI. (e) Box plots showing the distribution of HI values attributed to each PMF-resolved factor. (f) Percentage contribution of each factor to the total HI for Children and Adults. (HQ: hazard quotient; HI: hazard index; Factor 1: cation exchange-silicate weathering; Factor 2: agricultural fertilizers; Factor 3: human-enhanced silicate dissolution; Factor 4: organic matter degradation-carbonate; Factor 5: mixed coal mine drainage and industrial/domestic wastewater; Factor 6: domestic sewage and livestock manure).

5.4. Integrated source-risk assessment and management implications

5.5. Limitations and future research directions

This study has several limitations to be addressed in future work. First, the analysis is based on a single sampling campaign (May 2025). While this represents a stable hydrogeological period, seasonal variations (especially for NO₃^– affected by agriculture and recharge) cannot be captured. Detailed screen and sampling depths were not recorded, limiting vertical water quality characterization. Future studies should include multi-seasonal monitoring, record detailed well construction data, and follow standard protocols requiring parameter stabilization before sample collection to resolve temporal and vertical dynamics and validate source apportionment stability. Furthermore, the health risk assessment has two limitations: the additive HI model does not account for potential synergistic or antagonistic effects among co-occurring contaminants (F^–, NO₃^–, SO₄^2-), and dermal exposure parameters rely on default USEPA values lacking region-specific data. Future studies should explore mixture toxicology models and develop region-specific exposure parameters to refine risk estimates. Third, the PMF receptor model provides a static snapshot of source contributions, whereas the groundwater system is dynamic (preferential flow paths, drawdown cones, seasonal recharge). Future studies should integrate numerical groundwater flow and transport modeling (e.g., MODFLOW) with receptor models to resolve temporal dynamics and simulate source contribution evolution under changing hydrogeological conditions. Additional future work should include: compound-specific isotope analysis (e.g., ?¹⁵N-NO₃) to better constrain nitrate sources; inclusion of organic contaminants (e.g., PAHs) for more comprehensive risk assessment; and development of probabilistic risk-cost optimization frameworks to support cost-effective remediation.