Source Apportionment of Heavy Metals Based on Multiple Approaches for a Proposed Subway Line in the Southeast Industrial District of Beijing, China

Apportioning the sources of heavy metals (HMs) in soil is of great importance for pollution control. A total of 64 soil samples from 13 sample points at depths of 0–21 m were collected along a proposed subway line in the southeast industrial district of Beijing. The concentrations, distribution characteristics, and sources of eight HMs were investigated. The results showed that the concentrations of Hg, Cd, Cu, Pb, As, and Zn in the topsoil (0–2 m) exceeded the Beijing soil background values. Three sources were identified and their respective contribution rates calculated for each of the HMs using multiple approaches, including correlation analysis (CA), top enrichment factor (TEF), principal component analysis (PCA), and positive matrix factor (PMF) methods. As (63.11%), Cr (61.67%), and Ni (70.80%) mainly originated from natural sources; Hg (97.0%) was dominated by fossil fuel combustion and atmospheric deposition sources; and Zn (72.80%), Pb (69.75%), Cu (65.36%) and Cd (53.08%) were related to traffic sources. Multiple approaches were demonstrated to be effective for HM source apportionment in soil, whilst the results using PMF were clearer and more complete. This work could provide evidence for the selection of reasonable methods to deal with soils excavated during subway construction, avoiding the over-remediation of the soils with heavy metals coming from natural sources.


