The Dianchi Lake basin (102°29′–103°01′E, 24°29′–25°28′N), located on the Yunnan-Guizhou Plateau in southwestern China, is a densely populated and highly urbanized region with intensive socioeconomic activities (Figure 1). Dewatered sewage sludge samples were collected from 27 municipal WWTPs, (designated D1–D27) surrounding Dianchi Lake in Kunming City (Figure 1; Supplementary Table S1). The spatial distribution of the sampling sites was designed to reflect the regional population density gradient, which decreases clockwise from the north to the west, ensuring representative coverage. Specifically, WWTPs D1–D11 are situated in the densely populated northwestern area (encompassing Xishan, Wuhua, and Panlong Districts). D12–D16 are located in the moderately populated eastern sector (including areas near Changshui Airport and Guandu District). D17–D20 serve the less urbanized Chenggong District, while D21–D25 are in the southern Jinning District. Finally, D26 and D27 are found in the sparsely populated Haikou Subdistrict. Sampling was performed quarterly from May 2022 to February 2025, in accordance with the Chinese national standard method (GB 24188–2009).

The total concentrations of eight HMs mercury (Hg), arsenic (As), cadmium (Cd), lead (Pb), copper (Cu), nickel (Ni), chromium (Cr), and zinc (Zn) in the sewage sludge were determined following Chinese National Standard methods. Samples were digested using a microwave-assisted system with an HCl–HNO₃–HF–H₂O₂ mixture. After evaporation to near dryness at 130 °C, the residues were diluted to volume with 2% nitric acid, filtered, and analyzed by Inductively Coupled Plasma Mass Spectrometry (ICP-MS). MOS-grade reagents and ultrapure water were used throughout the process. Quality control included procedural blanks, parallel samples (≥10%), and spiked recovery samples (≥10%). Certified reference materials (Supplementary Table S2) were analyzed to ensure data accuracy, and all procedures complied with the “Environmental Water Quality Monitoring Quality Assurance Handbook” (He et al., 2016; Zhao et al., 2008).

The Nemerow Pollution Index model (Chen et al., 2024) integrates both the maximum and average values of single-factor pollution indices, thereby emphasizing the impact of the most contaminated element. The calculation formulas are shown in Supplementary Table S2. The pollution degree is classified into five levels based on the P _N value, as shown in Supplementary Table S3.

To more sensitively identify high-risk WWTPs, this study employed the Weighted Comprehensive Risk Score (WCRS) method (Jiang et al., 2025; Xia et al., 2024). Based on the single-factor pollution index (P _i , Supplementary Table S2), this method calculates a comprehensive risk score according to the degree of exceedance of HMs limits and their proximity to threshold values (the risk grading criteria and weight design are shown in Supplementary Table S4). It should be noted that the weights of WCRS are derived from empirical analysis, with core calculation basis including exceeding frequency, exceeding degree, and proximity to compliance limits. Toxicological weights, bioavailability, or uncertainty analysis are not included, making it a management-oriented tool for risk priority ranking and early warning. The assessment in this study was primarily based on the “Sewage Sludge Quality from Municipal Wastewater Treatment Plant” standard (GB24188-2009) (Supplementary Table S2).

The Potential Ecological Risk Index method, proposed by the Swedish scholar Hakanson (Hakanson, 1980) is based on sedimentological theory. It comprehensively considers the toxicity, migration potential, and environmental behavior characteristics of HMs to assess their pollution degree and potential ecological risk in sediments. The calculation formulas are shown in Supplementary Table S2. The prescribed toxic response factors (T _r ) for Cd, Zn, Cr, Cu, Ni, As, Hg, and Pb are 30, 1, 2, 5, 5, 10, 40, and 5, respectively. The reference values were set according to the Grade II standard (pH 6.5–7.5) of the “Environmental Quality Standard for Soils” (GB 15618–1995). The classification criteria for the ecological risk factor (E _r ) and index (RI) are provided in Supplementary Table S5.

