Contamination Characteristics and Source Apportionment of Heavy Metal in the Topsoil of a Small Watershed in South Taihang

Abstract

With the rapid development of industrialization and urbanization, the issue of soil environmental pollution is becoming more and more prominent, especially concerning heavy metal contamination, which has garnered significant scholarly attention. The surface watershed formed by waterline is influenced by various factors such as topography, industrial emissions, and agricultural runoff, resulting in a complex process of migration and accumulation of heavy metal elements from multiple sources. In this study, the pollution characteristics and sources of heavy metal elements Hg, As, Pb, Ni, Cd, Cr, Cu and Zn in 165 surface soil samples from the Manghe River watershed in Jiyuan City were comprehensively analyzed using a variety of methods, including statistics, geostatistics, enriched factor analysis and the Positive Matrix Factorization Model (PMF). The results showed that the concentrations of Hg, Cd, As, Cu, Pb and Zn exceeded their corresponding background values with varying degrees of enrichment. Notably, the average contents of Cd, Hg and Pb were 26.70 times, 3.69 times and 4.49 times higher than those in Chinese soils on average, respectively, showing obvious enrichment characteristics. Moreover, there were distinct spatial distribution patterns for each heavy metal element; Ni and Cr exhibited similar trends mainly controlled by the parent material, while human activities significantly affect the other six elements forming high-value areas around mining and related industries. It is noteworthy that Cu, Hg and Zn were influenced by dominant wind direction in autumn and winter, forming sub-high-value zones in southern forested areas; meanwhile, Cu and Zn were also influenced by agricultural fertilizer application as well as surface runoff, leading to secondary high-value areas in the dryland areas. Further analysis revealed a significant positive correlation among these heavy metal elements, suggesting that they may share common sources. Through the PMF Model, four main factors were identified, with factor 2 (36.25%), factor 1 (23.00%), factor 3 (21.20%) and factor 4 (19.55%) ranked in descending order of contribution rate. The heavy metal pollution in the study area was attributed to anthropogenic activities and natural factors, accounting for 63.75% and 36.25%, respectively. Coal mining, chemical industry smelting, vehicle emissions and excessive use of agrochemicals were identified as the main sources of heavy metal pollution. These pollutants entered the soil through direct emissions, atmospheric deposition, transportation and agricultural activities, exerting a significant impact on the soil environment. Therefore, delving into the spatial distribution pattern of soil heavy metal pollution and precise analysis of its sources are of great importance for effective treatment and remediation of soil heavy metal pollution in small watersheds, maintaining healthy soil ecology and safeguarding human health.

1. Introduction

Soil serves as the fundamental basis for agricultural production and plays a vital role in upholding the ecological stability of our planet [1]. However, the process of industrialization and urbanization has resulted in severe contamination of the soil environment, particularly with heavy metals [2]. Globally, heavy metal pollution affects approximately 5 million hectares of land and causes annual economic losses in excess of $10 billion [3,4]. Soil heavy metals are insidious, highly toxic and persistent substances [5] that, when surpassing the threshold levels, contaminate both soil and water resources [6]. This contamination detrimentally affects crop growth and compromises the safety of agricultural products [7], thereby posing a threat to human health [8]. Therefore, comprehending the spatial variability and sources of heavy metal pollution is crucial for effective soil conservation and remediation efforts [9,10].

Heavy metals are spatially correlated and their spatial variability is influenced by both natural environmental factors and human activities at different locations and intensities [11]. Geostatistical modeling has been extensively employed in pollution studies to accurately simulate its spatial structure and variability by fitting variogram functions and using kriging interpolation [12,13]. Source identification and attribution plays a pivotal role in the management of heavy metal pollution in soil. Various multivariate statistical analyses and mathematical models, such as the Chemical Mass Balance method, Positive Matrix Factorization (PMF) Model, etc. are extensively employed in research [14,15,16]. In particular, the PMF method is recommended by the US Environmental Protection Agency (EPA) for source contribution estimation [17]. This approach involves the quantitative calculation of the contributions of potential sources to soil heavy metal contamination at each data point, which performs excellently in practical applications [18].

