An Integrated Approach for Source Apportionment and Health Risk Assessment of Heavy Metals in Subtropical Agricultural Soils, Eastern China

: Unreasonable human activities may cause the accumulation of heavy metals (HMs) in the agricultural soil, which will ultimately threaten the quality of soil environment, the safety of agricultural products, and human health. Therefore, the accumulation characteristics, potential sources, and health risks of HMs in agricultural soils in China’s subtropical regions were investigated. The mean Hg, Cu, Zn, Pb, and Cd concentrations of agricultural soil in Jinhua City have exceeded the corresponding background values of Zhejiang Province, while the mean concentrations of determined 8 HMs were less than their corresponding risk-screening values for soil contamination of agricultural land in China. The spatial distribution of As, Cr, Ni, Cu, and Pb were generally distributed in large patches, and Hg, Zn, and Cd were generally sporadically distributed. A positive deﬁnite matrix factor analysis (PMF) model had better performance than an absolute principal component–multiple linear regression (APCS-MLR) model in the identiﬁcation of major sources of soil HMs, as it revealed higher R 2 value (0.81–0.99) and lower prediction error ( − 0.93–0.25%). The noncarcinogenic risks (HI) of the 8 HMs to adults and children were within the acceptable range, while the carcinogenic risk (RI) of children has exceeded the safety threshold, which needs to be addressed by relevant departments. The PMF based human health risk assessment model indicated that industrial sources contributed the highest risk to HI (32.92% and 30.47%) and RI (60.74% and 61.5%) for adults and children, followed by agricultural sources (21.34%, 29.31% and 32.94% 33.19%). Therefore, integrated environmental management should be implemented to control and reduce the accumulation of soil HMs from agricultural and industrial sources. PMF model better than APCS-MLR model in the identiﬁcation of major sources of soil HMs as it revealed higher R 2 value (0.81–0.99) and lower prediction error ( − 0.93–0.25%). PMF model also had higher accuracy than APCS-MLR model in the estimation of source contribution, and apportioned ﬁve sources, namely, agricultural natural sources trafﬁc sources industrial sources (13.42%), and atmospheric deposition (13.27%), respectively. Both RI and HI of adults are below the safety threshold, while the RI of children has exceeded the safety threshold, indicating that local children are already facing the threat of RI, and needs to be concerned by relevant departments. Moreover, a PMF based human health risk assessment model indicated that industrial sources contributed the highest risk to HI (32.92% and 30.47%) and RI (60.74% and 61.50%) for adults and children, followed by agricultural sources (21.34%, 29.31% and 32.94%, 33.19%). Therefore, integrated environmental management should be implemented to control and reduce the accumulation of soil HMs from these two sources. This study provided an effective approach to quantify the priority pollution sources of concern for human health risk, which is of great signiﬁcance for pollution control and risk reduction under limited resources.