To identify the sources and interrelationships of the HMs, the following multivariate statistical methods were comprehensively applied. Pearson correlation analysis (Xia et al., 2024; Wang et al., 2025) was performed based on the single-factor pollution indices (P _i = C _i /S _i ) to calculate correlation coefficients between the risk indices of different metals and their significance (p-values). Data normality was checked using the Shapiro-Wilk test. Local correlation analysis was conducted specifically for exceeding metals in high-risk WWTPs, with results corrected using the Bonferroni method.

Hierarchical cluster analysis (Q-mode clustering) was applied to classify the eight HMs based on their annual average concentrations. Data were standardized using Z-scores. The Euclidean distance was used as the measure, and the average linkage method was employed to construct the cluster structure, visualized via a dendrogram (Bux et al., 2023).

Principal Component Analysis (PCA) (Bux et al., 2023) is a statistical method for data dimensionality reduction and information extraction. PCA was performed on the monitoring dataset (27 WWTPs × 12 quarters × 8 HMs). Data were standardized using Z-scores. The suitability of the data for PCA was confirmed by the Kaiser-Meyer-Olkin (KMO) measure (0.72) and Bartlett’s test of sphericity (p < 0.001). For comparative validation, non-linear dimensionality reduction techniques, namely t-Distributed Stochastic Neighbor Embedding (t-SNE, perplexity = 30) and Uniform Manifold Approximation and Projection (UMAP), were also employed.

All data analysis, result processing, and visualization in this study were conducted on the Jupyter Notebook platform. The NumPy (v1.26.0) and Pandas (v2.1.1) libraries of Python were employed for raw data cleaning, standardization, and extraction of statistical characteristics. Principal Component Analysis (PCA), Pearson correlation analysis, and Shapiro-Wilk normality test were implemented using the scikit-learn (v1.3.0) and SciPy (v1.11.2) libraries. Q-mode hierarchical cluster analysis was accomplished with the hierarchical clustering module of the scikit-learn library, and the dendrogram visualization was realized by relying on the seaborn (v0.12.2) and matplotlib (v3.8.0) libraries. Nonlinear dimensionality reduction analyses (t-SNE and UMAP) were performed by invoking the TSNE module of the scikit-learn library and the umap-learn (v0.5.3) package, respectively. Spatiotemporal distribution maps of HMs concentrations, risk index classification maps, and statistical charts were collaboratively plotted using the matplotlib and seaborn libraries. The spatial distribution map of sampling sites and related spatial risk distribution maps were compiled with ArcGIS 10.2 software, and the base maps were uniformly projected in the Plate Carrée projection (WGS84 coordinate system). All analytical processes used a fixed random seed to ensure reproducibility.

Three Chinese national standards were selected as evaluation criteria in this study: “Sludge Disposal from Municipal Wastewater Treatment Plants - Quality of Sludge for Landscaping” (GB/T 23488–2009, special for acidic soils), “Quality of Sludge from Municipal Wastewater Treatment Plants” (GB/T 24188–2009), and “Control Standards for Pollutants in Sludge for Agricultural Use” (GB 4284–2018). The selection basis is as follows: these three standards not only comprehensively cover key application scenarios such as general benchmarks for sludge disposal, special standards for landscaping in acidic soils, and special standards for agricultural use, but also have clear legal authority and wide practical universality, which can accurately align with the core theme of this study focusing on sludge resource utilization quality control and environmental risk prevention and control.

Compared with international representative standards (US EPA “Standards for Land Application and Surface Disposal of Sludge”, EU “Directive on the Limitation of Pollutants in Sludge for Agricultural Use” (86/278/EEC) and its revised version, supporting standards of Japan’s “Sludge Reuse Promotion Act”), their core commonalities are: all take municipal sewage sludge as the analysis carrier, core control indicators focus on key dimensions such as heavy metals, pathogens, and physical and chemical properties, limit setting is based on environmental risk thresholds as the core logic, and the scope of application collectively covers global mainstream sludge resource utilization paths such as landscaping and agricultural use. Therefore, the domestic standards selected in this study have high comparability with international standards. In addition, the current mainstream municipal wastewater treatment processes (e.g., activated sludge process) are highly similar globally, resulting in small differences in sludge generation characteristics, and the research results can feed back the standard-setting work of countries in similar regions.

