Pollution Characteristics, Spatial Patterns, and Sources of Toxic Elements in Soils from a Typical Industrial City of Eastern China

Abstract

Soil pollution due to toxic elements (TEs) has been a core environmental concern globally, particularly in areas with developed industries. In this study, we sampled 300 surface (0–0.2 m) soil samples from Yuyao City in eastern China. Initially, the geo-accumulation index, potential ecological risk index, single pollution index, and Nemerow composite pollution index were used to evaluate the soil contamination status in Yuyao City. Ordinary kriging was then deployed to map the distribution of the soil TEs. Subsequently, indicator kriging was utilized to identify regions with high risk of TE pollution. Finally, the positive matrix factorization model was used to apportion the sources of the different TEs. Our results indicated that the mean content of different TEs kept the order: Zn > Cr > Pb > Cu > Ni > As > Hg ≈ Cd. Soil pollution was mainly caused by Cd and Hg in the soil of Yuyao City, while the content of other TEs was maintained at a safe level. Regions with high TE content and high pollution risk of TEs are mainly located in the central part of Yuyao City. Four sources of soil TEs were apportioned in Yuyao City. The Pb, Hg, and Zn contents in soil were mainly derived from traffic activities, coal combustion, and smelting. Meanwhile, Cu was mainly sourced from industrial emissions and atmospheric deposition, Cr and Ni mainly originated from soil parental materials, and Cd and As were produced by industrial and agricultural activities. Our study provides important implications for improving the soil environment and contributes to the development of efficient strategies for TE pollution control and remediation.

1. Introduction

Soil pollution by toxic metal (TEs) has attracted wide attention both domestically and internationally, with public concern mainly regarding the safety of agricultural products and human health [1,2,3,4,5,6]. Among the TEs, chromium (Cr), lead (Pb), cadmium (Cd), mercury (Hg), arsenic (As), copper (Cu), zinc (Zn), and nickel (Ni) were listed as priority monitoring pollutants regulated by the Chinese government. Owing to rapid advances in industrialization and urbanization during the last several decades, TE pollution in soil has become increasingly serious and widespread in China [7,8,9,10,11].

TEs are important natural elements of the Earth’s crust. Before the emergence of human civilization, soil TEs were mainly derived from natural sources, and the content of TEs in soil remained at a low level. However, after the emergence of human civilization, especially since the industrial revolution in the 1840s, anthropogenic inputs of TEs to soils remarkably exceeded natural inputs from pedogenesis [12]. Numerous studies have confirmed that anthropogenic activities, such as industrial activities, agrochemical product input, sewage irrigation, smelting operations, fossil fuel combustion, and atmospheric deposition are important sources of TEs in the soil [13,14,15,16,17,18,19,20].

Assessing the pollution status of soil TEs provides vital knowledge for understanding the soil quality and the potential health risk posed to humans due to exposure to soil TEs. Currently, the geo-accumulation index (I_geo), potential ecological risk index (ER and RI), single pollution index (PI), and Nemerow composite pollution index (NPI) are the most widely used to evaluate soil pollution caused by TEs [21,22,23,24,25,26]. Moreover, analysis of the source of soil TEs could serve as a basis for accurately reducing the input of soil sourced from anthropogenic activities and can contribute to protecting human health and the soil environment. As a multivariate receptor factor analysis model, the PMF has proven to be a feasible and useful method for apportioning sources of soil TE [27,28].

In addition, mapping soil TEs provides critical information for understanding the spatial patterns of soil TEs. Many spatial interpolation methods, such as the triangular interpolation network (TIN), inverse distance weighting (IDW), regression models (RM), thin plate splines (TPS), local trend surfaces (LTS), geographically weighted regression (GWR), and geostatistical methods have been extensively used to map soil properties [29]. Among these methods, ordinary kriging (OK) has been the most widely used for mapping the spatial pattern of soil properties characterized by optimal unbiased estimation [29]. Moreover, indicator kriging (IK) can use indicator (0 or 1) variables to generate probabilities that a critical value was exceeded or not at each location in the study area, which can help us identify areas with a high risk of soil TE pollution [30,31,32]. Identifying regions with high risk of soil TE pollution could save considerable time, as well as labor costs, improving the efficiency and accuracy of soil pollution remediation.