Introduction
Heavy metals (HMs) in soil are characterized by high toxicity, nondegradability, weak mobility, and strong bioaccumulation [1,2]. With rapid urbanization, HMs pollution in agricultural soil caused by anthropogenic activities is becoming increasingly prominent and has become a global issue [3][4][5][6]. The difference between soil parent material and soil acidity and alkalinity is an important factor affecting the natural source of HMs. The subtropical red soil region is an important agricultural production base, as well as the most serious soil acidification area in China [7]. Soil acidification can not only lead to the leaching of nutrients [8] but also greatly improve the mobility and bioavailability of HMs [9,10]. The accumulation of HMs in agricultural soil can not only degrade soil structure and function, and thus reduce the quality and quantity of crop yield, but also posing considerable risks to human health via ingestion, inhalation, and skin contact [11,12]. Generally, soil HMs mainly originate from natural sources and anthropogenic sources. The former is mainly derived from soil parent material and the soil formation process [13,14], while the latter are associated with industrial emissions and discharge, coal combustion, atmospheric deposition, vehicle exhaust, and unreasonable agricultural inputs [12,15,16]. Therefore, accurately identifying the potential sources of HMs in subtropical agricultural soil and determining the risks to human health is essential for local prevention and remediation of soil HMs pollution is important.
Receptor models are commonly used for identifying the sources of a specific pollutant and quantitatively estimating the contributions of each source at certain receptor locations [17,18]. Various receptor models (for instance, absolute principal component scores-multiple linear regression (APCS-MLR) and positive matrix factorization (PMF)) have been widely applied for source apportionment of soil HMs. Compared with the traditional principal component analysis, APCS-MLR introduces a sample point with zero content, and then the multiple linear regression of the source can identify specific emission sources and quantify each emission influences [18]. PMF model considers the uncertainty and non-negative constraints of the data, which not only identify the number of soil HMs pollution sources and the contribution rate of each pollution source but also quantitatively identify the contribution rate of each pollution source to each HMs [13,19,20]. Due to the complex sources of soil HMs, the accuracy of the source-analyzed receptor model is affected by both sample errors and modeling errors and the result of individual model is difficult to verify [21]. Although many studies have conducted in-depth analyses on the sources of soil HMs [11,22], few studies compared the effects of multiple receptor models for source apportionment of HMs in subtropical agricultural soils. Some HMs will carry human blood through breathing, skin contact, diet intake, etc., and enter various organs and tissues of the human body along with blood circulation. When the HMs in the human body accumulate to a certain level, it will cause respiratory diseases, myocardial infarction (and other diseases), and serious consequences such as various cancers [23]. Presently, the risks of HMs to human health (carcinogenic risk and noncarcinogenic risk) have been received with wide concern [24][25][26]. However, the constituents of HMs varied from different sources, and the toxicities of various HMs were different; therefore, the human health risk posed by the highest source of HM concentrations is not necessarily the greatest [4,25]. Nevertheless, it is limiting to quantitatively determine the contribution of pollution sources of HMs in agricultural soils to human health risks. Therefore, an integrated approach combining human health risk assessment model with PMF model was proposed to quantitative human health risks from different sources of HMs and determine priority control sources in this study.
The subtropical region of eastern China has a suitable climate and adequate water and heat conditions, and is considered one of the important food bases in the Yangtze River Delta. Over the past 20 years, unreasonable human activities have led to the accumulation of HMs in the agricultural soils, which have eventually grown to threaten the soil quality and safety of agricultural products [27][28][29][30][31]. Therefore, the identification of HM pollution sources and health risks in this region is important for the remediation of agricultural soils in subtropical regions. The objectives of this study were to: (1) characterize the concentration and spatial distribution of agricultural soil HMs in Jinhua City, subtropical China; (2) compare the performance of APCS-MLR and the PMF model in the identification of potential sources and quantitation of the contribution of each source to soil HMs; and (3) establish an integrated health risk model to quantitatively determine contribution of pollution sources of HMs to human health risks.

Study Area
Jinhua City is located in the middle of Zhejiang Province (28 • 32 -29 • 41 N, 119 • 14 -120 • 46 E) and the south wing of Yangtze River Delta, covering an area of 10,942 km 2 ( Figure 1). Jinhua City belongs to the eastern part of Jinqu Basin, with high terrain in the north and south and low terrain in the middle. It is characterized by a subtropical monsoon climate, with an average annual temperature of 17.5 • C and an average annual rainfall of 1424 mm. The vegetation type in study area is mainly subtropical evergreen broad-leaved forest, and the soil types are diverse, including red soil, yellow soil, purple soil, lithologic soil, coarse bone soil, limestone soil, fluvo aquic soil, and paddy soil. The main land use types are shown in Figure S1. Jinhua is the fourth largest metropolitan area in Zhejiang Province, with a gross domestic product (GDP) of 4.56 × 10 11 CNY and a population of 4.92 million. Among them, Yiwu City has the highest population density (1683 persons/km 2 ), followed by Jindong District (766 persons/km 2 ) and Wucheng District (688 persons/km 2 ). Jinhua has distinctive regional economic characteristics and several major industries, such as the medicine and chemical industry, food, automobile parts, building materials, the light industry, and textiles, all of which cluster in this region.