To ensure the reliability of risk models and dimensionality reduction analysis results, the following validation and uncertainty analysis were conducted:

5-fold stratified cross-validation was used to ensure the reliability of t-SNE and UMAP models. A total of 335 sludge heavy metal samples were stratified according to spatiotemporal characteristics to ensure the consistency of the distribution of training/test sets; key parameters were optimized through grid search (t-SNE: perplexity 5–50, learning rate 10–1,000; UMAP: neighborhood size 5–50, min distance 0.01–0.5), the optimal parameters were determined with the goal of maximizing the silhouette coefficient, and the generalization ability of the model was verified through the test set.

Uncertainty sources and propagation: ① Data-driven uncertainty: Derived from the spatiotemporal heterogeneity of samples and high-dimensional data noise, which affects the clustering boundary through sample input propagation; ② Parameter uncertainty: Differences in key parameter values change the distance weight, which is propagated through parameter optimization and has a more significant impact under low sample size; ③ Inherent algorithm uncertainty: Differences in the focus of t-SNE and UMAP on data structure lead to result differences, which affect the interpretation of pollution sources through core logic propagation. The above uncertainties propagate step by step along “sample input—parameter optimization—dimensionality reduction output”, which may affect the accuracy of HMs spatiotemporal characteristic identification and source apportionment.

Uncertainty control measures: Through strict sample quality control (e.g., standardizing sampling and pretreatment processes), multiple rounds of parameter iterative optimization, and cross-validation of multi-algorithm (PCA/t-SNE/UMAP) results, the impact of uncertainty on the accuracy of research conclusions is minimized as much as possible.

Statistical summaries of HMs concentrations in sewage sludge samples are provided in Supplementary Table S6 (see Supplementary Figures S1–S8). The median concentrations of metals followed the order: Zn > Cu > Cr > Pb > Ni > As > Hg > Cd. The Kolmogorov-Smirnov test showed that all elements except Cr significantly deviated from the normal distribution. Cd, Zn, Cu, As, and Pb exhibited strong positive skewness (skewness coefficients: 3.68–11.26); the extreme skewness of As (11.26) and Cd (5.33) was likely caused by intermittent industrial point-source discharges (e.g., electronics manufacturing and electroplating) (Liu and Wang, 2024; Durairaj, 2024; Shao et al., 2024). One-way ANOVA confirmed significant seasonal and spatial variability in HM concentrations (F = 12.36, p < 0.001 for seasonal variation; F = 8.72, p < 0.001 for spatial variation).

The single-factor pollution index (P _i ) assessment over 12 quarters identified As as the most widespread exceeding pollutant, with 20 exceedance events across 8 WWTPs (Figure 2a). WWTP D1 had the most severe pollution, with five consecutive quarters of As exceedance (Q2 2023–Q3 2024). Seasonal analysis showed As exceedance in six WWTPs (Figure 2b): D2, D15, and D21 had higher risks in the wet season, while D3, D12, and D13 showed elevated risks in the dry season.

Cd exhibited distinct seasonal and spatial characteristics, with exceedances occurring primarily in Q1 (dry season) (Figure 2c). Cd pollution was concentrated in D3, D12, and D15 (Figure 2d), with all exceedances recorded in the dry season—a pattern attributed to rainfall dilution in the wet season versus concentrated industrial wastewater discharge in the dry season. Zn pollution was minor, with isolated exceedances in D5 and D6 (Figure 2e,f); D5 had higher Zn concentrations in the dry season, while D6 showed the opposite trend, indicating potential non-point source pollution from rainfall runoff and galvanized pipe corrosion (Li et al., 2019; Zhang P. et al., 2024; Zhang L. et al., 2024).

