Heavy Metal Pollution and Its Prior Pollution Source Identification in Agricultural Soil: A Case Study in the Qianguo Irrigation District, Northeast China

Abstract

Heavy metals are the primary pollutants in agricultural soil and have hindered the sustainable development of agriculture. To control heavy metal pollution, it is essential to identify the pollution sources, particularly the prior source, in agricultural soils. In the current study, Qianguo Irrigation District, a typical agricultural region in Northeast China, was selected to be investigated for the source apportionment of soil heavy metals and identify the prior pollution source. The results showed that the study area was at a moderate pollution level with considerable ecological risk, while Hg and Cd were the main pollutants. Human-health risk assessment indicated that the non-carcinogenic risk for all populations was acceptable (HI < 1), and the carcinogenic risk was not negligible (10⁻⁶ < TCR < 10⁻⁴). The main pollution sources were concluded to be of lithogenic origin (35.5%), livestock manure (25.4%), coal combustion (21.5%), and chemical fertilizers (17.7%). Coal combustion was identified as the prior pollution source, accounting for 47.69% of the RI contribution. This study can provide scientific support for environmental management and pollution control of soil heavy metals in agricultural regions.

1. Introduction

Agricultural soils are facing tremendous production pressure with population growth and economic development [1,2,3]. It was reported that global agricultural production must double between 2005 and 2050 to fulfill human needs [3]. However, intensive anthropogenic activities have caused worldwide agricultural soil pollution, hindering the sustainable development of agriculture [4,5,6]. According to the Bulletin of National Soil Pollution Survey in China, the total exceedance rate of agricultural soil nationwide was 19.4% of which heavy metal pollution accounted for 82.4% [7]. Heavy metals (HMs) are revealed as the most important contaminants in agricultural soils due to their ubiquity, high toxicity, and bioavailability, and they may pose a serious threat to agroecosystem and human health [3,8,9]. HMs in agricultural soils may not only reduce the quality and yield of crops but also cause pathogenic and carcinogenic risks to humans via long-term exposure through ingestion, inhalation, dermal contact, and food chains [10,11,12,13,14]. Consequently, it is of great significance to evaluate pollution characteristics and environmental risk and, determine the sources of HMs in agricultural soils for the sustainable development of agriculture and human health.

Currently, many researchers have developed various methods to evaluate the pollution of HMs in agricultural soils. Classical assessment indices, such as contamination factor, geo-accumulation index, enrichment factor, pollution load index, Nemerow pollution index, improved Nemerow index, and risk assessment code, have been widely applied to assess the pollution levels of HMs in soil [15,16,17]. Meanwhile, considering the variable toxicity of different HMs, risk assessment methods, including potential ecological risk index and human-health risk assessment, are applied to evaluate the risk of soil HMs [11]. In this study, a combination of various methods was employed to acquire the comprehensive pollution status of HMs in agricultural soils.

To effectively control and reduce HM pollution, it is crucial to conduct source apportionment of HMs [4,18,19]. In recent years, the identification and quantification of pollution sources were collectively referred to as source apportionment [20]. Multivariate statistical analyses (correlation analysis and principal component analysis) are excellent tools for the identification of pollution sources by soil HMs [16,21,22]. Positive matrix factorization (PMF), an optimal receptor model, can quantify the source contribution of HMs at each site under non-negativity constraints and has been widely used to quantitatively determine the sources of HMs in soil [23,24,25]. However, multivariate statistical analysis cannot quantify the sources of HMs in soil, and the source identification by PMF model exhibits some uncertainty [17]. The combination of the PMF model and multivariate statistical analyses will improve the accuracy of source apportionment for soil HMs. Furthermore, previous studies focused on quantitative identification of the concentrations for different pollution sources, but ignored the toxicity differences of HMs [11,25]. The sources with the high contributions did not necessarily have high environmental risks. To identify prior pollution sources, the corresponding environmental risks should be considered [26,27]. Thus, an integrated method was developed based on the combination of pollution source information and ecological risk assessment in the current study, which will provide accurate information on the prior source of HMs in agricultural soils.

There is widespread concern on HM pollution of agricultural soils in China. However, previous studies mainly focused on the developed central and eastern coastal areas [28,29]. Northeast China, as the main commercial grain base, provides 40% of the national commercial grain, which is also facing the enrichment of HMs in agricultural soils. The Qianguo Irrigation District, located in Northeast China, is a main commercial grain production area exposed to various sources of HMs due to intense industrial and agricultural activities. Thus, this region was selected to (1) investigate the distribution of HMs concentrations; (2) comprehensively assess soil HM pollution; (3) quantitatively identify the pollution sources of soil HMs and their prior source in this region.