Yuyao City is one of the most developed cities and counties in China. It is well-known for its electronic appliances, mechanical hardware, plastics, and molds. In addition, Yuyao City has more than seven thousand years of rice planting history, and it is an important production area for food and agricultural products. Long-term and extensive industrial and agricultural actions result in the accumulation of TEs in the soil and pose a great health threat to residents [33,34,35,36,37]. For example, excessive intake of Pb can lead to adverse effects on the immune, nervous, and endocrine systems [3,4,38,39,40]. Bladder, kidney, and liver cancers can be induced by long-term exposure to As [3,4,41,42]. Massive exposure to high doses of Zn, Cr, and Cu can lead to impairments in cholesterol, fertility, and liver function [43,44]. However, until now, the knowledge related to the pollution status, spatial pattern, pollution risk, and related sources of soil TEs in Yuyao City is still scarce. Therefore, to fill these gaps, in this study, we analyzed the pollution status, spatial patterns, and origins of TEs in Yuyao City based on 300 surface soil samples from Yuyao City. The main aims of this study are: (1) summarizing the statistics of soil TEs and assessing the pollution status of soil in Yuyao City using the I_geo, ER, RI, PI, and NPI index; (2) analyzing the spatial pattern of soil TEs using the OK method and identifying regions with high risk of soil TE pollution in Yuyao City using the IK method; and (3) apportioning the sources of different TEs via the PMF model. Our study could provide useful implications for environmental protection in the survey region and can also provide a valuable foundation for similar studies in other regions.

2. Study Area and Materials

2.1. Study Area

Yuyao City is situated in the northeast of Zhejiang Province, eastern China (Figure 1) with area of around 1500 km². The Yuyao City climate has a subtropical monsoon climate with an average annual temperature of 18.1 ℃ and annual precipitation of 1835.2 mm in 2020. Yuyao City is one of the most developed cities or counties in China, featuring developed hardware products, home appliance manufacturing, and plastic product industries.

2.2. Soil Sampling and Analysis