Overall, As, Cd, and Zn risks were lower in the wet season than in the dry season. Spatially, the combined concentration of these three metals decreased clockwise along the lakeshore from D2 to D26. Hg concentrations were significantly higher in the wet season (p < 0.05), suggesting rainfall-facilitated migration; Pb concentrations increased from north to south, peaking at D24 (potentially linked to high population density and limited WWTP capacity in Jinning District). Cr and Ni showed no clear spatiotemporal patterns. Key spatiotemporal risk characteristics are summarized in Supplementary Table S7 (Ranking of Heavy Metal Risk Levels by Season and WWTP Location).

The Nemerow comprehensive pollution index P _N for each HMs was calculated with reference to three Chinese national standards (Olejnik, 2024) (Supplementary Table S1), with the results detailed in Supplementary Tables S8-S10. The assessed pollution characteristics varied considerably depending on the standard applied. Exceedances of Cd and As were most pronounced under the GB/T 23486-2009 standard. In contrast, only As at Plant D1 exceeded the limit under the GB 24188-2009 standard, largely due to its more lenient thresholds for metals such as Cd. Under the GB 4284-2018 standard, As at Plant D1 remained in exceedance. Notably, Cr, Cu, Ni, Hg, and Pb did not exceed the limits in any of the three standards, indicating a relatively low pollution risk from these elements. Plants D1, D2, and D3 exhibited significantly higher exceedance frequencies than others.

D1 had P _N values > 1.0 under all three standards (2.24, 1.62, and 1.66, respectively; Figure 3), indicating a consistent pollution risk. The GB/T 23486-2009 standard was the most stringent, identifying three plants (D1, D3, D12) as medium-to-high risk. Both the GB 24188–2009 and GB 4284-2018 standards classified only Plant D1 as lightly polluted (Level 3), despite the latter imposing stricter limits for Zn and Cr. The P _N values for Plants D2 and D3 exceeded 0.7 under all standards, placing them in a risk warning state. Their geographical proximity to Plant D1 suggests potential common pollution sources and cumulative risks, warranting priority investigation for source apportionment. Statistical summary of risk levels (Figure 3, inset) shows that the proportion of plants classified as safe (Level 1) increased from 74% under GB/T 23486–2009 to 89% under both the GB 24188–2009 and GB 4284-2018 standards, demonstrating that the stringency of the standard significantly influences the assessment outcome.

The contribution rates of individual HMs to the comprehensive P _N value also varied with the standard employed (Supplementary Figure S3). Under GB 24188–2009, Cd (0.7%) and As (0.6%) were the primary contributors, followed by Hg (0.3%). Under the stricter GB/T 23486-2009 standard, the contribution rate of Cd increased substantially to 11.0%, that of Hg rose to 8.2%, while As remained stable at 0.6%. Under GB 4284–2018, Cd (1.2%) and Hg (0.9%) remained the highest contributors. The contribution rate of As was stable (approximately 0.6%) across all standards, consistently ranking among the top three, whereas Zn and Cu consistently demonstrated the lowest contributions (<0.02%). These results robustly identify Cd, Hg, and As as the core pollutants in the study area.

These findings underscore that the Nemerow index assessment is highly dependent on the stringency of the standard limits. Under lenient standards (GB 24188–2009), the index tends to underestimate pollution levels, whereas under stringent standards (GB/T 23486–2009), it becomes more sensitive and effective at identifying potential risk areas. The Nemerow index is suitable for compliance screening based on sludge disposal pathways (e.g., agricultural use, landscaping) but has limitations: it reflects overall pollution levels mathematically but does not account for HMs biotoxicity or ecological hazards.