2. Materials and Methods

2.1. Research Area

The Qianguo Irrigation District is located in Songyuan City on the left bank of the Songhua River in Jilin Province, Northeastern China (Figure 1). The region is in the temperate continental semi-arid monsoon zone, with average annual precipitation and evaporation of 400–500 mm and 1500 mm, respectively. The special climate conditions lead to the invasion of soil salinization. The local soil type is mainly light chernozem and paddy soil, and the soil parent material is loess sediment.

The Qianguo Irrigation District is a key commercial grain production base in China. It was firstly cultivated in 1944, with a total area of about 716 km² of which paddy fields account for 2/3, i.e., 433 km². To guarantee grain production, a great number of chemical fertilizers (urea, phosphate, potassium sulfate, zinc sulfate, etc.) have been applied to the arable land since the late 1990s.

2.2. Sampling and Analysis

A total of 35 soil surface samples (0–15 cm) were collected in early May 2019, the sampling points were located using GPS (Figure 1), and the details on location and the reclamation years are presented in Supplementary Materials Table S1. Each sample was a mixture of five subsamples from 5 m² around the sample site. All samples were collected with a wooden spatula and kept in polyethylene bags.

Soil was dried naturally in the laboratory after removal of plant residues and gravels. The dry samples were ground and passed through a 100–mesh nylon sieve. Hereafter, the samples (0.1 g) were digested with a mixture of HCl–HNO₃–HF–HClO₄. In addition, the available fractions of HMs were extracted into the solution by a mixture of DTPA–CaCl₂–TEA. The contents of V, Cr, Co, Cu, Zn, Cd, Pb, Hg, and As were measured using inductively coupled plasma-mass spectrometer (Agilent 7500, Santa Clara, CA, USA) and atomic fluorescence spectrometry (AFS–933, Beijing, China). Soil pH was measured by potentiometry with a ratio of 1:2.5 (soil: water) using a pH meter (PHS–3G, Shanghai, China).

For quality control and quality assurance, the reference material (GBW07408) was determined using the same method, and the recovery rates were between 89.1% and 108.8% with a relative standard deviation lower than 10%.

2.3. Pollution Level Assessment

Contamination factor (CF), geo-accumulation index (I_geo), enrichment factor (EF), pollution load index (PLI), Nemerow pollution index (P_N), and improved Nemerow index (INI) were used to evaluate the pollution levels of HMs in soil. The calculation methods and the grade classification criteria (Supplementary Materials Table S2) are shown in the Supplementary Materials.

2.4. Ecological Risk Assessment

Both the risk assessment code (RAC) and potential ecological risk index (RI) were applied for evaluating the ecological risks of soil HMs [30,31]. Their equations are shown in the Supplementary Materials, and the grade classifications are shown in Supplementary Materials Table S3.

2.5. Human-Health Risk Assessment

Human-health risk assessment (HHRA) can be applied to evaluate the carcinogenic and non-carcinogenic risks (CR and NCR) of soil HMs to all populations (children and adults) through three pathways: ingestion and dermal and inhalation absorption [11,32]. The equations are shown in the Supplementary Materials, and the corresponding parameters and grade classifications are presented in Supplementary Materials Tables S4–S6.

2.6. Quantitative Source Apportionment

Spearman correlation analysis (CA) and principal component analysis (PCA) were conducted using SPSS v23.0 (IBM, USA). PMF model was applied to obtain quantitative source apportionment, which divided the HM concentration matrix into two matrices, i.e., factor contribution and factor profile. The specific principles and formulas are presented in the Supplementary Materials.

2.7. Source-Oriented Ecological Risk Assessment

An integrated method was established to quantify source-oriented risk assessment of soil heavy metals. For each sample, the HM contents of different sources were obtained according to the PMF model, and then RI was calculated based on the HM contents of each source; thus, the ecological risk of each source was evaluated. The source with the highest RI value would be the prior source of HMs. The calculation formulas are as follows:

$$E_r^i{=T}_r^i \times \frac{C_{\mathit{ij}}}{C_b}$$

(1)

$$RI={\displaystyle \sum _{i=1}^n{E}_r^i}$$

(2)