Sample Collection and Analysis
Soil samples were selected by grid sampling method as suggested by previous study due to its uniform distribution, so as to provide an optimal geostatistical fitting effect [32][33][34]. Therefore, establishing a density of 25 km 2 by using a 5 km × 5 km grid were established to determine the sampling sites and a total of 209 topsoil (0-20 cm) samples were collected from agricultural land (Figure 1). At each site, five subsamples within a 100 m radius from the sampling site were collected and mixed into a composite sample, placed in a polyethylene bag, brought back to the laboratory. The original weight of each sample was guaranteed to exceed 1 kg. Soils were air-dried, ground, and passed through a 0.149 mm sieve for testing the HMs before remove animal and plant residues. Soil pH was measured using a pH meter (Mettler Toledo, Switzerland) with a ratio of 1:2.5 soil to deionized water [35]. Soil Cr, Cd, Cu, Pb, Ni, and Zn were digested with HCl-HNO 3 -HClO 4 on a hot plate and then determined by inductively coupled plasma emission spectrometer (ICP-AES); soil Hg and As were digested by aqua regia-digestion in water bath and measured using an inductively coupled plasma fluorescence spectrometer (ICP-AFS) [36]. All reagents used in the analyses were of premium grade. Quality control was carried out with the national soil standard reference material (GSS-24, Center of National Reference Materials of China), and 20% of the samples were randomly selected to test duplicate samples, the error was kept within 5%, the recovery rate ranged from 81% to 109%, and the measurement results were within the allowable range of error [37].

Absolute Principal Component Score and Multiple Linear Regression
Absolute principal component multiple linear regression is based on the principal component by introducing a variable with zero concentration and subtracting the factor fraction of the principal component from the variable with zero concentration, thus converting the factor fraction of the principal component into an absolute principal component factor fraction (which in turn identifies the pollution sources [22]), and then applying a multiple linear regression model to further determine the contribution of each source [18]. The equation are as follows: where X ik and Z ik represent the concentration and standardized concentration of the element i at sampling site k, respectively; X i and σ represent the average concentration and standard deviation of element i, respectively; j is the number of factors; P jk is the factor j affecting the total HMs in sample site k; P j0 is the influence of factor j on the concentration of sample site k with zero concentration; (APCS) jk is the factor score of transformed sample k; b 0i represents the constant term of multiple linear regression; and b ij represents the coefficient of linear regression of factor j on element i. The average value of b ij (APCS) jk at all sampling sites is used to represent the contribution of HMs sources.

Positive Definite Matrix Factor Analysis
The PMF model is a receptor model based on factor analysis for inferring the number and contribution of each source based on the situation of the pollution source [18,21]. The equation of the PMF 5.0 model is described in the USEPA PMF 5.0 guidance manual [38]. PMF model decomposes the original matrix into the source component spectrum matrix and source contribution rate matrix [12], expressed as: where x ij is the concentration of the jth element in the i samples; G ik is the contribution of kth source to the ith sample; F kj is the concentration of the jth element in the source k; p is the number of source factors; e ij is the residual matrix obtained by minimizing the object function Q, which can be expressed as: where Q is the square sum of the difference (i.e., e ij ) between the original dataset (X ij ) and the PMF output (G ik F kj ); U ij represents uncertainty; RSD is the relative standard deviation; and MDL is the detection limit of the species-specific method.

Health Risk Assessment
Human health risk assessment is estimated based on the exposure factor manual of USEPA, in which individuals are exposed to pollutants mainly through ingestion, respiratory inhalation and dermal contact [39,40]. According to the manual, the daily intake doses for adults and children under these three exposure routes are expressed as follows: where ADD ing , ADD inh , and ADD derm represent the daily intake doses under three different exposure routes: ingestion, respiratory inhalation, and dermal contact, respectively; C soil represents the concentration of soil HMs. The definitions and values of exposure parameters of adult and child intake are shown in Table S1. The noncarcinogenic risk (HI) and carcinogenic risk (RI) are calculated as follows: where HQ is the ratio of the daily intake dose to the corresponding reference dose, which represents the non-carcinogenic risk effect of pollutants under a certain route; HI represents the noncarcinogenic effect of pollutants under all routes; CR is the daily dose multiplied by the corresponding slope factor, which represents the carcinogenic effect of pollutants in a certain way; RI represents the carcinogenic effect of pollutants in all routes. The reference dose (RfD) and slope factor (SF) of the HMs applied in health risk assessment are shown in Table S2 [4,37].

Statistical Analyses
The descriptive statistics, K-S test, and APCS-MLR model were performed using SPSS 25.0 for Windows. The semivariogram modeling of the HMs was analyzed using GS + 9.0 for Windows, and the spatial distribution was analyzed using ordinary kriging interpolation in ArcGIS 10.2 for Windows. USEPA PMF 5.0 was conducted to quantify the HMs sources and contributions.