To enhance the sensitivity of identifying high-risk WWTPs, this study applied the Weighted Comprehensive Risk Score (WCRS) method, with all assessments conducted against the GB24188-2009 standard. The spatial distribution of pollutant exceedances revealed distinct patterns for key HMs (Figures 4a–d). As was identified as the most severe contaminant, with an exceedance rate of 6.0% (20 out of 335 samples). These exceedances were geographically concentrated at Plants D1 (5 exceedances), D2 (6 exceedances), D3 (4 exceedances), and D15 (2 exceedances), suggesting the potential for persistent industrial discharges or other significant As sources in their catchment areas. In contrast, Cd demonstrated a much lower overall exceedance rate (0.9%, 3/335). Its exceedances were isolated, occurring solely at Plants D3, D12, and D15, which is indicative of localized point source pollution that requires targeted attention. Zn exhibited a similarly low exceedance rate of 0.6% (2/335), with incidents located at Plants D5 and D6. These occurrences may be linked to discharges from electroplating or metallurgical activities. Notably, Hg did not exceed the standard limit in any sample, implying relatively effective control of this element in the incoming wastewater to the WWTPs. Furthermore, as supported by Supplementary Figures S16-S19, no exceedances were detected for Cr, Pb, Cu, or Ni. The median risk coefficients for these metals were substantially below 1, reflecting well-controlled input concentrations throughout the monitoring period. Synthesizing these findings (Figure 4d; Supplementary Table S11), the pollutants were categorized hierarchically: As was unequivocally the primary pollutant, with Cd and Zn constituting secondary pollutants of concern.

The application of the WCRS method, utilizing the weighting criteria in Supplementary Table S4, enabled a refined risk classification for plants with recorded exceedances (detailed results in Supplementary Table S12). The risk categorization was as follows. Medium Risk: Plant D3, with a score of 5 points due to concurrent exceedances of As and Cd; Plants D12 and D15, each scoring 4 points for exceeding two HMs; and Plant D1, scoring 3 points primarily for As exceedance coupled with the proximity of other metals to their thresholds. Low Risk: Plants D5 and D6, scoring 2 points for Zn exceedance only; and other plants, such as D2 and D21, which also scored 2 points for a single As exceedance. The key advantage of the WCRS method lies in its ability to identify “critical pollution states.”

Unlike the traditional Nemerow index, which focuses primarily on pollutants that have already exceeded standards, the WCRS method incorporates points for metals approaching their thresholds. This provides an early warning for potentially escalating risks, shifting the management focus towards prevention. Kappa consistency test confirmed good agreement between WCRS and Hakanson risk grades (Kappa = 0.76, p < 0.001). It should be noted that, as a management-oriented tool, the WCRS has inherent limitations: it does not consider toxicological differences between heavy metals, bioavailability of heavy metals in sludge, or uncertainties in monitoring data. These factors may restrict its ability to reflect the actual ecological toxicity risks of heavy metals, and thus it should be used in conjunction with ecotoxicity-oriented indices such as the Hakanson index rather than as a standalone risk assessment method.

Sewage sludge that meets the criteria specified in the standard Disposal of Sludge from Municipal Wastewater Treatment Plants - Quality of Sludge for Agricultural Use (GB 24188–2009) can be utilized for land-based resource recovery. In this study, the Hakanson potential ecological risk index method was adopted, based on the similarity of toxicological interaction mechanisms between HMs and biological receptors in sludge and those in sediments (Lei, 2023; Yu et al., 2023; Geng et al., 2021; Buthnoo et al., 2025). The model identified Cd and Hg as the predominant contributors to the overall ecological risk (Figure 5a). Cd posed a substantial risk across the study area. The ecological risk factor (Er) for Cd far exceeded the threshold for “very high risk” (Er > 320) at several plants, including D1 (1,065.5), D3 (819.5), D2 (486.8), D12 (551.8), D15 (461.8), and D8 (395.2). Eleven additional plants were classified as “high risk” (Er = 170.1–291.3), and ten plants posed a “considerable risk” (Er = 100.3–137.0).

In contrast, Hg exhibited a clear spatial gradient in risk. Most plants from D1 to D14 presented a “high risk” (Er = 192.2–268.8), whereas those from D15 to D27 were predominantly at “considerable risk” (Er = 53.9–155.3), with the exception of D16, which showed a “moderate risk” (53.9). Apart from a “moderate risk” (Er = 64.3) at Plant D1, As was of low risk elsewhere. Zn, Cr, Cu, Ni, and Pb consistently posed a low risk (Er ≤ 40) across all plants.