where C_ij is the concentration of HM i for factor j, C_b is the background value for the corresponding element, and the other parameters are consistent with the RI (Supplementary Materials).

3. Results and Discussion

3.1. HM Concentrations

The statistical analysis of soil pH and HM concentrations in Qianguo Irrigation Districts are presented in

Table 1. Soil pH and heavy-metal concentrations.

Parameter	pH	Heavy-Metal Concentrations in the Soil (mg/kg)
		V	Cr	Co	Cu	Zn	Cd	Pb	Hg	As
Mean	7.89	50.97	28.42	8.05	12.05	53.25	0.28	17.55	0.17	6.36
Minimum	6.86	32.26	17.06	5.71	6.96	31.69	0.10	15.11	0.05	4.29
Maximum	8.85	72.56	63.83	12.98	19.76	93.39	0.64	21.81	0.26	10.75
SD a	0.45	8.95	10.77	1.75	3.10	14.31	0.10	1.69	0.06	1.30
CV b		0.18	0.38	0.22	0.26	0.27	0.38	0.10	0.34	0.20
BV c	8.20	60.50	34.26	9.58	13.00	40.84	0.08	19.74	0.02	6.1
RSV d	pH > 7.5		250		100	300	0.6	170	3.4	25
Nexceed e (%)		17.14	17.14	17.14	34.29	77.14	100.00	11.43	100.00	57.14

^a: standard deviation; ^b: coefficient of variation; ^c: the background values of soil environmental in Jilin province [34]; ^d: the screening value of soil-pollution risk of agricultural land [33]; ^e: the rate of the number of soil samples exceeding the background value to the total number of samples.

. The mean value of soil pH was 7.89 with an overall alkaline appearance. The mean contents of V, Cr, Co, Cu, Zn, Cd, Pb, Hg, and As were 50.92, 28.42, 8.05, 12.05, 53.25, 0.28, 17.55, 0.17, and 6.36 mg/kg, respectively, which were lower than the corresponding risk-screening values (RSVs) [33]. Nevertheless, the mean concentrations of Zn, Cd, Hg, and As surpassed the local background concentrations by a factor of 1.40, 3.34, 6.99, and 1.04, respectively [34]. Compared with the median values of HMs in European agricultural soils, the mean contents of Cr, Co, Zn, Pb, Hg, and As in the study area were higher, while those of Cd and Cu were lower [35,36]. The exceedance rate (N_exceed) was ranked as Hg (100%) = Cd (100%) > Zn (77.14%) > As (57.14%) > Cu (34.29%) > V (17.14%) = Cr (17.14%) = Co (17.14%) > Pb (11.43%), indicating that Hg, Cd, Zn, and As had been highly enriched in local agricultural soil. The coefficients of variation (CV) for Cd (0.38) and Hg (0.34) exceeded 0.3, indicating higher spatial variability. Based on their high concentrations, it was suggested that the concentrations of Cd and Hg had been mainly influenced by anthropogenic sources [20].

3.2. Pollution Levels of Soil HMs

The results of CF are shown in Figure 2, and the mean values of Zn, Cd, Hg, and As were 1.30, 3.34, 6.99, and 1.04, respectively. Hg and Cd were at high and considerable contamination levels, respectively. Zn and As were at moderate contamination levels. The PLI_n values (Supplementary Materials Table S7) of all samples were greater than 1, and the PLI_zone value was 1.31, which indicated that the study area had been contaminated. The P_N values of all samples (Supplementary Materials Table S8) were greater than 3, indicating that the study area was at a serious contamination level. In general, the study area had been seriously contaminated, with Hg and Cd as the most polluting HMs.

I_geo, EF, and INI were utilized to assess the pollution levels of soil HMs by anthropogenic influence. The results of I_geo are shown in Figure 2. For Hg, 62.9% and 33.3% of the samples were moderately to heavily polluted and moderately polluted, respectively, 68.6% moderately or moderately to heavily polluted for Cd, and 28.6% unpolluted to moderately polluted for Zn. The other HMs were at the unpolluted level. Similarly, EF results (Figure 2) indicated that 60.0% and 20.0% of Hg in all samples were at the level of moderate and significant enrichment, respectively, and 37.1% of Cd were moderate enrichment. This suggested that Hg and Cd in soil exhibited the highest enrichment levels, and they were more likely to originate from human activities. The results of INI showed that the study area was moderately contaminated, with 91.4% of the samples at moderately contaminated levels (Supplementary Materials Table S9). The spatial distribution of INI (Figure 3) showed that southeastern part of the Qianguo Irrigation District exhibited the highest pollution level, which may be due to intensive industrial activities and the longer period of cultivation. Meanwhile, the concentrations of Cd and Zn increased significantly (p < 0.05) along with the reclamation year (Figure 4), which also proved that the longer cultivation period in the southeast of the study region resulted in the high enrichment of HMs [37,38].