Introduction
Heavy metal pollution in urban soil has become an important threat to both the soil environment and human health in China [1][2][3][4]. This situation is more serious in the old industry districts because of the intense human activities [5][6][7]. In recent years, many old industry districts around the city have been shut down and relocated. Before site redevelopment, soil pollution status must be investigated, and corresponding risk control measures or remediation activities must be taken [8]. As is well known, HMs in the soil originate both from human activities (e.g., mineral resource development, metal processing and smelting, fossil fuel combustion, vehicle exhaust emissions, and wastewater irrigation), and from natural processes (e.g., weathering and pedogenesis processes, and precipitation) [9][10][11][12]. It is unreasonable and expensive to blindly carry out remediation activities regardless of the pollution sources. Thus, distinguishing various sources of HMs in soil is crucial for site remediation and environmental management.
The existing methods of analyzing sources of heavy metal pollution include geostatistical analysis, multivariate statistical analysis, and receptor model analysis methods. Geostatistical and multivariate analyses can only qualitatively identify the source of pollution [12][13][14]. Receptor models, such as chemical mass balance (CMB) [15][16][17], positive matrix factorization (PMF) [3,10,[18][19][20], absolute principal component analysis-multivariate linear regression (APCA-MLR) [21][22][23][24], and the UNMIX model [25][26][27][28][29][30] can identify and quantify source contributions based on the chemical and physical characteristics of pollutants in sources and receptors. The CMB model is mainly used for source apportionment of atmospheric particulate matter [15][16][17]. The PMF model decomposes the original dataset into a contribution matrix and a profile group for contribution calculation and source identification [31], and has been widely used to determine the sources of pollutants in soil [3,10,[18][19][20]. Guan et al. [18] used multivariate statistical analysis, including GIS-based geostatistical methods and PMF, to apportion sources of HMs in agricultural soil in the Hexi Corridor. Wu et al. [19] observed that the partition computing-based PMF approach can, to a certain extent, overcome the uncertainty of the spatial heterogeneity for source apportionment. APCA-MLR identifies the main pollution sources by reducing the dimension of the PCA data, and then determining the contribution rate of each pollution source using its absolute factor score and multiple linear regression of pollutant content [24]. Jin et al. [21] employed the APCA-MLR receptor model to identify pollution sources of heavy metals in farmland soil. Zhang et al. [22] suggested that the PMF approach could provide a more physically plausible source apportionment in the study area, and a more realistic representation of groundwater pollution, than the solutions from the APCA-MLR model. The UNMIX model was first applied to the source apportionment of atmospheric particulate matter [30], and has also been used to analyze heavy metal pollution sources in soil in recent years [25][26][27][28][29]. In addition, tracer techniques [32][33][34] and binding forms analysis [35,36] have also been applied to source apportionment. Wang et al. [32] revealed that the Bayesian stable isotope analysis mixing model (MixSIAR) could reduce the uncertainty of the PMF model and allocate the contributions of different sources more accurately. Sungur et al. [35,36] believe that HMs in uncontaminated soil are mainly immobile in forms bound to silicates and primary minerals, and that the anthropogenic and lithosphere sources of HMs are distinguished by analyzing different geochemical fractions. However, most of these studies are concerned with agricultural soil at a regional scale, and are mainly focused on the comparison of different approaches. The suitability of multiple approaches for source apportionment of soil pollution at site scale remains a challenge. Meanwhile, little attention has been devoted to the application of source apportionment results to site remediation and environmental management.
The Beijing southeast industrial district used to be one of the oldest industrial areas. Many highly polluting industries, such as coking and chemical plants, were scattered here. Severe pollution levels have been observed in this area [37,38]. According to the development plan, a new subway line will pass through this area. It is important to investigate the soil pollution conditions along the subway corridor and to identify the potential sources and their contributions to the soil contamination. The specific objectives of this study were (1) to investigate the content and distribution characteristics of eight HMs in this study area, (2) to explore the sources of different HMs using multiple approaches and to quantify the source contributions through PMF modeling, and (3) to provide a reference for soil remediation and environmental management of this subway line.

Study Area and Sample Collection
The proposed subway line, with a total length of 8.8 km, is located in the old chemical industrial area of Chaoyang district, Beijing ( Figure 1). It consists of eight stations and seven depots. Many heavily polluting industrial enterprises were distributed along the subway corridor; including pharmaceutical, coking, chemical manufacturing, dyeing, and glassmaking plants, among others. All of these industrial enterprises have been shut down and relocated. Thirteen sampling points were chosen along the planned subway line, spaced as evenly as possible with an approximate interval of 600 m, as shown in Figure 1. line, spaced as evenly as possible with an approximate interval of 600 m, as shown in Figure 1. A total of 64 soil samples were collected, based on the pollution location and formation lithology, including 13 topsoil samples (0-2 m) and 51 subsoil samples (2-21 m). Soil profile samples were collected using a cable percussive drilling method. Percussive drilling makes use of the gravity and downward impulse of drilling tools to cause the drill to hit the bottom of the hole, so as to break up the soil layer and realize drilling. During drilling, the soil samples and cores taken by the core pipe were stored in a core box in sequence, and soil profile samples were taken from the fill layer (0.7-1.9 m), the upper weak permeable layer (2.2-6.6 m), the weak permeable layer (7.1-10.7 m), the aquifer (9.7-17.6 m), and the impermeable layer (16.6-20.9 m) within each borehole.

Sample Preparation and Analysis
All soil samples were naturally air-dried in a ventilated and dark laboratory, sieved into 60-mesh size particles after removing gravel, residual roots, and other unwanted materials, and then sealed in brown glass bottles and conserved in a refrigerator at 4 ℃. As, Be, Cd, Cr, Cu, Ni, Pb, V, and Zn were analyzed using ICP-MS (Agilent 7500, USA). A 1.0 g aliquot of soil sample was weighed into a 50 mL polyethylene digestion tube, 4 mL (1 + 1) nitric acid and 10 mL (1 + 4) hydrochloric acid were added, and the mixture heated at 95 °C for 45 min. It was then slightly cooled and made up to volume with deionized water, shaken well and stood still before measurement.
Hg was determined using an atomic fluorescence photometer (AFS; Beijing Jitian Instruments Co., Ltd. Production, Beijing, China, AFS-820). A 0.250 g soil sample was weighed into a 25 mL colorimetric tube, 10 mL of aqua regia (HNO3:HCl = 3:1) was added, A total of 64 soil samples were collected, based on the pollution location and formation lithology, including 13 topsoil samples (0-2 m) and 51 subsoil samples (2-21 m). Soil profile samples were collected using a cable percussive drilling method. Percussive drilling makes use of the gravity and downward impulse of drilling tools to cause the drill to hit the bottom of the hole, so as to break up the soil layer and realize drilling. During drilling, the soil samples and cores taken by the core pipe were stored in a core box in sequence, and soil profile samples were taken from the fill layer (0.7-1.9 m), the upper weak permeable layer (2.2-6.6 m), the weak permeable layer (7.1-10.7 m), the aquifer (9.7-17.6 m), and the impermeable layer (16.6-20.9 m) within each borehole.

Sample Preparation and Analysis
All soil samples were naturally air-dried in a ventilated and dark laboratory, sieved into 60-mesh size particles after removing gravel, residual roots, and other unwanted materials, and then sealed in brown glass bottles and conserved in a refrigerator at 4 • C. As, Be, Cd, Cr, Cu, Ni, Pb, V, and Zn were analyzed using ICP-MS (Agilent 7500, Santa Clara, CA, USA). A 1.0 g aliquot of soil sample was weighed into a 50 mL polyethylene digestion tube, 4 mL (1 + 1) nitric acid and 10 mL (1 + 4) hydrochloric acid were added, and the mixture heated at 95 • C for 45 min. It was then slightly cooled and made up to volume with deionized water, shaken well and stood still before measurement.
Hg was determined using an atomic fluorescence photometer (AFS; Beijing Jitian Instruments Co., Ltd. Production, Beijing, China, AFS-820). A 0.250 g soil sample was weighed into a 25 mL colorimetric tube, 10 mL of aqua regia (HNO3:HCl = 3:1) was added, then it was shaken well and cold digested overnight. After placing in a boiling water bath for 2 h the next day, the volume was made up with deionized water, shaken well and stood still before measurement. In quality assurance and quality control (QA/QC) studies, the recovery of a spiked blank ranged from 80% to 112% within the control range of 80-120%, and the relative percentage differences (RPD) of duplicate samples ranged from 1% to 16% within the control range of 0-20%. The RPD of spike matrix duplicates was 0-16% within the control range of 0-20%. These results established that the precision and accuracy of our analytical method were acceptable.

Top Enrichment Factor
The TEF value is the concentration ratio between the topsoil and subsoil [12]. It is a method to differentiate the topsoil metals originating primarily from anthropogenic or natural sources. The TEF is calculated as per Equation (1) [13]: where C i-topsoil is the concentration of target metal i in topsoil, and C i-subsoil is the concentration of target metal i in subsoil.

Principal Component Analysis
PCA was carried out to identify the HM sources. Varimax and Kaiser normalization rotation was applied because orthogonal rotation minimized the number of variables with a high loading of each component and facilitated the interpretation of results. The data analyses were performed using SPSS Statistics 22.0 (IBM Corp., Armonk, NY, USA). The biplot was drawn with OriginPro 2018C (OriginLab Corporation, Northampton, MA, USA).

Positive Matrix Factorization Method
PMF relies on a receptor model that quantifies the contribution of sources in samples based on the composition or fingerprints of the sources. In this study, PMF version 5.0 (U.S. EPA, Washington, DC, USA) was used for source apportionment. This model decomposes the original matrix x ij into two factor matrices, g ik and f kj , and a residual matrix e ij ; the basic equation is as Equation (2): where x ij is the content of the jth HM in sample i; g ik is the contribution of the kth source for sample i; and f kj is the source profile of the jth HM in source k. The residual error matrix e ij is calculated as the minimum value of the objective function Q: where u ij refers to the uncertainty of the jth HM in i number of samples. The main feature of PMF is that the model requires the HM content of samples and the uncertainties of the content that are used to analyze the quality of the content values individually. The uncertainty can be calculated through various methods [39]. In this study, the uncertainty was calculated using Equations (4) and (5) as follows: where c ij is the content of the element, mg/kg; θ is the error fraction; and MDL is the method detection limit.

Descriptive Statistics Analysis
The basic descriptive analysis statistics are presented in Table 1. The coefficients of variation (CVs) also indicated that the regional differences in contaminant content of the topsoil (0. 16-2.42, with low to high variation) were stronger than those in the subsoil (0.32-0.83, with median to high variation), especially for Hg. Therefore, the HM contamination in topsoil was heavier than that for the subsoil.

Distribution Characteristics
The vertical distributions of the eight HMs in the soil profile are presented in Figure 2. Concentrations of Cd, Cu, Pb, Zn, and Hg were significantly accumulated at the depth of 0-2 m, indicating that they had been affected by human activity [12]. The concentrations of Cr and Ni appeared to show no significant differences across the entire soil profile. As had higher concentrations at depths of approximately 10 m, which was exactly at the location of the groundwater level in this area. Accornero et al. [40] also found an accumulation of As in deep soil within the scope of our study area, and various geochemical processes were believed to be the main causes, such as the variation in soil composition, percentage and mineralogy of clays, organic carbon content, variation in the pH, exchange capacity of the soil, presence of either saturated or unsaturated layers, and compositional characteristics of the groundwater. Therefore, the HM content in the subsoil was unlikely to have been affected by human activity, and source apportionment was only carried out for the topsoil.
2. Concentrations of Cd, Cu, Pb, Zn, and Hg were significantly accumulated at the depth of 0-2 m, indicating that they had been affected by human activity [12]. The concentrations of Cr and Ni appeared to show no significant differences across the entire soil profile. As had higher concentrations at depths of approximately 10 m, which was exactly at the location of the groundwater level in this area. Accornero et al. [40] also found an accumulation of As in deep soil within the scope of our study area, and various geochemical processes were believed to be the main causes, such as the variation in soil composition, percentage and mineralogy of clays, organic carbon content, variation in the pH, exchange capacity of the soil, presence of either saturated or unsaturated layers, and compositional characteristics of the groundwater. Therefore, the HM content in the subsoil was unlikely to have been affected by human activity, and source apportionment was only carried out for the topsoil.  Figure 3 shows the horizontal distributions of the six HMs with contents exceeding the background values in the topsoil. For Hg, the samples showing concentrations higher than the background levels were evenly distributed across the study area, while those for the other five HMs were mainly concentrated at Points 1, 6, 9, and 12. Points 1, 6, and 12 were located near the former industrial enterprises, indicating that they might have been affected by production activities. However, Point 9 is far away from the production enterprises and located in the parking lot of a large playground. The organic carbon content in the topsoil was 15.3%, which was significantly higher than the average level (3.4%). Heavy metals in soils can complex and chelate with organic matter to fix them in soil [41], which might result in the enrichment of HMs at Point 9.  Figure 3 shows the horizontal distributions of the six HMs with contents exceeding the background values in the topsoil. For Hg, the samples showing concentrations higher than the background levels were evenly distributed across the study area, while those for the other five HMs were mainly concentrated at Points 1, 6, 9, and 12. Points 1, 6, and 12 were located near the former industrial enterprises, indicating that they might have been affected by production activities. However, Point 9 is far away from the production enterprises and located in the parking lot of a large playground. The organic carbon content in the topsoil was 15.3%, which was significantly higher than the average level (3.4%). Heavy metals in soils can complex and chelate with organic matter to fix them in soil [41], which might result in the enrichment of HMs at Point 9.

Correlation Analysis
Spearman's rank correlation coefficient was adopted to explore the relationships among the eight HMs in topsoil. As shown in Table 2, Cu, Pb, Zn, and Hg were significantly correlated with each other; Cr, Ni, and As were significantly correlated with each other; and Cd was significantly correlated both with As and with Cu, Pb, Zn, and Hg. The significant correlations indicate that these metals probably originated from the same source [12].

Top Enrichment Factor Analysis
As shown in Figure 4 and Table A1, the TEF values for As, Cr, and Ni were lower than 2 in all soil samples. Generally, a TEF value below 2 indicates a natural enrichment process [42]. A few sampling points for Cd, Cu, Pb, Zn, and Hg had TEF values greater than 2; mainly concentrated at Points 1 and 6 (both located near the chemical plant), and Point 12 (located near the dye plant), indicating that the topsoil might have been affected by anthropogenic processes [13].

Correlation Analysis
Spearman's rank correlation coefficient was adopted to explore the relationships among the eight HMs in topsoil. As shown in Table 2, Cu, Pb, Zn, and Hg were significantly correlated with each other; Cr, Ni, and As were significantly correlated with each other; and Cd was significantly correlated both with As and with Cu, Pb, Zn, and Hg. The significant correlations indicate that these metals probably originated from the same source [12].

Top Enrichment Factor Analysis
As shown in Figure 4 and Table A1, the TEF values for As, Cr, and Ni were lower than 2 in all soil samples. Generally, a TEF value below 2 indicates a natural enrichment process [42]. A few sampling points for Cd, Cu, Pb, Zn, and Hg had TEF values greater

Principal Component Analysis
PCA was applied to analyze the sources of the HMs in topsoil. The first two principal components (PCs) explained 89.22% of the original information, with eigenvalues greater than 1 (Table A2).
As shown in Figure 5 and Table A3, Hg (0.927), As (0.852), Ni (0.897), Cd (0.819), and Cr (0.725) had higher contributions to PC1, which was mainly expressed at Point 1-that located near the chemical plant. The concentrations of Hg, As, and Cd exceeded the background values in 76.92%, 15.38%, and 38.46% of the topsoil samples, respectively. However, Ni and Cr showed the opposite, and their TEF values were all below 2. Therefore, PC1 is considered to be related to mixed sources.
Zn (0.988), Cu (0.933), and Pb (0.920) had higher contributions to PC2, which was mainly expressed at Points 12 (located near the dye plant) and 6 (located near the chemical plant). The concentrations of Cu, Pb, and Zn exceeded the background values in 30.77%, 23.08% and 7.69% of the topsoil samples, respectively, and were significantly accumulated in the topsoil. Therefore, PC2 was considered to be related to anthropogenic source. than 2; mainly concentrated at Points 1 and 6 (both located near the chemical plant), and Point 12 (located near the dye plant), indicating that the topsoil might have been affected by anthropogenic processes [13].

Principal Component Analysis
PCA was applied to analyze the sources of the HMs in topsoil. The first two principal components (PCs) explained 89.22% of the original information, with eigenvalues greater than 1 (Table A2).
As shown in Figure 5 and Table A3, Hg (0.927), As (0.852), Ni (0.897), Cd (0.819), and Cr (0.725) had higher contributions to PC1, which was mainly expressed at Point 1-that located near the chemical plant. The concentrations of Hg, As, and Cd exceeded the background values in 76.92%, 15.38%, and 38.46% of the topsoil samples, respectively. However, Ni and Cr showed the opposite, and their TEF values were all below 2. Therefore, PC1 is considered to be related to mixed sources.
Zn (0.988), Cu (0.933), and Pb (0.920) had higher contributions to PC2, which was mainly expressed at Points 12 (located near the dye plant) and 6 (located near the chemical plant). The concentrations of Cu, Pb, and Zn exceeded the background values in 30.77%, 23.08% and 7.69% of the topsoil samples, respectively, and were significantly accumulated in the topsoil. Therefore, PC2 was considered to be related to anthropogenic source.

Positive Matrix Factorization
The PMF model was employed to further determine and quantify the pollution sources and their contribution rates to the contaminants in the soils. The start seed number was chosen randomly, and the number of runs was set to 20. When the factor number was set to 3, the residual error was between −3 and 3, R2 was between 0.66 and 1.00, and the best fit results were obtained. Thus, the factor number was set to 3. The results of the factor profiles and source contributions of the HMs are shown in Figure 6.

Positive Matrix Factorization
The PMF model was employed to further determine and quantify the pollution sources and their contribution rates to the contaminants in the soils. The start seed number was chosen randomly, and the number of runs was set to 20. When the factor number was set to 3, the residual error was between −3 and 3, R2 was between 0.66 and 1.00, and the best fit results were obtained. Thus, the factor number was set to 3. The results of the factor profiles and source contributions of the HMs are shown in Figure 6.
The first factor accounted for 31.63% of the total variance, featuring high loadings of As (63.11%), Cr (61.67%), and Ni (70.80%). Many studies have reported that Cr, Ni, and As in soil would originate from soil parent materials [7,43,44]. Weathering processes of parent soil materials have been demonstrated to be a significant source for Cr, Ni, and As in both topsoil (0-20 cm) and subsoil (120-180 cm) in the Beijing urban area [45]. Sun et al. [7] found that parent materials and pedogenic processes were major factors in the levels and distribution of Ni, Cr, and As in soils in Tangshan. Ni and Cr are iron group elements; they can substitute for each other and create symbiosis enrichment in minerals [46]. Multiple studies have shown that Ni and Cr are derived from parent materials [47,48]. Accornero et al. [40] reported that the site-specific upper baseline concentration (UBC) of As in the study area is 10.4-12.6 mg/kg, which is higher than the wider scale background value for the Beijing area. In our study, As, Cr, and Ni were significantly correlated with each other, and their TEF values were all below 2. Only one sample had an As concentration (13.24 mg/kg) slightly greater than the UBC; this was located at a depth of 1.1 m at Point 1. It seems reasonable to attribute Factor 1 as a natural source.

Positive Matrix Factorization
The PMF model was employed to further determine and quantify the pollution sources and their contribution rates to the contaminants in the soils. The start seed number was chosen randomly, and the number of runs was set to 20. When the factor number was set to 3, the residual error was between −3 and 3, R2 was between 0.66 and 1.00, and the best fit results were obtained. Thus, the factor number was set to 3. The results of the factor profiles and source contributions of the HMs are shown in Figure 6. The first factor accounted for 31.63% of the total variance, featuring high loadings of As (63.11%), Cr (61.67%), and Ni (70.80%). Many studies have reported that Cr, Ni, and As in soil would originate from soil parent materials [7,43,44]. Weathering processes of parent soil materials have been demonstrated to be a significant source for Cr, Ni, and As in both topsoil (0-20 cm) and subsoil (120-180 cm) in the Beijing urban area [45]. Sun et al. [7] found that parent materials and pedogenic processes were major factors in the levels and distribution of Ni, Cr, and As in soils in Tangshan. Ni and Cr are iron group elements; they can substitute for each other and create symbiosis enrichment in minerals [46]. Multiple studies have shown that Ni and Cr are derived from parent materials [47,48]. Accornero et al. [40] reported that the site-specific upper baseline concentration (UBC) of As The second factor accounted for 16.04% of the total variance, featuring high loading of Hg (97.0%). Similar results have been reported in which Hg exhibits a different pattern from the other HMs and is classified solely as a principal component in the source analysis [7,45,49,50]. About 1960 tons of Hg have been discharged to the environment every year from anthropogenic sources such as gold mining, coal combustion, and the production of non-ferrous metals [51,52]. Hg is a volatile component in coal; it is emitted in its oxidation state, particle state, and elemental state forms during the coal combustion process. The dust collector only had a high removal effect for the oxidation and particle state forms, while the elementary Hg would escape into the atmosphere and has entered into the soil through atmospheric transport and deposition [53,54]. A high mercury content caused by long-term coal combustion has also been confirmed in the urban soils of Beijing [55]. The area studied in this work was the earliest chemical industrial district in Beijing, with large-scale industrial enterprises such as coking and chemical plants which used coal as their main fuel, gathered in this region. The long-term production processes would inevitably lead to the escape and deposition of mercury. The topsoil in the studied area had a significantly enriched Hg content, with 76.92% of the topsoil samples having a Hg content exceeding the background value ( Figure 2). Therefore, coal combustion was likely responsible for the enrichment of Hg in this area, and Factor 2 was attributed to the source of fossil fuel combustion and atmospheric deposition.
The third factor accounted for 52.33% of the total variance, featuring high loadings of Zn (72.80%), Pb (69.75%), Cu (65.36%), and Cd (53.08%). Previous studies have shown that soil Cd, Zn, Pb, and Hg levels are tightly associated with human activities [7,45,[56][57][58]. Agricultural production and vehicular transport were the major factors in the concentrations of Cd, Zn, Pb, and Hg in soil from Shanxi Province [58]. van Bohemen [59] found that traffic emissions are the main source of Pb (fuel combustion), Zn (tire wear), and Cu (brake wear and radiator corrosion). Tetraethyl lead was used as an anti-knock agent in automobile gasoline; Pb would be emitted in the exhaust gas during the combustion process of gasoline, and enriched around the roadsides [49]. Zn was often used as a hardness additive in tires, and Cu was a component of brake pads [60]. High quantities of Cd and Zn have been observed in tire rubber [61]. Cd has been proven to be the primary pollutant of highway soil in Beijing, and has mainly come from motor vehicle exhaust and dust produced by tire wear and braking [56]. Zhang [62] also found high concentrations of Cd in smelter dust and automobile exhaust. In this study, Zn, Pb, Cu, and Cd were significantly enriched in the topsoil, and the soil boreholes were located along the main traffic roads, which indicates that traffic sources are major factors. Consequently, Factor 3 may represent traffic source.

Comparison of Multiple Approaches
The source apportionment results obtained by multiple approaches were basically consistent. They all distinguished Cu, Pb, and Zn versus As, Cr, and Ni as two different sources, but were controversial in relation to Hg and Cd. CA classifies sources based on the degree of correlation between elements, but it was difficult to distinguish the elements with complex correlations. In this study, CA mixed the source of Cd with Cu, Pb, Zn, and As. TEF distinguishes natural and anthropogenic sources according to the accumulation of elements in topsoil, but it cannot subdivide anthropogenic sources and is not applicable to a situation where the subsoil contamination is heavier than the topsoil, such as contamination due to leakage of an underground tank. PCA and PMF are both factor analysis methods. PMF imposes non-negative limits on pollution sources and contributions, and also takes into account the uncertainties of the original data. Therefore, PMF can correct data deviations to make results more reliable. In this study, PCA failed to distinguish Hg and Cd from As, Cr, and Ni. PMF apportioned Hg as an isolated factor and classified Cd with Cu, Pb, and Zn as being from the same source-which agreed exceptionally well with existing research conclusions [10,36]. Overall, the results of CA and TEF were based on quantitative analysis, which could be used as a support for PCA and PMF, and the results of PMF were clearer and more complete compared to PCA.

Conclusions
The HM contamination status in the topsoil was heavier than that of the subsoil. The topsoil samples where Hg was over background levels were distributed evenly in the study area, but those for Cd, Cu, Pb, As, and Zn were mainly concentrated at Points 1, 6, 9, and 12. Three sources and their respective contribution rates for each HM were identified using CA, TEF, PCA, and PMF. The natural source contributed the most in relation to As (63.11%), Cr (61.67%), and Ni (70.80%); the source of fossil fuel combustion and atmospheric deposition contributed the most for Hg (97.0%); and the traffic source contributed the most for Zn (72.80%), Pb (69.75%), Cu (65.36%), and Cd (53.08%). The combined use of multiple methods was proven to be effective for HM source apportionment in soil. The results from different methods could be verified against each other, and the results from PMF were clearer and more complete compared to those from other approaches. The HMs originating from anthropogenic sources need to be remediated or appropriately controlled, but no action is required for the HMs originating from natural sources. It is of great importance to fully consider the sources of pollutants during soil remediation in order to avoid over-remediation.