Jiyuan’s mining industry is renowned for its remarkable scale of lead and zinc smelting, which stands as the largest in Asia while also achieving a prominent position in national production. As a prominent hub for high-quality and specialized steel in the central and western regions, Jiyuan boasts an assemblage of over 160 steel and processing enterprises, encompassing a vast array of bar and wire production capabilities. Additionally, Jiyuan holds significant stature as a pivotal silver production base within China, showcasing unparalleled output levels. The advancement of the mining industry has contributed to economic growth; however, it has also led to significant heavy metal pollution in the soil environment. A watershed is a relatively independent natural geographic system, which is the main source of accumulation and integration of heavy metals through rivers or highlands. It exhibits heightened susceptibility to heavy metal pollution. Taking the Manghe River watershed in Jiyuan City as the study area, this study conducted a comprehensive analysis of the spatial pattern and sources of eight soil’s heavy metal contaminants (Cr, Hg, As, Pb, Ni, Cd, Cu, Zn) by employing geostatistical interpolation and the PMF Model. The main objectives of this study are as follows: (1) to explore the content and enrichment characteristics of heavy metals in surface soils by applying statistical analysis and pollution models; (2) to reveal spatial heterogeneity of heavy metal pollution by using geostatistical analysis; and (3) to quantitatively estimate the potential sources of heavy metals and their contributions by employing the PMF Model. The findings will provide decision support for addressing issues related to soil pollution control, restoration, and farmland protection in the Manghe River watershed in Jiyuan City.

2. Materials and Methods

2.1. Study Area

The Manghe River watershed is situated in the northern part of Jiyuan City, nestled within the southern foothills of the Taihang Mountains and adjacent to the bank of the Yellow River. Its geographical coordinate is 112°23′37″~112°33′2″ E, 35°3′4″~35°9′40″ N (Figure 1). The total area of the watershed is about 112.82 km², with a population of about 150,000 people. The terrain is ‘horseshoe’ shaped, with a gradual slope from north-west to south-east, transitioning from mid-to-high mountains to hills, ridges and plains. The climate is warm temperate monsoon, with a stable average annual temperature of 14.5 °C and average annual rainfall of about 567 mm. The soils are mainly brunisolic soil and fluvo-aquic soils, with a pH of 6.07–8.16. Furthermore, this region experiences four distinct seasons along with ample sunlight exposure and abundant heat and water resources, providing unique favorable conditions for both industrial and agricultural development. After years of development, the industrialization level in the region has surpassed 70%, establishing itself as a leading smelting center for lead and zinc as well as a prominent silver production hub in China. The region contributes a quarter of China’s total lead production and 12% of its silver production, while also serving as a pivotal hub for steel, energy, chemicals and equipment manufacturing. However, as the first national demonstration zone for industrial and urban integration, the soil heavy metal pollution problem is particularly prominent.

2.2. Soil Sampling and Laboratory Analysis

A sampling grid of 1 km × 1 km was established in the watershed based on remote sensing and land change survey data. In the concentrated mining industry region, a finer grid of 0.5 km × 0.5 km was utilized to enhance the interpretation. The final sampling points were optimized by considering various factors such as transportation, water, residential areas, industrial enterprises, terrain and field investigation results. Thus, we obtained a total of 180 sampling points, including 59 encrypted sampling points. Soil sampling was conducted in March 2023, using handheld GPS for accurate positioning. Five samples from top layer (0–20 cm) were collected using the plum blossom sampling method, and a 1 kg mixed sample was obtained for experimental purposes following the ‘four-part method’. Simultaneously, an attribute table was completed to record surrounding environmental information. After impurity removal, drying and grinding, the sample was homogenized and sieved through 100 mesh nylon for storage.