3.3. Risk Assessment of HMs in Soil

3.3.1. Potential Ecological Risk

The results of RAC (Figure 5) indicated that Cu (19.12) and Cd (12.82) exhibited high values and were at a medium risk level. The relatively high proportion of available Cu and Cd indicated that they were likely to migrate into surrounding ecosystems or agricultural products, which may pose a high ecological risk [39,40]. The rest followed the order as Pb (7.65) > Zn (2.59) > As (1.60) > Co (1.49) > Cr (0.97) > V (0.58), all of which were at light and no risk levels. The average value of RI indicated a considerable ecological risk, with a mean value of 408.28, and most of the samples (77.1%) were at a considerable level (Supplementary Materials Table S10). Hg and Cd showed high and considerable potential risk, respectively, and they are the main contributors to RI, with a mean contribution rate of 67.2% and 25.5%, respectively. The spatial analysis (Figure 6) showed that the southeastern area exhibited the highest ecological risk level, which was in line with the results of INI assessment, with intensive industrial and agricultural activities leading to the higher ecological risk. Therefore, the southeast was the most polluted region in the study area with the highest ecological risk and needed to be controlled as a priority.

3.3.2. Human-Health Risk

The results of the HHRA are presented in

Table 2. The assessment of health risk for heavy metals in soil.

Parameter	Heavy Metals	Mean ± (SD)				Median
		Children		Adult		Children	Adult
HQ	V	5.86 × 10−2	1.01 × 10−2	1.56 × 10−2	2.70 × 10−3	5.73 × 10−2	1.53 × 10−2
	Cr	6.93 × 10−2	2.59 × 10−2	1.76 × 10−2	6.57 × 10−3	6.65 × 10−2	1.69 × 10−2
	Co	2.90 × 10−3	6.20 × 10−4	8.76 × 10−4	1.87 × 10−4	2.80 × 10−3	8.45 × 10−4
	Cu	1.99 × 10−3	5.04 × 10−4	4.68 × 10−4	1.19 × 10−4	1.87 × 10−3	4.41 × 10−4
	Zn	1.18 × 10−3	3.11 × 10−4	2.77 × 10−4	7.35 × 10−5	1.10 × 10−3	2.59 × 10−4
	Cd	2.24 × 10−3	8.34 × 10−4	6.00 × 10−4	2.23 × 10−4	2.16 × 10−3	5.78 × 10−4
	Pb	3.33 × 10−2	3.16 × 10−3	7.89 × 10−3	7.48 × 10−4	3.24 × 10−2	7.68 × 10−3
	Hg	3.78 × 10−3	1.25 × 10−3	9.06 × 10−4	3.00 × 10−4	4.03 × 10−3	9.66 × 10−4
	As	1.40 × 10−1	2.82 × 10−2	3.28 × 10−2	6.63 × 10−3	1.37 × 10−1	3.23 × 10−2
HI	Total	3.13 × 10−1	5.74 × 10−2	7.71 × 10−2	1.42 × 10−2	3.06 × 10−1	7.52 × 10−2
CR	Cr	8.72 × 10−6	3.26 × 10−6	9.46 × 10−6	3.53 × 10−6	8.36 × 10−6	9.07 × 10−6
	Cd	8.66 × 10−8	3.22 × 10−8	9.33 × 10−8	3.47 × 10−8	8.35 × 10−8	8.99 × 10−8
	As	5.39 × 10−6	1.09 × 10−6	5.49 × 10−6	1.11 × 10−6	5.30 × 10−6	5.41 × 10−6
	Co	1.22 × 10−9	2.61 × 10−10	5.32 × 10−9	1.14 × 10−9	1.18 × 10−9	5.13 × 10−9
	Pb	8.37 × 10−8	7.94 × 10−9	8.49 × 10−8	8.05 × 10−9	8.14 × 10−8	8.26 × 10−8
TCR	Total	1.43 × 10−5	3.86 × 10−6	1.51 × 10−5	4.14 × 10−6	1.38 × 10−5	1.47 × 10−5