Descriptive Statistics of Heavy Metals in Soil
The descriptive statistics of 8 HMs in agricultural soil from Jinhua City were summarized in Table 1. The mean soil pH value was 5.34, which indicated moderate acidity in study area. The mean concentrations of As, Cd, Cr, Cu, Hg, Ni, Pb, and Zn in the agricultural soils of Jinhua City are 6.37, 0.18, 36.73, 18.6, 0.11, 11.95, 32.48 and 76.05 mg kg −1 , respectively. Among them, the concentration of Cd, Hg, Pb, and Zn exceeded their corresponding background values in Zhejiang and China, whereas the concentration of As, Cr, and Ni did not. Moreover, the concentration of Cu was higher than its corresponding background values in Zhejiang, while lower than those in China. Particularly, the mean concentrations of 8 HMs in Jinhua were less than their corresponding critical values in comparison to the risk control standard for soil contamination of agricultural land according to the Chinese Soil Quality Standard [41]. The coefficient of variation (CV) of soil HMs in the study area revealed that Hg (152.24%) had strong variability, which demonstrated higher discontinuity and a wider potential for anthropogenic effects [4]. The rest of the HMs were presented at a moderate variability and followed the order of Ni (98.55%) > Cr (72.12%) > As (57.56%) > Cu (56.29%) > Cd (53.07%) > Zn (34.39%) > Pb (25.59%). Compared with agricultural soils in the other regions of the Yangtze River Delta [27][28][29][30][31], the overall concentrations of HMs in the study area are lower ( Figure 2). It is shown that it pollutes the soil environment lightly in the process of economic development in the region. Among them, Cr, Cu, Ni, and Zn have the lowest concentration in the whole delta region. Although the Cd concentration in the study area is lower than that in other cities in Zhejiang Province, it is significantly higher than that in the cities in Jiangsu and Anhui. The Pb concentration is higher than Hangzhou, Ningbo, Suzhou, and Yangzhou, while it is much lower than Shanghai and Nanjing. Therefore, it is worth mentioning that the economic development of the study area belongs to the upper-middle level in this area, but the degree of pollution is relatively low, indicating that the soil environmental planning and governance in this area can provide some reference for other areas.  Figure 3 shows the spatial distribution of HMs concentration in agricultural soils in Jinhua City. The high concentration areas of As is mainly concentrated in the central, northern and western parts of Jinhua City, including the center of Lanxi, Jindong, and Yiwu, eastern Pujiang, and northwestern Wucheng. Generally, As is closely related to industrial emissions [26,43,44]. Many factories were built near the river to facilitate the sewage treatment of factories and enterprises. As may be related to industrial activities such as industrial discharge and sewage sludge in this area, so it can be found that As elements are mainly distributed along the river [45,46]. Similarly, the spatial distribution of Hg is relatively uniform, with sporadic high concentration areas along with the Wujiang River. Studies have shown that coal-burning activities are prone to producing a large amount of Hg, and Hg has the characteristics of volatility and high fluidity, making it easy for it to enter the atmosphere for migration and diffusion [15,47]. The flow of Hg in the atmosphere is less obstructed, and it only enters the soil under weather such as rainfall and has little impact on the soil within a certain range, so the spatial distribution of Hg brought by atmospheric deposition is more uniform [48]. In addition, similar to the distribution of As element, most of the factories in the study area are located close to the river, and the exhaust gas also causes the accumulation of Hg in the area [15]. The highest concentrations of Cr, Cu, and Ni are concentrated in the central and western parts of Jinhua City, including Lanxi, Wucheng, and Jindong. Alongside this, the highest concentrations of Pb and Cd are mainly centered in Wucheng and the junction of Lanxi-Wucheng, respectively (which is the downtown of Jinhua and contains large-scale development zones such as Jinhua Technology and Economic Development Zone, Jinxi Development Zone, etc.). A previous study has demonstrated that the unreasonable discharge of domestic tannery and electroplating enterprises may cause Cr pollution [49], and the discharge of waste liquid from chemical plants usually causes the accumulation of Cu [26,42], which may be one of the reasons for the accumulation of Cr and Cu in this area. Along with this, the Pb and Cd used in dyes of textile printing and dyeing enterprises will eventually enter the soil through sewage treatment [50]. Simultaneously, the agricultural land in this area is widely distributed, and unreasonable agricultural activities, such as excessive use of fertilization, pesticides and herbicides, and sewage irrigation may also lead to the accumulation of Cd and Cu in the soil [51]. Unlike other HMs, the high concentration areas of Zn are mainly distributed in Wuyi. Moreover, Ni, Cu, and Zn are also significantly enriched in the southern part of Wuyi. According to field investigations, Wuyi has large factories and enterprises such as the Hardware Industrial Zone and Huashan Industrial Park, which will inevitably emit a lot of smoke and dust during the production process. Additionally, Wuyi is a well-known fluorite mining base in China, and the industrial sites, transportation roads, and underground mining areas will all have an impact on the surrounding ecological environment, especially the groundwater and soil environment [11]. Firstly, the HMs contained in the waste slag and dust will migrate to the soil surface under the rain leaching of the waste ore [23]; secondly, the waste gas and dust generated in the transportation process of vehicles will inevitably produce the accumulation of zinc and other elements in the surrounding agricultural soil [52]. Soil pH ranges between 4.65-6.46 in the study area (Figure 3i) and the soil pH in Jindong District is weakly acidic, while the soil pH in the southwest of the study area (mainly including Wucheng, Wuyi, and Yongkang) is strongly acidic.  (Table S3). However, there is a lack of information in the method of principal component analysis. Thus, we introduced a zero-concentration variable to calculate the absolute principal component (Table 2). Then, by taking the factor scores after the absolute principal component as the independent variable, taking the standardized pollutant concentration as the dependent variable, and substituting it into the multiple linear regression program, the calculated contribution rates of the four main pollution sources are 43.1%, 27.8%, 15.4%, and 13.7%.