The spatial relationship between the single-factor pollution index and the potential ecological risk index was visualized in a bubble plot (Figure 5b), revealing a significant positive correlation between Cd and the comprehensive Risk Index (RI) (r = 0.82, p < 0.001). Plants D1, D3, D2, and D12 formed a distinct cluster of extreme risk, classified at the “very high” pollution level, with Plant D1 exhibiting the peak value. Plant D15 was also associated with a strong ecological hazard. Densely populated areas (encompassing Plants D4–D11, D13–D14, D20–D21) constituted a zone of concentrated “strong risk,” where Cd levels and the RI were strongly coupled.

A clear spatial pattern emerged, with the core area around Plants D1–D3 exhibiting the highest risk, which decreased radially towards the periphery. Plants D16–D19 and D22–D25 were characterized by “moderate risk,” while a localized increase observed at Plants D26–D27 suggests potential emerging pollution sources. This risk gradient was visually pronounced, as the bubble diameters in high-risk areas (yellow) were approximately 1.8 times larger than those in moderate-risk areas (green), underscoring a strong spatial consistency between pollution intensity and risk level.

The fundamental characteristic of the Hakanson index is its transformation of pollution concentrations into ecological risks through toxic response factors. Meanwhile, it should be noted that the toxic response factors in the Hakanson index are calibrated based on the properties of soils and sediments, and no specific adjustments have been made for the unique characteristics of sludge (such as high organic matter content and high water content). Therefore, the risk estimation results may be conservative. This study positioned the Hakanson index as a tool for relative risk comparison and spatial prioritization. The assessment results showed that cadmium (Cd) and mercury (Hg) were the dominant ecological risk drivers, Cd and mercury (Hg) accounted for 55.4% and 34.6% of the total ecological risk, followed by arsenic (As); spatially, the sludge from WWTPs D1, D3, D12, and D15 exhibited relatively high ecological risks, which was consistent with the spatial distribution characteristics of HMs concentrations. This further confirms that the application of the Hakanson index on the basis of clear premise assumptions can effectively identify high-risk areas and key risk factors in sludge, providing a reliable basis for the spatial prioritization of sludge risk management.

Comparison of risk assessment results from three indices summarizes the complementary advantages and mutually verified assessment of the Nemerow index, WCRS, and Hakanson index, constructing a comprehensive risk assessment framework from regulatory, managerial, and ecological perspectives: the Nemerow index and WCRS identify pollutants of concern (As, Cd) and risk areas (northeast), while the Hakanson index further reveals which of these pollutants (Cd, Hg) pose the greatest potential threat to the ecosystem.

To systematically identify the sources and interrelationships of HMs in sewage sludge from WWTPs in the Dianchi Lake basin, this study integrated multiple statistical methods, including correlation analysis, cluster analysis, and principal component analysis (PCA). Correlation analysis revealed significant co-variation among the elements (Figure 6a). Four element pairs—Cd-Cu, Cd-Ni, Cd-As, and Cu-Ni—exhibited highly significant positive correlations (r = 0.43–0.46, p < 0.001). These specific combinations serve as strong evidence for composite pollution sources. The Cd-Cu association is potentially linked not only to electroplating wastewater but also to agricultural activities. According to Wang (2014), pesticides (Wang et al., 2018), fertilizers, and organic fertilizers (Yang et al., 2021) used in the intensive agricultural greenhouses along the eastern coast of Dianchi Lake are significant sources of Cd and Cu, which enter the sewage system via surface runoff. The Cu-Ni combination is a typical marker for industrial wastewater discharge from electroplating (Peng et al., 2023), the Cd-Ni combination is commonly associated with the production and disposal of batteries, and the Cd-As and Cu-Ni combinations are closely linked to non-ferrous metal smelting activities (Luo et al., 2023; Peng et al., 2023). Historically concentrated paper mills and smelters in northern Dianchi are potential sources of such pollution (Yu et al., 2013). In contrast, the correlation between Cr and Ni, Pb may suggest another industrial source, such as stainless-steel production or leather tanning, while the weak correlation between Cr and Cd/As indicates independent input pathways for these elements. These inferences regarding industrial origins are highly consistent with the findings of Wang (2014) on the sector-specific characteristics of HMs in the Dianchi Lake basin and are strongly corroborated by the spatial distribution pattern: the WWTPs with the highest concentrations of these elements (D1, D3, D12) are all concentrated in northern Dianchi, an area identified in official planning documents (Li et al., 2018) as a traditional agglomeration zone for Kunming’s machinery manufacturing, metal processing, and electronics industries.