HQ: hazard quotient; HI: hazard index; CR: carcinogenic risk; TCR: total carcinogenic risk (Supplementary Materials).

. For NCR, the mean hazard index (HI) values for all populations were lower than 1, indicating that these HMs posed little potential adverse effect on human health. Children were more susceptive due to the relatively high HI values in comparison with adults. As contributed 44.73% and 42.54% of the HI values for children and adults, respectively, which indicated that more attention should be given to As. For CR, the mean TCR values for children and adults were 1.43 × 10⁻⁵ and 1.51 × 10⁻⁵, respectively, with a range of 10⁻⁶ to 10⁻⁴, indicating a non-negligible carcinogenic risk for all populations. The higher TCR values in adults indicated that adults suffered a higher carcinogenic risk. The mean CR values for all populations followed the order as Cr > As > Cd > Pb > Co, and Cr and As were the main contributors of TCR. The mean CR values of Cr and As were greater than 10⁻⁶ for all populations, indicating that Cr and As were identified as the dominant factors in carcinogenic risk. Therefore, for human-health risks, the prior control elements were Cr and As.

3.4. Source Apportionment

3.4.1. Multivariate Statistical Analysis

The results of the Spearman correlation analysis (CA) of HM contents were calculated and are presented in Figure 7. Significant positive correlations among HMs were observed, indicating that they may originate from similar sources, while negative or nonsignificant correlations may represent different sources of HMs [11,17]. In the current study, a significant positive correlation was observed among all the HMs (p < 0.05) other than Cd, except for Hg-Cu and Hg-Cr, indicating that they may share a similar origin. The were no significant correlations between Cd and the other HMs, suggesting that Cd may have a unique source.

The results of the principal component analysis (PCA) are presented in Supplementary Materials Table S11. The results of the KMO test (0.800) and Bartlett star test (< 0.01) suggested that the HM data were suitable for PCA [41]. To guarantee a sufficiently large total variance (≥85%), four principal components (PCs) were extracted in which PC1 accounted for 34.63% of the total variance, exhibiting a high loading (≥50%) of V, Cr, Co, Cu, Zn, and Pb. Considering the low enrichment levels of V, Cr, Co, and Pb, PC1 may be attributed to soil parent material; the total variance of PC2 was 25.67%, which was dominated by As (87%), and Cu (76%); PC3 and PC4 accounted for 18.52% and 11.45% of the total variance dominated by Hg (96%) and Cd (99%), respectively. Due to the high enrichment level of Hg, Cd, and As, all of the PC2, PC3, and PC4 were associated with different anthropogenic sources [4,21].

3.4.2. Source Apportionment by PMF Model

The number of factors was set to four based on the above multivariate statistical analysis. The most suitable results were obtained after 150 operations. The fitting coefficients (R²) of V, Cr, Co, Cu, Zn, Cd, Pb, Hg, and As were 0.793, 0.979, 0.909, 0.866, 0.880, 0.997, 0.769, 0.985, and 0.714, respectively. The results of the PMF model are shown in Figure 8.

The first factor (F1) exhibited a total contribution of 25.4% and was dominated by Cu (39.3%), As (35.2%), V (30.1%), Co (26.2%), and Zn (21.7%). As mentioned above, Zn and As had been enriched in the soil due to anthropogenic activities. Previous studies had indicated that Zn and As were commonly used as additives to livestock feed, and they entered into agricultural soils via application of livestock manure as fertilizers [42,43]. This was further verified by the fact that the concentration of Zn also increased along with cultivation year (Figure 4). Thus, F1 was assigned to livestock manure.

The second factor (F2) with a total contribution of 35.5%, was dominated by Cr (53.4%) and contained moderate loadings of V (31.3%), Co (30.1%), Cu (30.4%), Zn (36.4%), Pb (25.6%), and As (25.6%). The mean contents of V, Co, Cu, Zn, and Pb were below the background value, and the CV values for V, Cr, Co, Cu, and Pb were below 0.3, indicating that they were uniformly distributed with less anthropogenic perturbation. Furthermore, V, Cr, and Co were usually representative elements of natural sources, and they were widely presented in the pedogenic process and soil parent material [20,44,45]. Therefore, F2 stood for the lithogenic origin.