Spatial Distribution of Heavy Metals in Soil
The loadings of the four factors of principal component were listed in Table 1. The variance contribution rate of APCS-MLR_1 is 35.42%. Cr, Ni, and Cu have larger loads on factor 1, and their contribution rates are 89.2%, 86.0%, and 83.4%, respectively. Compared with the descriptive statistical results (Table 1), the contents of these elements are lower than or close to the background value (especially for Cr and Ni, whose contents do not exceed the soil background value of Zhejiang Province and are less affected by human activities). Simultaneously, earlier studies found that these elements are mainly controlled by the geological background [23,48]. In addition, the PMF_1 is consistent with the results obtained by APCS-MLR_1. It also has a higher contribution rate to Ni, Cr, and Cu, which are 70.7%, 66.9%, and 40.0%, respectively. Therefore, the APCS-MLR_1 and PMF_F1 can be regarded as a "natural sources factor". The variance contribution rate of APCS-MLR_2 is 22.82%. Pb and Cd have larger loads on APCS-MLR_2, with contributions of 84.6% and 78.3%, respectively, which were higher than the corresponding background values of Zhejiang Province, and affected by human factors. An earlier study has shown that soil pH has an important effect on the migration and transformation of Cd elements in the soil. Usually, the solubility of Cd decreases with the increases of pH value, which makes it difficult for the Cd in the soil to migrate and deposit in situ [10]. Figure 3 also shows that the distribution of high concentration area of Cd basically consistent with the high value area of soil pH. The use of fertilizers, herbicides, pesticides, and the discharge of wastewater from factories and enterprises are important ways to artificially change the pH value of the soil environment, which in turn affects the migration and accumulation of corresponding soil heavy metals in the soil. Therefore, Cd is indicated to mainly come from the application of the chemical fertilizers. Although the application of chemical fertilizers can effectively supplement the nutrients required for crop growth and improve the quality and yield of food, excessive input will cause the accumulation of Cd in the soil [11]. Besides, Cd may also come from the "three wastes" emissions from electroplating and metallurgy industries [51]. It has been demonstrated that Pb mainly comes from automobile exhaust and coal combustion [24]. Therefore, the APCS-MLR_2 can be regarded as a "mixed pollution sources of industrial, agriculture, and transportation". The PMF_2 has a higher contribution rate to Cd (69.7%) and Cu (53.0%), and the mean concentration of Cd and Cu are higher than the corresponding background values of Zhejiang Province and their CV revealed a moderate variability, indicating that they are affected by human factors. Previous studies have shown that chemical fertilizers, pesticides, and plastic films used for crop insulation in agricultural production are important reasons for the accumulation of Cu and Cd in agricultural soil [48,52]. It has been proven that the application of nitrogen and phosphate fertilizers may cause the accumulation of Cd in agricultural soil [51]. Most pesticides contain large amounts of HMs. For example, herbicides and pesticides contain a large amount of Cu [5]. Moreover, the excrement of livestock and poultry are also the important sources of HMs in agricultural soil [21]. Cu and other elements are usually added as additives to the feed for raising livestock and eventually penetrate into the soil through the residual infiltration of the food chain and cause soil pollution [3]. Therefore, PMF_2 can be regarded as the "agricultural sources factor".
The variance contribution rate of APCS-MLR_3 is 12.67%, with the highest contribution to As (81.2%), and its mean concentration is higher than the background value of Zhejiang Province. Earlier studies have proven that As mainly comes from industrial pollution emissions [26,42]. There are 3586 industrial enterprises in the study area, and As is closely related to industrial activities. PMF_3 also has a high contribution rate to As, and its contribution rate is 68.7%. Thus, here the same explanation applies. Therefore, the PCAS-MLR_3 and PMF_3 can be characterized as the "industrial source factor".
The variance contribution rate of APCS_4 is 11.25%, and has the higher contribution rate to Hg (78.7%). Previous studies have shown that the burning of fossil fuels such as coal will release a large amount of Hg [15,47]. According to the environmental survey of the study area, there are many coal-fired powers plants and other industrial areas in the study area. During the coal-burning operation of these factories, a large amount of Hg will inevitably be released into the atmosphere and finally sink into the soil under the weather such as rainfall [47]. The result of PMF_4 is consistent with the result of APCS-MLR_4, The contribution rate is 84.2%. Therefore, all fourth factors can be represented as the "atmospheric deposition sources factor".
The PMF_5 has a higher contribution rate to Pb (74.9%) and Zn (63.1%). Studies have shown that Pb and Zn pollution are closely related to traffic [53]. In order to improve the antiknock performance of gasoline, "lead water" (tetraethyl lead) is usually added to gasoline and the leaded exhaust pollutants emitted by vehicles will eventually enter the soil [14]. Zn is an important part of antioxidants in automobile tires [53]. Dust generated by the wear and tear of automobile tires will also fall into the surface soil, making Zn enter the soil [4]. The closer it is to the road, the higher corresponding HM accumulation. The highway passenger traffic volume of Jinhua City is 107.73 million people, ranking second in Zhejiang Province, and the freight volume is 99.11 million tons. The huge traffic volume inevitably releases a large amount of HMs. Therefore, it can be concluded that PMF_5 can represent the traffic pollution sources factor.