The sample underwent microwave digestion with a combination of HNO₃, HF, HCL and HCLO₄ prior to analysis. Cr, Cu, Zn and Ni were quantified using flame atomic absorption spectrophotometry, Pb and Cd were determined using graphite furnace atomic absorption spectrophotometry, Hg was measured via cold atomic absorption spectrophotometry and As was analyzed using atomic fluorescence method. All measurements were conducted in triplicate, with the standard deviation remaining within ±5% of the mean value. The quality assurance and quality control (QA/QC) program was conducted using the certified soil reference material GBW07403 (GSS-3) provided by the National Center for Standard Reference Materials of China. The results demonstrated a relative standard deviation (RSD) of less than 5% and a recovery rate within ±10%, indicating that the method employed for detecting soil heavy metals was relatively accurate and reliable. The detection limits of Cr, Hg, As, Pb, Ni, Cd, Cu and Zn were 2 mg/kg, 0.0003 mg/kg, 0.05 mg/kg, 1 mg/kg, 1 mg/kg, 0.02 mg/kg, 1 mg/kg and 2 mg/kg, respectively. The Grubbs method was employed for outlier identification and removal resulting in retention of data from 165 soil samples.

2.3. Enrichment Factors

Enrichment factor (EF), a method proposed by Zoller et al. [19] in their study of the sources of chemical elements in pollutants over Antarctica, was used to assess the impact of anthropogenic activities on soil heavy metal concentrations. The calculation formula is as below:

$$EF=\frac{(Ci/Cref)sample}{(Bi/Bref)background}$$

(1)

where sample $\left(Ci/Cref\right) sample$ is the ratio of target metal to reference metal in soil samples and $\left(Bi/Bref\right)background$ is the ratio of target metal to reference metal in the background. Elements that are more stable and less affected by human activities (e.g., Fe, Mn and Al, etc.) are usually chosen as reference elements [20,21]. In this study, Fe was chosen as the reference element. According to reference [22], $EF$ was categorized into six levels: not enriched ($EF$ < 1), mildly enriched (1 < $EF$ < 2), moderately enriched (2 < $EF$ < 5), highly enriched (5 < $EF$ < 20), intensely enriched (20 < $EF$ < 40) and very highly enriched ($EF$ > 40).

2.4. Positive Matrix Factorization Model

PMF is a receptor modeling technique that utilizes sample components or fingerprints to quantify pollution sources. The model decomposes the original matrix $X_{ij}$ into two factor matrices $g_{ik}$ and $f_{kj}$ and a residual matrix $e_{ij}$. The basic equation is as follows:

$$X_{ij}={\sum }_{k=1}^p{g}_{ik}{f}_{kj}+e_{ij}$$

(2)

where $X_{ij}$ represents the concentration of the j th element in the I th sample; $g_{ik}$ is the contribution of source $k$ to sample $i$; $f_{kj}$ is the concentration of the $j$th element in source $k$; and $e_{ij}$ is the residual matrix.

The model is calculated iteratively by the method of least squares in order to minimize the objective function $Q$. The value of $Q$ is calculated as follows:

$$Q={\sum }_{i=1}^n{\sum }_{j=1}^m{(\frac{e_{ij}}{u_{ij}})}^2$$

(3)

where $u_{ij}$ is the uncertainty of the $j$ th element in the sample $i$.

The calculation of uncertainty depends on the relationship between the detection quantity and the detection limit, and the formula is as follows:

$$U_{\mathrm{ij}}=\left\{\begin{array}{c}\sqrt{(Y_{RSD}{X}_{ij}{)}^2+B_{MDL}} C>B_{MDL} \\ \frac{5}{6}\times B_{MDL} C\le B_{MDL} \end{array}\right.$$

(4)

where C is the elemental concentration; Y_RSD is the relative standard deviation; and B_MDL is the method’s limit of detection with a 5% margin of error.

2.5. Statistics and Processing of Data

One-way analysis of variance (ANOVA) was used to compare the significant differences in soil’s heavy metal content under different land use types. Pearson’s correlation analysis was employed to assess the relationship between heavy metals. Additionally, the PMF Model was utilized to dissect the sources of each heavy metal and determine their respective contributions. Prior to conducting source apportionment using PMF 5.0, a three-fold standard deviation method was applied to test and exclude any values that could potentially impact the results. The Kolmogorov–Smirnov (K-S) test was employed to assess the normal distribution of heavy metal content before conducting the one-way ANOVA and geostatistical analysis. Data deviating from normality were subjected to logarithmic transformation and normalization. Statistical analyses, Spearman’s correlation analysis and one-way ANOVA were performed in IBM SPSS Statistics 27, a semi-variogram model was used in GS+ 9.0, and ordinary Kriging interpolation was performed in ArcGIS 10.8.

3. Results

3.1. Heavy Metal Concentrations

Table 1. Descriptive statistics of heavy metals and their EF (mg/kg).

	Cd	Hg	Pb	As	Cu	Zn	Cr	Ni
	mg/kg	mg/kg	mg/kg	mg/kg	mg/kg	mg/kg	mg/kg	mg/kg
Mean	2.59	0.131	121	19.5	33.9	137	56.1	27.6
Maximum	16.2	2.54	924	109	179	534	87.1	49.6
Minimum	0.551	0	0	0.301	8.43	32.2	21.2	9.87
Standard deviation	2.03	0.361	150	13.5	20.7	79.1	8.74	5.53
Coefficient of variation	0.781	1.48	1.24	0.691	0.610	0.581	0.161	0.201
Background of Hena 1	0.0741	0.0340	19.6	11.4	19.7	60.1	63.8	29.9
Average of China 2	0.0971	0.0650	27.1	11.1	23.1	74.1	61.1	27.1
Chinese soil criteria (Grade II) 3	0.601	1.01	240	20.1	200	300	350	190
Enrichment factor (n = 165)	35.0	7.10	6.19	1.71	1.72	2.27	0.881	0.921

¹ Data came from The Chinese Environmental Monitoring Centre (CNEMC, 1990); ² Data came from Teng [23].; ³ Data came from National Environmental Protection Agency of China (CNEPA, 1995).

shows the descriptive statistics and EF results for the eight metals in the topsoil.

The average content of each element in descending order was Hg (0.13 mg/kg) < Cd (2.59 mg/kg) < As (19.46 mg/kg) < Ni (27.60 mg/kg) < Cu (33.89 mg/kg) < Cr (56.02 mg/kg) < Pb (121.24 mg/kg) < Zn (136.66 mg/kg). The contents of Ni and Cr met the criteria for normal distribution (p > 0.05), while those of Cd, As, Cu, Pb and Zn met the criteria for log-normal distribution. According to the CV values, elements Cr (15.60%) and Ni (20.05%) showed low variability. In contrast, elements Zn (57.84%), Cu (61.09%), As (69.32%) and Cd (78.35%) demonstrated medium variability. Notably, elements Pb (123.71%) and Hg (147.70%) displayed high variability.

When compared with the Chinese National Standard for Soil Environmental Quality (CNEPA, 1995), we found that the average contents of other heavy metals fell below their corresponding risk screening values, except for Cd. However, it remained lower than the risk control value, indicating a potential soil pollution risk [24]. In addition, the mean concentrations of Cr and Ni were lower than their corresponding background values (CNEMC, 1990), while the contents of Hg, Cd, As, Cu, Pb and Zn exceeded their respective background values. It was noteworthy that the mean contents of Cd, Hg and Pb were 26.70, 3.69 and 4.49 times higher than those found in Chinese soils, respectively [25]. The average EF for heavy metals in the study area followed this order: Cr < Ni < As < Cu < Zn < Pb < Hg < Cd. The average EF of Cr and Ni did not exhibit any significant enrichment, whereas the remaining heavy metals demonstrated varying degrees of enrichment. Notably, the EF values for Cd and Hg were remarkably high at 218 and 74.7, respectively, indicating strong enrichment characteristics.

3.2. Spatial Distribution of Heavy Metals

Statistical analysis indicated that the concentrations of Cr and Ni in soil followed a normal distribution, whereas those of Hg, As, Cu, Pb, Zn and Cd followed a log-normal distribution. The exponential model was used to describe Cr, Ni and As and the spherical model was employed for Hg and Cu, while the Gaussian model provided the best fit for the Pb, Zn and Cd data. The r² values of these models ranged from 0.731 to 0.997, effectively elucidating the spatial structural characteristics of soil heavy metals. The spatial distribution of each heavy metal content was depicted in Figure 2 using ordinary kriging interpolation.

(1) The spatial distribution of Cr, Ni and Cu exhibited a similar pattern, characterized by a gradual decrease in concentration from the central high-value area towards the east and west. The high-value regions for Cr and Ni were oriented in a north-west–south-east direction, while those for Cu were oriented in a north-east–south-west direction. (2) The distribution patterns of Zn and As exhibited similarities, characterized by high-value centers situated in the north-west, north-east, and south-east, with values decreasing outward. The irregular decrease in the shape of Zn contrasted with the regular decrease observed in As. (3) The spatial distribution of Cd, Hg and Pb also exhibited some similarities, characterized by high-value regions extending in a north-east–south-west direction and displaying a decreasing concentration gradient. However, their distributions differ. The high-value areas for Cd were located in the south-east, while those for Hg were situated in the south-west and central regions, with the highest concentration found in the south-west corner; the high-value areas for Pb were predominantly concentrated in the north-east and central regions.

3.3. Relationship between Heavy Metal Content and Environmental Factors

The distribution of land use in the Manghe River watershed is illustrated in Figure 1. After categorizing and analyzing metal content across land uses, we ran statistical tests like Tamhane’s T2 and Tukey’s variance test for one-factor variance analysis, as detailed in

Table 2. Results of ANOVA for the mean heavy metal contents under different land uses (mg/kg).

		Ns	Hg	As	Cr	Cu	Ni	Pb	Zn	Cd
Nature	Forest	41	0.193	10.7	67.9	24.1	33.9	88.2	92.79	1.62
	Orchard	35	0.112	12.4	66.9	24.9	32.1	88.1	94.42	1.78
Farmland	Dryland	24	0.176	17.0	56.6	28.8	29.9	99.7	137	2.05
	Paddy	17	0.154	16.9	56.1	28.3	28.5	98.0	135	1.94
Constructions	Residential	8	0.179	20.2	49.8	35.8	22.7	113.1	147	2.69
	Industrial and mining	13	0.237	23.2	51.8	39.0	24.1	126.73	152	2.74
	Transportation	11	0.211	20.7	50.4	40.1	23.7	117.36	158	2.91
	Other construction	16	0.172	20.4	52.2	38.3	23.3	115.78	149	2.83
F			8.61 **	2.18 *	1.13	2.26 *	1.30	3.25 **	2.02 *	1.29 **

*. Significant at the 0.05 level; **. Significant at the 0.01 level.

.

Different land use types varied in heavy metal content. As, Cu, Pb, Zn and Cd levels in construction land exceeded natural vegetation and cultivated land. Residential areas had lower pollution than industrial/mining sites. Transportation facilities accumulated Cu, Zn and Cd. Dryland had higher pollution than paddy fields due to soil erosion. Forests and orchards showed varying accumulation patterns.

3.4. Source Apportionment

3.4.1. Correlation Analysis

The Spearman correlation coefficients of the eight heavy metals in soil are presented in

Table 3. Spearman correlations matrix for the heavy metal contents in topsoil.

	Hg	As	Cr	Cu	Ni	Pb	Zn	Cd
Hg	1.00
As	0.16 *	1.00
Cr	−0.05	0.16 *	1.00
Cu	0.32 **	0.63 **	0.12	1.00
Ni	−0.06	0.20 **	0.83 **	0.14	1.00
Pb	0.25 **	0.67 **	0.08	0.61 **	0.11	1.00
Zn	0.21 **	0.46 **	0.03	0.52 **	0.02	0.60 **	1.00
Cd	0.13	0.61 **	0.05	0.46 **	0.04	0.59 **	0.45 **	1.00

*. Significant at the 0.05 level; **. Significant at the 0.01 level.

. A significant positive correlation (p < 0.01) was observed between Cr and Ni (r = 0.830), As and Pb (r = 0.673), As and Cu (r = 0.633) and As and Cd (r = 0.612), as well as As and Zn (r = 0.459), indicating possible homology among these elements. Additionally, strong correlations were observed among Pb, Cu, Cd and Zn, but specific sources need further investigation.

3.4.2. Positive Matrix Factorization (PMF)

Initially, 165 samples’ content data for eight heavy metals and related uncertainty concentration data served as input files for the PMF Model. The best FOUR-factor solution was identified based on minimum and stable Q-value. Notably, most residuals fell within the −3 to 3 range at this stage, indicating a good fit to the data. Meanwhile, all heavy metals exhibited signal-to-noise ratios (S/N) greater than 2, demonstrating strong explanatory power and model validation. Figure 3 shows the contribution of the four-factor components extracted from the PMF Model to the load of each heavy metal. Based on the factor fingerprints of each heavy metal, we accurately quantified the respective contributions of individual pollutant sources to heavy metal contamination. The following could be drawn:

(1) Factor 1 contributed the most to Hg (60.7%) and Cu (50.2%). Hg and Cu displayed a strong positive correlation (r = 0.317), indicating a common source. Coal combustion, specifically Chinese coal with high Hg levels (~0.17 mg/kg), was a major Hg release pathway [26,27]. Furthermore, Cu accumulation may be attributed to its usage in pesticides, herbicides and fertilizers in agriculture [28]. Notably, Cu, a common livestock feed additive, was strongly linked with livestock manure, especially pig manure [29]. Despite the fact that Factor 1 contributed only 31.6% to Zn, Zn showed a strong positive correlation with Cu (r = 0.52). In fact, Zn was also a major additive in agricultural fertilizers, so it could be concluded that Factor 1 might be associated with industrial and agricultural activities.

(2) Factor 2 contributed the most to Cr (74.2%) and Ni (73.8%). Ni and Cr were usually considered as indicators of natural sources, predominantly originating from the soil matrices [30]. Hence, it could be inferred that Factor 2 was of natural source.

(3) Factor 3 contributed the most to Pb (88%), As (33.2%) and Cd (27.8%). Transport emissions, such as fuel combustion, engines and catalyst use, were major Pb sources [31]. The correlation between As and Cd was found to be strong, with a correlation coefficient of 0.612, and they were closely related to agricultural activities. Phosphorus content in the study area averaged 0.741 g/kg, peaking at 3.23 g/kg. This met moderate richness criteria set by the Second National Soil Survey. Thus, excessive phosphorus fertilizer usage may lead to As and Cd accumulation [32]. Furthermore, the extensive use of inorganic arsenic compounds in pesticides and herbicides has also contributed to As accumulation in soil [33]. Therefore, Factor 3 is inferred as a mixed source of traffic and agriculture.

(4) Factor 4 contributed the most to Zn (62.3%) and Cd (34%). Correlation analysis showed a positive link between Zn and Cd (r = 0.452), suggesting a common source. Zn and Cd were commonly found in minerals as sulphides and disulphides, which were closely associated with mining activities [34,35]. Additionally, industrial processes such as electroplating, metal smelting, and the chemical industry were the main contributors to Cd pollution; however, under certain conditions, Cd can also enter the soil via fly ash generated from coal combustion or iron ore processing [36]. Therefore, it can be concluded that Factor 4 might be strongly associated with industrial activities.

(5) The contribution of the four factors in descending order was Factor 2 (36.25%) > Factor 1 (23.00%) > Factor 3 (21.20%) > Factor 4 (19.55%). Anthropogenic activities and natural factors accounted for 63.75% and 36.25%, respectively.

4. Discussion

4.1. Spatial Distribution of Heavy Metals and Their Influencing Factors

Soil serves as an environmental receptor for heavy metals. The surface soil (0–20 cm) is a crucial component of the soil system and exhibits the highest sensitivity to environmental changes, such as temperature fluctuations, precipitation variations, and SOM accumulation. This heightened sensitivity results in significant heterogeneity in soil’s heavy metal accumulation. Therefore, mixed soil samples from this depth for heavy metal pollution analysis and source appointment are representative. In the present study, higher concentrations of heavy metals were observed in constructed land, particularly in industrial and mining areas as well as transportation zones, which is consistent with the findings reported by Liu et al and Li et al. [3,13]. Furthermore, previous studies have demonstrated severe Pb contamination in the sediment of the Shi River (upstream of the Manghe River Basin) near the smelter within this study area [37], indicating that mining and smelting operations significantly affect water quality in nearby rivers due to topographic relief and runoff [38]. Drylands, due to their steeper slopes, are more prone to soil erosion than paddy. This causes heavy metal accumulation and higher contamination levels [39]. Autumn–winter winds impact forested areas located downwind from industrial and mining facilities, causing higher heavy element concentrations via atmospheric deposition [40]. In summary, heavy metal contamination of soils in small watershed is influenced by a combination of natural and anthropogenic factors, including topography, soil erosion, vegetation cover, mining activities, industrial discharges and agricultural practices.

4.2. Appointment of Heavy Metal Sources

When two elements correlate, there are several possible reasons: they have similar soil-forming matrices; they have similar chemical properties and enrichment pathways; they have the same sources; or they come from different sources but operate in the same spatial regions. The correlation coefficient for Cr and Ni is 0.83, while those among As, Cu, Pb, Zn and Cd are 0.46–0.63, indicating natural sources and anthropogenic activity sources, respectively. The results of this paper are in agreement with the findings from stable isotope measurements of Pb, validating the applicability of the PMF in analyzing sources of heavy metal contamination in soils [41,42]. Compared to alternative identification methods, the PMF Model enables accurate quantification of the contribution rate from each source without requiring prior knowledge of the source composition spectrum. Moreover, this model effectively handles missing and imprecise data and avoids negative values in both source composition spectra and contribution rates by imposing non-negative constraints on factor decomposition matrices [43,44]. The combination of the PMF Model and the geostatistical method provides a practical way of tracing pollution sources and targeted solutions for managing regional soil’s heavy metal pollution.

4.3. Countermeasures against Heavy Metal Pollution

Soil’s heavy metal pollution accumulates rapidly within the Manghe River watershed; therefore, there is an urgent need to implement effective measures for stringent control and reductions in the input of heavy metal pollutants into the soil and water environment, aiming to safeguard both the soil ecosystem and human health and safety. On the one hand, a scientifically based zoning approach will be employed to ensure safe use of construction land and agricultural land in the study area, taking into account the extent of heavy metal contamination. Emphasis will be placed on enhancing protection measures for lightly contaminated areas, while implementing soil remediation strategies for areas with moderate to severe contamination. On the other hand, a differentiated control strategy including source reduction, process control, and final treatment is proposed based on an analysis of spatial distribution driving factors as well as pollution sources.

4.4. Shortcomings and Prospects

To characterize the spatial distribution of heavy metal concentrations, this study selected the semi-variogram function with the highest coefficient of determination as the best-fit model and used ordinary kriging interpolation to simulate the spatial concentration distribution. However, due to their trace and transport properties (physical transport and redox reactions), soil heavy metals are usually associated with soil organic matter and ferromanganese compounds, exhibiting significant non-linear relationships with covariates such as topography and climate. Consequently, linear fitting in kriging interpolation is inadequate. Therefore, establishing effective predictive relationships incorporating multiple covariates is crucial for achieving high precision interpolation and represents a key focus of future research in this thesis. In terms of quantitative resolution, the PMF Model used in this study demonstrated superior capability in quantifying the contribution of each source. However, subjective judgments regarding specific sources of heavy metals could introduce errors in traceability results. To address this limitation, future research on soil heavy metal source appointment could focus on integrating multiple methods and improving quantitative assessment of pollutant source contributions. Moreover, a single sampling depth cannot fully reveal the response relationship between soil heavy metals and external factors. This not only negatively impacts the precision of interpolation simulation but also leads to insufficient information for source apportionment. Therefore, it is recommended to subdivide the 0–20 cm soil layer into two layers: the surface layer (0–10 cm) and the subsurface layer (10–20 cm). By distinguishing the content differences between these two layers, we can more accurately elucidate the correlation between heavy metals and other environmental factors, as well as migration and transformation information, which could improve the accuracy of spatial pattern simulation and source apportionment.

5. Conclusions

The contamination level, spatial distribution, and potential sources of eight soil heavy metals in the Python River sub-watershed of Jiyuan City were systematically analyzed in this study using enrichment factors, geostatistics, correlation analyses and PMF Models. The main findings are as follows.

(1) The descending order of the average heavy metals concentrations was Hg (0.13 mg/kg) < Cd (2.59 mg/kg) < As (19.46 mg/kg) < Ni (27.60 mg/kg) < Cu (33.89 mg/kg) < Cr (56.02 mg/kg) < Pb (121.24 mg/kg) < Zn (136.66 mg/kg). The average concentrations of Cr and Ni were lower than background values, exhibiting no significant enrichment and low spatial variability; however, the other six heavy metals exceeded background values, with Zn, Cu, As and Cd showing medium variability, while Hg and Pb showed high variability. Notably, the average contents of Cd, Hg and Pb were 26.70, 3.69 and 4.49 times higher than those in Chinese soils on average, respectively, indicating a highly significant enrichment.

(2) The concentrations of Cr and Ni in the soil exhibited a normal distribution, while those of Hg, As, Cu, Pb, Zn and Cd followed a lognormal distribution. The semi-variogram analysis revealed that the nugget coefficients (C0/C0 + C) of different heavy metals ranked as follows in descending order: Hg (0.873) > Cd (0.824) > Pb (0.797) > Zn (0.679) > Cu (0.472) > As (0.422) > Ni (0.237) > Cr (0.225). The spatial distribution of heavy metals was influenced by various factors including topography, land use type, surface runoff, wind direction and anthropogenic activities, resulting in significant spatial heterogeneity.

(3) Cr exhibited a highly significant positive correlation with Ni, while As demonstrated a significant positive correlation with Pb, Cu, Cd and Zn. Additionally, there was a strong correlation among Pb, Cu, Cd and Zn. The contribution of anthropogenic activities to heavy metal pollution in the study area accounted for 63.75%,whereas natural factors contributed 36.25%. Mining and related chemical smelting activities, vehicle exhaust emissions, and improper use of fertilizers and pesticides in agriculture were identified as the primary sources of heavy metal pollution in this region. These pollutants enter the soil environment through various pathways such as direct emissions, atmospheric deposition, transport mechanisms and agricultural practices, consequently resulting in substantial adverse impacts on soil ecology.