Cluster analysis further classified the eight HMs into three groups with distinct sources (Figure 6b). Cluster 1 (Cu, Cd, As) exhibited the smallest within-cluster distance, indicating a highly homogeneous source, which can be strongly attributed to “Metal Processing and Smelting Industrial Sources”. Cluster 2 (Ni, Cr, Pb) likely represents “Mixed Industrial Sources”. Previous studies (Luo et al., 2023) have found that Pb primarily originates from vehicle emissions and industries involved in smelting, manufacturing, and using Pb products. Wang (2014), combining spatial distribution characteristics, also inferred that Pb likely stems mainly from industrial production activities. Cluster 3 (Zn, Hg) was geographically distinct, with its high-value zone (D5-D7) located in the old city center of Kunming. The association of Zn and Hg is widely recognized as an indicator of “Combustion Sources” such as coal burning and municipal waste incineration (Zeng et al., 2021), consistent with the area’s history and current status as a densely populated center of energy consumption.

Principal Component Analysis (PCA) quantitatively resolved the contribution of each source (Figure 6c). PC1 (41.8% variance explained) showed high loadings on Cd, Cu, and As, clearly corresponding to the industrial source represented by Cluster 1. PC2 (22.4% variance explained) was dominated by Zn and Hg, confirming the combustion source represented by Cluster 3. The separation of pollution grades identified by PCA was consistently validated by non-linear dimensionality reduction techniques, namely t-SNE (Figure 6d) and UMAP (Figure 6e), demonstrating the robustness of the source apportionment results. The tight clustering of Grade 1 (clean) WWTPs suggests their catchment areas are dominated by domestic wastewater with single and low-intensity HMs inputs. In contrast, the dispersion of Grade 3 and 4 (moderately and heavily polluted) WWTPs reflects their exposure to multiple, high-intensity industrial or mixed pollution inputs.

Spatial analysis provided a geographical anchor for the above source apportionment (Figure 7). The study area exhibited a significant “high in the northeast, low in the southwest” pollution pattern. This macroscopic gradient is highly coupled with the industrial layout and urbanization development level of Kunming City. The high-risk cluster in the northeast (D1-D3, D12) is directly downstream of the main industrial zones; their extremely high Cd and As pollution indices are a direct manifestation of point source discharges from non-ferrous metal-related industries. Furthermore, the observed decreasing pollution trend along the Baoxiang River (D12 to D13, D14) and the Panlong River-Mingtong River systems (D4 to D11) confirms the location of point sources (near D12) and also reveals the self-purification capacity of the river systems or dilution effects from tributaries. Overall, except for the northeastern hotspot, WWTPs closer to Lake Dianchi itself generally showed lower sewage sludge pollution indices, proving that the existing sewage treatment network acts as an effective barrier, playing a crucial role in intercepting HMs pollutants within the basin and protecting the water quality of Lake Dianchi.

In summary, by coupling multivariate statistics, spatial analysis, and regional anthropogenic activity data, this study advanced the understanding of HMs sources in sewage sludge from statistical inference to geographical substantiation. The results indicate that HMs in Kunming’s WWTPs sewage sludge primarily originate from two core processes: (1) Metal processing and smelting industrial activities in the northern and northeastern areas, dominating the input of Cd, Cu, and As; (2) Historical and current energy consumption (coal burning, waste incineration) in the central urban area, contributing most of the Zn and Hg. This refined source apportionment provides a solid scientific basis for implementing a targeted source control strategy characterized by “Zoning, Classification, and Grading”.