Comparison of APCS-MLR and PMF
Compared with PMF model, the number of sources identified by the APCS-MLR is reduced (Table 3). In addition, APCS-MLR does not have a good analysis of the source of APCS-MLR_2, which contains mixed sources of agriculture, transportation, etc. This might be attribute to the APCS-MLR excludes some variables with poor publicity in the calculation process [22]. There are also certain differences in the contribution rate between two models. For example, in terms of natural sources, the contribution rate obtained by APCS-MLR is 43.1%, while the contribution rate obtained by PMF model analysis is 24.29%. The reason for the difference in the contribution rate between APCS-MLR and PMF model is that the principal component does not consider the uncertainty of the data in the calculation process, and there is a small amount of information loss [54]. In the PMF model, the uncertainty data of soil HMs and non-negative constraints were added in the calculation process, which made the result more reasonable [54]. Therefore, PMF can provide a better result than APCS-MLR.
Some previous studies have shown that the sources apportioned by the PMF model are more reliable than those apportioned by the APCS-MLR [55][56][57], while others are not. The results of this study show that APCS-MLR cannot accurately identify some similar sources. For example, the pollution source 2 resolved by APCS-MLR is a mixed source, while the PMF model can identify more pollution sources. Therefore, it may be concluded that the PMF model has better explanatory power for areas with few pollution sources and light pollution and has stronger credibility for areas with more complicated pollution, which is basically consistent with the research of Wu et al. [58]. Overall, in the specific soil pollution investigation, appropriate receptor models were chosen based on the size of the research area and the soil pollution situation. For example, the PMF model was more suitable in a large-scale urban soil pollution source survey, while APCS-MLR model was more suitable for a small-scale survey of pollution sources such as factories and enterprises.
We evaluated the suitability of each soil HM with the model by comparing the measured concentration with the predicted concentration, as well as the coefficient of determination (R 2 ) and prediction error [18]. Overall, both models show a higher fitting ability, and that of the PMF is higher than that of APCS-MLR (Table 4). The R 2 of HMs in APCS-MLR are all greater than 0.71, and the R 2 of PMF model are greater than 0.81. The prediction errors of the two models are less than 1, except for Pb in APCS-MLR model. For the APCS-MLR model, all elements are overestimated. In comparison, the PMF model overestimates the As, Cd, and Ni, and underestimates the Cr, Cu, and Pb.  In general, the analysis results of the APCS-MLR and PMF models are roughly the same, and both can be applied to the source analysis of soil HMs (Figure 4). However, in terms of the contribution rate of source analysis, the PMF model has greater advantages than APCS-MLR in this study.