The third factor (F3) exhibited a total contribution of 21.5% and was dominated by Hg (64.1%). According to the above analysis, Hg was the most serious pollutant, exhibiting the highest ecological risk, mainly associated with anthropogenic sources. In the study area, long heating period and industrial activities were accompanied by a large amount of coal combustion. Previous studies had shown that the Hg content in Chinese coal was 0.17 mg/kg of which a relatively large amount entered the atmosphere with coal combustion and eventually entered the agricultural soil [38,46]. Therefore, F3 represented coal combustion.

The fourth factor (F4), with a total contribution of 17.7%, was dominated by Cd (64.1%). The high CV value and pollution levels indicated that Cd was influenced by anthropogenic activities. Cd was an indicator element for chemical fertilizers, which was positively correlated with the application of chemical fertilizers in China [47,48,49]. The result also indicated that the concentration of Cd also increased along with cultivation year (Figure 4). Thus, F4 represented the application of chemical fertilizers.

3.5. Ecological Risk Assessment of Pollution Sources

The mean RI values of different pollution sources (Figure 9) were as follows: coal combustion (196.28) > chemical fertilizers (74.02) > livestock manure (71.04) > lithogenic origin (70.24), with the proportions of 47.69%, 17.98%, 17.26%, and 17.07%, respectively. Coal combustion was a prior pollution source for RI; consequently, more attention should be paid to coal combustion in the prevention and control of ecological risk for soil HMs. Figure 9 showed that the spatial distribution of ecological risk to soil HMs was different among pollution sources. For coal combustion, the southeastern portion of the study area exhibited the highest ecological risk, which was attributed to the presence of the Guoerluosi Industrial Park in this region. For livestock manure and fertilizer application, the southern (Hongguang Farm) and southeastern portions (Hongqi Farm) of the study area showed the highest ecological risks, which were attributed to their long cultivation periods. Therefore, to mitigate the ecological risk of soil HMs, monitoring and management of coal combustion and agricultural fertilizers need to be strengthened in the southern and southeastern parts of the area.

To summarize, the prior pollutants for HMs in Qianguo Irrigation District were Hg and Cd, and the prior pollution source was coal combustion. The coal combustion plant is essential for electric power generation and heating in Northeast China, whereas China has already committed to peak carbon dioxide emissions before 2030 and achieve carbon neutrality before 2060. Thus, it can be expected that the consumption of coal will reduce in the near future accompanied by the reduction of HM emission. According to our investigation, clean energy, such as wind-generated electricity and solar power stations, has been applied in the research area, which may facilitate the reduction of not only caron emissions but also HM pollution. In addition, livestock manure and chemical fertilizers were the main anthropogenic sources. Currently, the direct application of livestock manure as fertilizers has been banned in China since 2020; thus, it may not contribute largely to HM pollution for agricultural soil in the future. Chemical fertilizers were the most intractable pollution source. It is difficult for farmers to reduce the application amount spontaneously, even though certain slow-release fertilizers have been developed to improve the utilization efficiency, which may cut down the contribution to soil HM pollution. Therefore, it is essential for the government to enhance the management and monitoring of agricultural fertilizer applications.

4. Conclusions

In conclusion, the environmental risk of soil HMs were comprehensively assessed in Qianguo Irrigation District, Northeast China, and the sources and source-oriented ecological risk were quantitatively identified, and the main conclusions are as follows: the mean concentrations of Zn, Cd, Hg, and As exceeded the local background values by factors of 1.40, 3.34, 6.99, and 1.04, respectively. According to the results of CF, PLI, I_geo, EF, INI, and RI, the HM pollution in the study area was classified into a moderate contaminated level with a considerable ecological risk, and Hg and Cd were the main contributors to it. In addition, the non-carcinogenic risks of soil HMs were acceptable (HI < 1) for all the populations, whereas the carcinogenic risks were non-negligible (10⁻⁶ < TCR < 10⁻⁴). Based on the combination of the PMF model and multivariate statistical analyses, the sources of soil HMs, including lithogenic origin (35.5%), livestock manure (25.4%), coal combustion (21.5%), and chemical fertilizers (17.7%), were quantitatively identified. Livestock manure was the main anthropogenic source in the study area. The ecological risk assessment of pollution sources for soil HMs indicated that coal combustion accounted for 47.69% of the RI contribution and was identified as the prior pollution source. To reduce the enrichment and ecological risk of HMs in soil, it is necessary to seek green energy substitutes. This study provided a scientific basis of HM pollution control in agricultural soils in Northeast China.