In this study, the process of soil sampling and analysis strictly refer to the national standard issued by the Ministry and Ecology and Environment of People’s Republic of China (The Technical Specification for soil Environmental monitoring, HJ/T 166-2004, https://www.mee.gov.cn/ywgz/fgbz/bz/bzwb/jcffbz/200412/t20041209_63367.shtml, accessed on 21 October 2021). We collected 300 surface soil samples (0–0.2 m) from Yuyao City in 2012 (Figure 1) based on spatial coverage, land use, and soil groups. In some areas which have been known as heavily polluted regions or an important production base for rice and others agricultural crops, we conducted a more extensive survey, and thus more samples were collected from these regions. The coordinates of the sampling sites were recorded with a portable GPS. The soil samples were then transported to the laboratory and air-dried for several days at ambient temperature. After that, all soil samples were sieved with a 2 mm nylon sieve, meanwhile, plant residuals and gravels were removed. The concentrations of soil Cr, Pb, Cd, Cu, Zn, and Ni were measured by inductively coupled plasma mass spectrometry (ICP-MS). The contents of Hg and As in the soil were determined by atomic fluorescence spectrometry (AFS). Reagent blanks and standard reference materials (GBW07405, supplied by IGGE at the Chinese Academy of Geological Sciences) were used for quality assurance and quality control with recovery rates varying between 90% and 110%.

2.3. Geo-Accumulation Index (I_geo)

The $I_{geo}$ is calculated to evaluate the enrichment status of soil TEs by comparing soil background values. $I_{geo}$ is determined by the following formula [45]:

$$I_{geo}={\log }_2(\frac{C_i}{1.5*C_0})$$

(1)

here, $C_i$ means the content of soil TE $i$, $C_0$ represents the corresponding soil background value [46]. The $I_{geo}$ is classified into seven classes: $I_{geo}$ < 0 (unpolluted), 0 < $I_{geo}$ ≤ 1 (unpolluted to moderately polluted), 1 < $I_{geo}$ ≤ 2 (moderately polluted), 2 < $I_{geo}$ ≤ 3 (moderately to heavily polluted), 3 < $I_{geo}$ ≤ 4 (heavily polluted), 4 < $I_{geo}$ ≤ 5 (heavily to extremely polluted), and $I_{geo}$ > 5 (extremely polluted).

2.4. Potential Ecological Risk Assessment (ER and RI)

The ER and RI combine the concentration of each TE with its toxicity efficiency [47]. The ER is calculated using the following formula below [48]:

$$RI={\displaystyle \sum _{}^{}E{R}^i;E{R}^i=T{R}^i*C_f^i};C_f^i=\frac{C_0^i}{C_n^i}$$

(2)

where $C_f^i$ indicates the contamination factor for TE i, $C_0^i$ denotes the content of soil TE i, $C_n^i$ on behalf of the background value of TE i issued by the China National Environmental Monitoring Center [46], ERⁱ means the single potential ecological risk of TE i l, and TRⁱ is the toxic coefficient of TE i [47]. The toxic coefficients of Cr, Pb, Cd, Hg, As, Cu, Zn, and Ni were 2, 5, 30, 40, 10, 5, 1, and 3, respectively. RI is the composite potential ecological risk index, which is an integrated value of several TEs. ER is classified into five grades [48]: ER < 40 (low risk), 40 < ER < 80 (moderate risk), 80 < ER < 160 (considerable risk), 160 < ER < 320 (high risk), and ER > 320 (very high risk). The RI is classified into four categories [48]: RI < 150 (low risk), 150 < RI < 300 (moderate risk), 300 < RI < 600 (considerable risk), and RI > 600 (high risk).

2.5. Single Pollution Index (PI) and Nemerow Composite Pollution Index (NPI)

The PI and NPI were also computed to evaluate the pollution status of TEs by comparing the content of soil TEs with national screening values regulated by the Chinese government [49]. The PI was calculated using the formula:

$$PI=\frac{C_i}{S_i}$$

(3)

here, $C_i$ on behalf of the concentration of soil TE i, and $S_i$ means the national screening value of TE $i$ [49].

$$NCPI=\sqrt{\frac{{(P_{i\max })}^2+{(\overline{PI})}^2}{2}}$$

(4)

where $P_{i\max }$ represents the maximum PI of TEs and $\overline{PI}$ means the averaged PI of TEs [50]. PI could be classified into five categories: safety (PI ≤ 1), slight pollution (1 < PI ≤ 2), mild pollution (2 < PI ≤ 3), moderate pollution (3 < PI ≤ 5), and severe pollution (PI > 5). The NPI value is classified as: safe (NPI ≤ 0.7), precaution (0.7 < NPI ≤ 1), slightly polluted (1 < NPI ≤ 2), moderately polluted (2 < NPI ≤ 3), and seriously polluted (NPI > 3).

2.6. Ordinary Kriging

The spatial variability of soil TEs is a core concern for soil TE pollution control and remediation. As one of the most classic interpolation algorithms, the OK method has been widely used to estimate the spatial variation of soil properties around the world [29]. In this study, the OK was deployed to map the content of TEs in the soil around the survey region. In OK, the experimental semi-variogram was fitted to reveal the spatial dependence of soil TEs using formula as follows [51]:

$$Z*(X_0)={\displaystyle \sum _{i=1}^n{\phi }_{i}Z(x_i)}$$

(5)

where $Z(x_i)$ is the observed TE content at sampling site i, and n is the number of observed soil samples. φ_i denotes the optimal weight value that could contribute to the unbiased estimation of the TE content with minimum variance. Detailed information on OK is available from Webster and Oliver (2007) [52].

2.7. Indicator Kriging

Indicator kriging (IK) is an extension of OK that can estimate the probability that the target variable exceeds a critical threshold at unvisited points [32,53,54]. The IK was deployed on a binary-transformed sample dataset. First, the indicator is defined using the following formula [52]:

$$I(x;z)=\{\begin{array}{cc}1& Z(x)\\ 0& Z(x)\end{array}\begin{array}{cc}\ge & z\\ \le & z\end{array}$$

(6)

Subsequently, using Equation (6), the original variable is transformed into a new set of binary indicator variable. Then, the experimental semi-variogram is as follows:

$${\gamma }_1^{*}(h)=\frac{1}{2N(h)}{\displaystyle \sum _{i=1}^{N(h)}{[I(x_i;Z_k)-I(x_i+h;Z_k)]}^2}$$

(7)

where ${\gamma }_1^{*}(h)$ denotes the experimental semi-variogram of the newly transformed binary indicator. $N(h)$ means the number of pairs of indicator transformations, $I({\chi }_i;Z_k)$ and $I({\chi }_i+h;Z_k)$, with a distance lag of h.

The conditional cumulative distribution function is produced by the IK at unvisited sites:

$$F({\chi }_0;Z_k|(n))=I^{*}(x_0;Z_k)={\displaystyle \sum _{i=1}^n{\lambda }_{i}I({\chi }_i;Z_k)}$$

(8)

here, $I^{*}({\chi }_0;Z_k)$ on behalf of the predicted indicator at the unvisited sites ${\chi }_0$, and ${\lambda }_1$ represents the weight assigned to the indicator in sites $x_i$.

2.8. US-EPA Positive Matrix Factorization (PMF)

The PMF (EPA PMF version 5.0) model is a multivariate factor analysis method [55], which was used in current research to identify the potential sources of soil TEs in Yuyao City. In PMF, the soil sample data are decomposed into two matrices: factor contributions and factor profiles. After that, the users could appoint the sources of TEs that may contribute to the samples based on source profile information and emissions inventories. A detailed description of the PMF method is available in the study organized by Bhuiyan et al. (2021), Wang et al. (2021), Wu et al. (2021) [56,57,58].

2.9. Data Analysis

In this study, descriptive statistical analysis and correlation analysis were finished by R software (R Development Core Team, 2017). OK and IK were carried out with ArcGIS 10.3, (ESRI ArcGIS 10.3, Redlands, CA, USA). The PMF model was conducted using EPA PMF software (version 5.0).

3. Results and Discussion

3.1. Summary Statistics of Soil TEs

The descriptive statistics of the TEs in Yuyao City are listed in

Table 1. Summary statistics of TEs in soil of Yuyao City (mg/kg) (n = 300).

TEs	Min	Max	Mean	SD	Skewness	Kurtosis	CV(%)	BV	TC	NV
Cr	9.21	246.00	61.80	22.12	2.05	16.26	35.79	52.90	2	61
Pb	17.30	66.40	34.43	1.58	0.58	−0.49	4.59	23.70	5	27
Cd	0.05	1.69	0.20	0.11	8.74	114.43	55.00	0.07	30	0.097
Hg	0.04	0.76	0.20	0.14	1.15	1.13	70.00	0.09	40	0.065
As	1.90	11.40	6.13	1.80	0.26	0.30	29.36	9.2	10	11
Cu	7.14	134.00	32.75	16.20	2.99	12.06	49.47	17.60	5	23
Zn	47.70	252.00	97.61	27.89	2.13	6.76	28.57	70.60	1	74
Ni	4.86	144.00	26.52	11.08	3.96	40.94	41.78	24.60	3	27

Note: CV means coefficient of variation, %; BV means background value of TEs in Zhejiang Province; TC means the toxic coefficient of different TEs [47]; NV means the national arithmetic value reported by Teng et al. (2014) [59].

. The averaged concentration of different TEs kept the order: Zn (97.61 mg/kg) > Cr (61.80 mg/kg) > Pb (34.43 mg/kg) > Cu (32.75 mg/kg) > Ni (26.52 mg/kg) > As (6.13 mg/kg) > Hg (0.20 mg/kg) ≈ Cd (0.20 mg/kg). The averaged concentration of all TEs were greater than corresponding soil background values in Zhejiang Province, in addition to As. In comparison with the national arithmetic averaged concentrations of different TEs in China reported by Teng et al. (2014) [59], the levels of Pb, Cd, Hg, Cu, and Zn in Yuyao City were clearly higher than their national values. Cr content was slightly higher than the corresponding national value (Table 1). Moreover, the contents of As and Ni were lower than their national values (Table 1).

The coefficient of variation (CV, %) was used to measure variation in the target variable [60,61]. The CV of different TEs followed a decreasing order: Hg (70.00%) > Cd (55.00%) > Cu (49.47%) > Ni (41.78%) > Cr (35.79%) > As (29.36%) > Zn (28.57%) > Pb (4.59%), indicating that Hg and Cd exhibited strong variation and suggested significant spatial heterogeneity of Hg and Cd in Yuyao City [61]. Our results are consistent with the results obtained in many other similar studies, such as those by Xia et al. (2020), Ren et al. (2021), and Wang et al. (2021) [62,63,64].

3.2. Correlation Analysis

In this study, we used the Spearman’s correlation method to analyze the associations among the different TEs. As presented in Figure 2, most TEs under study showed a significant positive correlation with each other in Yuyao City. In particular, Cr had a strong relationship with Ni (0.94), Pb had a strong relationship with Hg (0.70), and Cu had a strong relationship with Zn (0.63). TEs are highly related to each other and may share similar sources which could contribute to the identification of sources of soil TEs and more efficient control of soil TE contamination [47].

3.3. Pollution and Ecological Risk Assessment

3.3.1. The Geo-Accumulation Index

The I_geo was calculated to quantitatively evaluate the accumulation of soils in the study area. A boxplot of I_geo is shown in Figure 3. The average I_geo values for different TEs followed the order: Cd (0.84) > Hg (0.27) > Cu (0.19) > Pb (−0.13) > Zn (−0.17) > Cr (−0.46) > Ni (−0.59) > As (−1.24). The results revealed that evaluated by the geo-accumulation index, the pollution status for Pb, Zn, Cr, Ni, and As was at an unpolluted level, while the pollution status for Cd, Hg, and Cu was unpolluted to moderately polluted.

Table 2. Pollution status assessed by geo-accumulation index (n = 300).

Grade		Cr	Pb	Cd	Hg	As	Cu	Zn	Ni
Unpolluted $I_{geo}$ ≤ 0	n	266	158	210	209	300	278	297	299
	Percent (%)	88.67	52.67	70	69.67	100	92.67	99	99.67
Unpolluted to moderately polluted 0 < $I_{geo}$ ≤ 1	n	33	142	0	0	0	0	0	0
	Percent (%)	11	47.33	0	0	0	0	0	0
Moderately polluted 1 < $I_{geo}$ ≤ 2	n	1	0	85	84	0	19	3	1
	Percent (%)	0.33	0	28.33	28	0	6.33	1	0.33
Moderately to heavily polluted 2 < $I_{geo}$ ≤ 3	n	0	0	4	7	0	3	0	0
	Percent (%)	0	0	1.33	2.33	0	1	0	0
Heavily polluted 3 < $I_{geo}$ ≤ 4	n	0	0	0	0	0	0	0	0
	Percent (%)	0	0	0	0	0	0	0	0
Heavily to extremely polluted 4 < $I_{geo}$ ≤ 5	n	0	0	1	0	0	0	0	0
	Percent (%)	0	0	0.33	0	0	0	0	0
Extremely polluted $I_{geo}$ > 5	n	0	0	0	0	0	0	0	0
	Percent (%)	0	0	0	0	0	0	0	0

shows the proportions of different geological accumulation grades of different TEs in the whole survey region, referring to the I_geo classification standard [35]. Obviously, the soil pollution caused by Cd and Hg is much more serious than that caused by Cr, Pb, As, Cu, Zn, and Ni. Regarding As, all soil samples were unpolluted, while most of the soil samples were at unpolluted levels for Cr, Cu, Zn, and Ni (Table 2). In terms of Pb, 47.33% of the soil samples belonged to the unpolluted to moderately polluted grade, while the remaining soil samples were unpolluted. As for Cd, 28.33%, 1.33%, and 0.33% of the soil samples were moderately polluted, moderately to heavily polluted, and heavily to extremely polluted, respectively. With regard to Hg, 28.00% and 2.33% of soil samples were moderately polluted and moderately to heavily polluted, respectively.

3.3.2. The Potential Ecological Risk Assessment

The ER values of different TEs and RI were used to quantitatively evaluate the potential ecological risk due to TE accumulation. The mean value of ER for different TEs followed the order: Hg (89.57) > Cd (86.32) > Cu (9.30) > Pb (7.26) > As (6.67) > Ni (3.23) > Cr (2.34) > Zn (1.38) (Figure 4). The results showed that soil in Yuyao City had considerable potential ecological risk results from Hg and Cd accumulation, while the potential ecological risk caused by other TEs was low (Figure 4). The mean RI value was 206.08, which is at the moderate risk grade.

The proportions of soil samples at different potential ecological risk levels for different TEs are listed in

Table 3. Pollution status assessed by potential ecological risk assessment index (n = 300).

Grade		Cr	Pb	Cd	Hg	As	Cu	Zn	Ni
Low risk ER ≤ 40	n	300	300	9	73	300	300	300	300
	Percent (%)	100	100	3	24.33	100	100	100	100
Moderate risk 40 < ER ≤ 80	n	0	0	147	89	0	0	0	0
	Percent (%)	0	0	49	29.67	0	0	0	0
Considerable risk 80 < ER ≤ 160	n	0	0	137	103	0	0	0	0
	Percent (%)	0	0	45.67	34.33	0	0	0	0
High risk 160 < ER ≤ 320	n	0	0	6	34	0	0	0	0
	Percent (%)	0	0	2	11.33	0	0	0	0
Very high risk	n	0	0	1	1	0	0	0	0
ER > 320	Percent (%)	0	0	0.33	0.33	0	0	0	0

. The potential ecological risk of all soil samples was low for Cr, Pb, As, Cu, Zn, and Ni. With regard of Cd, the proportions of ER value of soil samples at low, moderate, considerable, high, and very high risk were 3%, 49%, 45.67%, 2%, and 0.33%, respectively. For Hg, the proportions of ER value of soil samples at low, moderate, considerable, high, and very high risk were 24.33%, 29.67%, 34.33%, 11.33%, and 0.33%, respectively. For the RI value, the proportions of soil samples at low, moderate, considerable, and high risk were 35.67%, 49.33%, 14.67%, and 0.33%, respectively.

3.3.3. Pollution Index

Moreover, we calculated the PI and NPI to evaluate the contamination status of TEs in the soil of Yuyao City. As presented in Figure 5, the averaged PI decreased by the following order: Cd (0.48) > Zn (0.42) > As (0.34) ≈ Hg (0.33) ≈ Pb (0.33) > Ni (0.29) > Cr (0.23) ≈ Cu (0.23). The averaged PI of all TEs were at a safe level, but it is worth noting that the PI of some soil samples were at polluted levels for Cd, Hg, Cu, Zn, and Ni. Additionally, the average NPI value was 0.46, which was also at a safe level.

The proportions of PI values of TEs in soil samples at different pollution levels are provided in

Table 4. Pollution status assessed by pollution index (n = 300).

Grade		Cr	Pb	Cd	Hg	As	Cu	Zn	Ni
Safety PI ≤ 1	n	300	300	284	291	300	298	296	299
	Percent (%)	100	100	94.67	97	100	99.33	98.67	99.67
Slight pollution 1 < PI ≤ 2	n	0	0	15	9	0	2	4	0
	Percent (%)	0	0	5	3	0	0.67	1.33	0
Mild pollution 2 < PI ≤ 3	n	0	0	0	0	0	0	0	1
	Percent (%)	0	0	0	0	0	0	0	0.33
Moderate pollution 3 < PI ≤ 5	n	0	0	1	0	0	0	0	0
	Percent (%)	0	0	0.33	0	0	0	0	0
Severe pollution PI > 5	n	0	0	0	0	0	0	0	0
	Percent (%)	0	0	0	0	0	0	0	0

. Similar to the I_geo and ER indices, the PI of Cr, Pb, and As in all soil samples were at a safe level, while only a small part of the soil samples had PI values at polluted levels (slight pollution or more serious) for Cd (5.33%), Hg (3%), Cu (0.67%), Zn (1.33%), and Ni (0.33%). With regard to the NPI value, the proportion of soil samples with NPI values at safe, precautionary, slightly polluted, moderately polluted, and seriously polluted levels were 86.33%, 11%, 2.33%, 0%, and 0.33%, respectively.

Overall, Cd and Hg clearly accumulated in the soil in Yuyao City. Generally, the pollution was mainly caused by Cd and Hg, while the content of other TEs was kept at a safe level. This result is consistent with previous surveys conducted in Yuyao City [65,66,67]. Yuyao City is ranked 12th among the top 100 counties in China according to the level of economic development. Extensive industrial activities may lead to the accumulation of soil Cd and Hg [68,69].

3.4. Identify Regions with High Risk of TEs Contamination

Traditional spatial interpolation algorithms, such as OK and IDW, can only map the specific value of the target variable in unobserved locations. However, the predictions made by all interpolation methods are biased, which means that uncertainty exists in the predicted results. To fill this gap, IK was developed to map the probability that the target variable exceeds a specific threshold. In this study, we used IK to predict the probability of the occurrence of soil TE pollution in Yuyao City. For the geo-accumulation status, we focused on Cd and Hg because Cd and Hg are the most seriously polluted TEs in Yuyao City. Meanwhile, for the potential ecological risk and pollution index, we focused on NPI and RI because it can represent the composite pollution status and ecological risk caused by different TEs in the soil.

As shown in Figure 6a, the high pollution risk caused by Cd accumulation (probability of I_ego > 1) was mainly expected in the central part of Yuyao City. In some regions, the probability of Cd pollution could exceed 90%. In terms of Hg, regions with a high pollution risk were situated in the center of Yuyao City, but most of the regions had a pollution risk lower than 50%. For the NPI value, areas with high composite pollution risk were mainly situated in the central part of Yuyao City, while the high value of probability of potential ecological risk (RI > 150) occurred mainly in the central and middle parts of Yuyao City.

In summary, the central part of Yuyao City had a relatively higher risk of soil TE contamination, especially for Cd and Hg. This indicates that these regions should be listed as priority regions for soil TE contamination control and remediation. There is an urgent need to take some measures to protect residents from the health risks caused by soil TE contamination. For example, the crop planting system could be adjusted in regions with high pollution risk, and economic plants could be planted take the place of food crops to reduce the accumulation of TEs in the human body via the food chain. Planting crops with low accumulation of TEs is another efficient way to reduce the risk of TE accumulation.

3.5. Spatial Pattern of TEs

Obtaining accurate information on the spatial variability of soil TEs plays a critical role in controlling TE pollution. In current research, OK interpolation was conducted to estimate the spatial variation of TE concentration in the soil of Yuyao City (Figure 7a–h). As presented in Figure 7, the concentrations of Cr, Cd, Hg, As, Zn, and Ni showed similar spatial patterns, with high values detected in the central and southern parts of Yuyao. Pb and Cu also had similar spatial patterns, with high values majorly detected in the central part of Yuyao City. According to the map of land use of Yuyao City presented in Figure 7i, the land use types in the central part of Yuyao are mainly urban, rural, industrial, mining, and residential land. Anthropogenic activities, such as industrial activities (emission of solid waste and wastewater), traffic activities (vehicle emissions), and agricultural activities (application of chemical fertilizer and pesticide) discharge large amounts of waste in which large amounts of TEs exist and then accumulate in soil [8,70,71,72,73,74].

3.6. Source Apportionment by PMF

PMF has been commonly used to apportion the source of soil TEs [57,75,76,77,78]. In current research, the PMF method was deployed to compute the contribution of different sources to the accumulation of eight soil TEs in Yuyao City. As displayed in Figure 8 and Figure 9, four sources were identified by the PMF model.

Factor 1 was heavily loaded with Pb (59.16%), Hg (75.96%), and Zn (35.05%). As presented in Figure 2, Pb was significantly and highly correlated with Hg and Zn, which also verified that Pb, Hg, and Zn were derived from similar sources in the soil of Yuyao City. Numerous studies have verified that vehicle exhaust emissions are a critical source of Pb accumulation in the soil [79,80]. As listed in Table 1, Hg had the highest CV value, which means that anthropogenic activities are the major origin of soil Hg [21]. Industrial activities such as smelting, coal mining, and combustion are regarded as the main sources of Hg pollution in the soil [21] In addition, coal mining and combustion emissions are important sources of Pb in the soil. Moreover, previous studies have also reported that Zn is closely associated with anthropogenic actions, such as burning fossil fuels, smelting, and coal combustion emissions [81,82]. The high values of Pb, Hg, and Zn were mainly situated in the central part of Yuyao City, which overlays with the urban land. Thus, Factor 1 could represent the source of traffic activities, coal combustion, and smelting (Figure 6).

Factor 2 was dominated by Cu, with a very high loading of 71.69%. Generally, Cu is closely related to the emission of industrial activities like metal smelting, petrochemical, cement production, and so on [83,84]. In addition, some researchers have also found that atmospheric deposition is another significant contributor to Cu accumulation in the soil [57,85]. Therefore, Factor 2 was mainly derived from industrial emissions and atmospheric deposition.

Factor 3 was dominated by Cr and Ni with high factor loading values of 74.28% and 76.06%, respectively. As presented in Figure 2, the content of Cr was extremely high, related to the content of Ni in the soil of Yuyao City, with a correlation coefficient of 0.94. This demonstrates that Cr and Ni in the soil originated from the same source. The mean content of Cr and Ni in Yuyao City was slightly higher than the corresponding background value, revealing that the I_geo and ER values of Cr and Ni were also at the safe or low risk level, which means that the Cr and Ni in the soil of Yuyao City may result from the weathering of soil parental materials. Previously, many researchers, such as Salonen and Korkka-Niemi (2007) [86], Chen et al. (2016) [87], Wang et al. (2020) [88], and Bhuiyan et al. (2021) [57] also reported that the contents of Cr and Ni were mainly sourced from soil parental materials. Therefore, Factor 3 represented natural sources, such as soil parent materials.

Factor 4 was mainly loaded with As and Cd with loadings of 50.84% and 40.56%, respectively. The enrichment of As and Cd has been shown to be associated with agricultural activities, like applying chemical fertilizers, herbicides, and pesticides [89]. In addition, the utilization of wastewater for irrigation could also give rise to the accumulation of As in the soil [89,90]. Some studies also reported that As could probably be sourced from poultry farming because organic As is widely used as a feed additive for livestock to promote animal growth, which could then result in a high content of As in livestock manure [91,92]. In terms of Cd, numerous studies have proved that Cd originates from industrial activities such as smelter plant emission, discharge of industrial waste, and so on. Cd is also usually used as an input for industrial production of various products such as plastics [93,94]. The plastic industry is a pillar industry in Yuyao. Thus, the plastic industry may also contribute to Cd accumulation in the soil of Yuyao City. Therefore, Factor 4 represents the source related to industrial and agricultural activities.

In summary, the main sources of TEs in the soil of Yuyao City were industrial activities, industrial emissions and atmospheric deposition, natural sources, and agricultural activities. The results of source apportionment of soil TEs were also verified by the results of summary statistics listed in Table 2, correlation analysis (Figure 2), and the map of the content of soil TEs in Yuyao City produced by OK (Figure 6) interpolation and the map of pollution risk outputted by IK (Figure 7). This confirmed the feasibility of using PMF for quantitatively analyzing the sources of soil TEs. However, as reported by previous researchers, uncertainty still exists in the results produced by the PMF model, which highlights the importance of introducing and combined with other source apportionment methods such as cluster analysis [12], principal component analysis [90], machine learning methods [17], the Unmix model [95], the stable isotope ratio technique [96], and the finite mixture distribution model [97] to reduce the uncertainty of the results produced by PMF.

3.7. Implication and Policy Recommendation

The results we obtained in from presented research could provide important implication for decision-makers to make effective management measures to control soil TEs contamination and thus improve environmental quality. Firstly, regions with high pollution risk should be listed as priority areas for soil TEs contamination control and remediation. Some management practices, such as lime application and bioleaching, could be employed in acidic soils to regulate soil pH and reduce the availability of HMs [98,99,100,101]. Secondly, we could consider establishing a spatial buffer around industrial park or regions featured by high density of industrial companies. We should also keep a safe distance between polluted land or industrial park and residents. Thirdly, breeding crop cultivars with low accumulation is another promising way to mitigate soil TEs contamination in food crops. Crops with different cultivars showed strong variation in ability of absorbing TEs from soil [7]. Finally, according to the results presented in Section 3.6, burning fossil fuels, traffic emissions, and coal combustion emissions contributed to Pb, Hg and As accumulation in the study area. Therefore, the government could encourage adopting green travel methods like bike or public traffic as well as promote use of clean energy such as wind energy, hydropower, tidal energy, and solar energy. These measures are expected to reduce accumulation of TEs in soil through atmospheric deposition, emission of wastewater, and solid waste. In our further work, comparison with previous similar studies which made surveys in the same regions as well as conducting another soil sampling in the same sites is a promising way to detect the temporal change of TEs in soil.

4. Conclusions

Obvious accumulation of Cr, Pb, Cd, Hg, Cu, Zn, and Ni was detected in the soil of Yuyao City. The concentrations of Cr, Pb, As, Cu, Zn, and Ni in most of the soil samples were at safe levels, whereas clear pollution caused by Cd and Hg accumulation in soil was detected. The high contents of the different TEs were mainly situated in the central and middle parts of Yuyao City. In addition, there was a low pollution risk caused by Cr, Pb, As, Cu, Zn, and Ni in Yuyao City, while a high risk of Cd and Hg contamination was detected. The regions with a high risk of TE pollution were mainly distributed in central Yuyao City. In terms of the sources of different TEs in the soil, Pb, Hg, and Zn were mainly derived from traffic activities, coal combustion, and smelting, Cu mainly sourced from industrial emissions and atmospheric deposition, Cr and Ni mainly came from soil parental materials, and Cd and As were attributed to industrial activities and agricultural activities. This study could provide critical information for improving soil environmental quality and contribute to the development of more efficient policies for reducing TE accumulation in soil.