Carcinogenic and Non-Carcinogenic Risk
The RI and HI assessment results of adults and children under different exposure routes are shown in Tables S4 and S5. The RI of the three different exposure routes are quite different and followed the rank order: nondietary ingestion > dermal contact > respiratory inhalation. The RI of different HMs are also significantly different, and the levels of RI were ranked as follows: As > Cd > Cr. Among them, the RI of Cr element is less than 10 -6 , indicating that Cr will not pose a threat to the health of residents for the time being. The RI of As and Cd are between 10 −4 -10 −6 , which is within an acceptable range. However, the total RI of multiple HMs to children is 1.30 × 10 -4 , which exceeds the RI threshold, and relevant environmental departments thus need to pay attention. The HI of the three exposure routes in study area are significantly different. Among them, the value of nondietary ingestion is the largest, followed by dermal contact, and then respiratory inhalation as the smallest. The HI of the 8 HMs followed the rank order: As > Cr > Pb > Ni > Cu > Hg > Zn > Cd. Among them, the HI of children is much greater than that of adults, but both are less than 1, indicating that the health risk of HMs in agricultural soil in the study area is within an acceptable range and will not cause HI to surrounding residents. A comparative analysis of the RI and HI of HMs in agricultural soils in the study area shows that children have a higher health risk than adults and are more susceptible to soil pollution hazards. Combined with previous studies, it can be found that although the health risks of soil HMs varied greatly in different regions, the RI and HI of children are usually greater than those of adults [59][60][61][62]. The main reason is that children's behavior patterns determine that they are more likely to be exposed to the soil environment, especially dermal contact [39], and children are more sensitive to environmental pollution and have lower immunity, which make them more vulnerable to soil pollution hazards [63]. Therefore, it is necessary to instruct children to wash their hands frequently and correctly, maintain hand hygiene, and change laundry frequently [4]. Secondly, it can be seen that whether it is children or adults, the RI and HI of nondietary ingestion are greater than the two exposure routes of respiratory inhalation and dermal contact, which is basically consistent with the results of previous studies [64,65].