To elucidate the relationship between different wastewater treatment processes and the characteristics of HMs accumulation in sewage sludge, this study evaluated 13 technologies across 27 WWTPs. It is critical to emphasize that the HMs concentration in sewage sludge primarily reflects the transfer and enrichment capacity of the process from the aqueous to the solid phase, rather than directly indicating purification efficiency of the effluent. Comprehensive analysis (Figure 8; Supplementary Figure S20) revealed significant differences in the overall HMs accumulation levels in sewage sludge produced by different processes (F = 9.42, p < 0.001). The Sequencing Batch Reactor (SBR) yielded sewage sludge with the lowest comprehensive HMs accumulation (comprehensive P _i = 0.084). Conversely, the Membrane Bioreactor (MBR) and Intermittent Cycle Extended Aeration System (ICEAS) produced sewage sludge with the highest accumulation levels (comprehensive P _i = 0.280 and 0.250, respectively). The highly efficient solid-liquid separation capability of MBR is a key reason for its significant retention and enrichment of metals like Zn. Modified A²/O processes (e.g., A²/O-M, A²/O-MBR) also performed better than the conventional A²/O process.

Regarding individual metals (Figure 8a-h), As accumulation was highest in sewage sludge from ICEAS and MBR processes, while modified processes reduced As enrichment by 48%–74% compared to the conventional A²/O. Cd enrichment was significant in sewage sludge from A²/O-PC and MBR processes but remained at very low levels in various modified processes. The prominent enrichment of Zn in MBR sewage sludge further confirms its highly effective retention capability.

Discussion and Implications: This study establishes that process type is the predominant factor determining the HMs accumulation characteristics in sewage sludge. This provides direct guidance for process selection aimed at sewage sludge safe disposal and resource recovery: for scenarios targeting sewage sludge resource utilization, processes with low HM accumulation (e.g., SBR) are optional, while for scenarios prioritizing effluent water quality control, processes with high HM retention efficiency (e.g., MBR, ICEAS) are more suitable. for new or upgrading WWTPs where sewage sludge recycling is a priority, SBR or modified A²/O processes should be considered; for existing plants employing high-accumulation processes like MBR, the strategic focus must be on enhancing source control and pre-treatment to manage the incoming HMs load at its origin. However, a critical trade-off is highlighted: processes like MBR, while ensuring excellent effluent quality, lead to high HMs accumulation in sewage sludge, thereby increasing subsequent disposal challenges and potential environmental risks. Although this field-based study reveals strong correlations, future controlled pilot-scale studies are needed to isolate the influence of influent loadings and precisely quantify the inter-phase transfer efficiency of HMs across different processes.

Our finding that SBR achieves the lowest sludge HM accumulation aligns with EU research demonstrating that optimized treatment processes reduce sludge Hg and Cd concentrations under strict regulations (European Environment Agency, 2025; Kirchmann et al., 2017). However, EU studies focus on plain river basins, while our work reveals the unique dry-season high-risk and “northeast-southwest” spatial gradient of Cd pollution in the Dianchi Lake plateau lacustrine watershed—filling the research gap in plateau lake sludge HM risk. A study in Haryana, India also identified industrial activities as the main driver of sludge HM spatial differences (Maheshwari et al., 2008), consistent with our source apportionment results. Notably, As was not detected in the Indian study, whereas our study recorded 100% As detection rate—attributed to intensive non-ferrous metal mining in the Dianchi Lake basin. Globally, most studies rely on single concentration analysis for risk assessment (Dai et al., 2007; Scancar et al., 2001); our use of the Hakanson index to quantify Cd/Hg risk contributions (55.4% and 34.6%) further validates the necessity of comprehensive ecological risk assessment. These comparisons highlight the global universality of industrial pollution impacts and the regional specificity of the Dianchi Lake basin, enhancing the study’s international academic value.