Health Risk Assessment from Different Sources
The PMF model was performed to quantify the different HMs sources to RI and HI ( Figure 5). It can be found that RI and HI of different HMs sources have the similar trend in the risk probability of adults and children. Industrial sources have the highest HI to adults (32.92%) and children (30.47%), followed by agricultural sources (21.34% and 29.31%), natural sources (26.26% and 18.56%), traffic sources (15.47% and 17.04%), and atmospheric deposition sources (4.02% and 4.62%). Therefore, the control of Cd, Cu, and As from industrial and agricultural sources is the key to preventing HI for adults and children. Regarding the RI, industrial sources (60.74% and 61.50%) and agricultural sources (32.94% and 33.19%) contributed the highest RI to adults and children. It should be pointed out that in this study, industrial sources (13.42%) accounted for a small proportion of the contribution rate of soil HM pollution sources. However, the health risks caused by it are high whether it is RI or HI. This may be that the main HMs from industrial sources such as have strong toxicity and high slope coefficient [65]. Therefore, we can find that the health risks caused by higher pollution sources are not necessarily high, which is basically consistent with the previous conclusions [4,24,25]. Therefore, the control of industrial sources and agricultural sources are the key to prevention for adults and children. The noncarcinogenic risk values for adults and children based on different pollution sources are less than 1, indicating that the noncarcinogenic risk caused by each pollution source is basically negligible; the carcinogenic risk for adults varies among different pollution sources, with the carcinogenic risk from natural sources, agricultural sources, atmospheric deposition sources, and traffic sources less than 10 −6 , and that of industrial sources between 10 −6 and 10 −4 (1.05 × 10 −5 ). There are also large differences in the carcinogenic risks to children from different pollution sources, with natural sources having carcinogenic values less than 10 −6 , atmospheric deposition sources and traffic sources having values of 2.14 × 10 −5 and 4.40 × 10 −6 , respectively, and agricultural and industrial sources having values of 4.31 × 10 −5 and 7.98 × 10 −5 , respectively. The USEPA guidelines indicate that when CR < 10 6 , the carcinogenic risk is basically negligible; many scholars also consider that when 10 −6 < CR < 10 −4 , the carcinogenic risk is within the acceptable range; when CR > 10 −4 , it indicates that there may be a greater carcinogenic risk [59][60][61][62]. In this study, a CR value between 10 −6 and 10 −5 was defined as low risk, a CR value between 10 −5 and 10 −4 was defined as medium risk, and a CR value > 10 −4 was defined as high risk based on the actual local situation. Therefore, the carcinogenic risk to children from industrial and agricultural sources in the study area is medium, and the risk from atmospheric deposition sources and traffic sources is low. However, the carcinogenic risk to adults from industrial sources is low, the risk from agricultural sources is medium, and the carcinogenic risk from other pollution sources is basically negligible ( Table 5). The industrial agglomeration in the study area is obvious, and there are many industrial zones, such as apparel factories, rubber factories, and plastic industries. Therefore, the overall cancer risk from industrial sources in the study area is medium to low, but since the spatial distribution of HMs pollution is mostly concentrated in densely populated areas, the subsequent supervision and management of factory enterprises with high pollution and discharge should also be strengthened to reduce the possible health risks caused by soil pollution in the study area.

Conclusions
This study demonstrates the spatial variability, potential sources of HMs in agricultural soils of subtropical China. Then, an integrated approach combining human health risk assessment model with PMF model was proposed to quantitative human health risks from identified sources of soil HMs and determine priority control sources in this study. The mean concentrations of Cd, Cu, Hg, Pb, and Zn in the soil were higher than the corresponding background values in Zhejiang Province, while the mean concentrations of determined 8 HMs were less than their corresponding risk-screening values for soil contamination of agricultural land. The high concentration of As and Hg are roughly distributed along the Wujiang River, while the Cu, Ni, Cd, and Pb accumulated in large patches in the central and western parts of the Jinhua City. Besides, Cr, Ni, Cu, and Zn are also accumulated in the southern of Wuyi County. PMF model performed better than APCS-MLR model in the identification of major sources of soil HMs as it revealed higher R 2 value (0.81-0.99) and lower prediction error (−0.93-0.25%). PMF model also had higher accuracy than APCS-MLR model in the estimation of source contribution, and apportioned five sources, namely, agricultural sources (30.06%), natural sources (24.29%), traffic sources (18.9%), industrial sources (13.42%), and atmospheric deposition sources (13.27%), respectively. Both RI and HI of adults are below the safety threshold, while the RI of children has exceeded the safety threshold, indicating that local children are already facing the threat of RI, and needs to be concerned by relevant departments. Moreover, a PMF based human health risk assessment model indicated that industrial sources contributed the highest risk to HI (32.92% and 30.47%) and RI (60.74% and 61.50%) for adults and children, followed by agricultural sources (21.34%, 29.31% and 32.94%, 33.19%). Therefore, integrated environmental management should be implemented to control and reduce the accumulation of soil HMs from these two sources. This study provided an effective approach to quantify the priority pollution sources of concern for human health risk, which is of great significance for pollution control and risk reduction under limited resources.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/ 10.3390/land10101016/s1, Figure S1: land use map of the study area, Table S1: values of exposure parameters of the health risk assessment model, Table S2: reference doses and carcinogenic heavy metal slope factors for different exposure routes, Table S3: the eigenvalues of principal component analysis and their contribution rate, Table S4: different exposure routes and total carcinogenic health risk assessment results in the study area, Table S5: different exposure routes and total noncarcinogenic health risk assessment results in the study area.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this paper are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.