Risk Assessment of Heavy Metals in Soils from Four Different Industrial Plants in a Medium-Sized City in North China

Laboratory experiments were carried out to analyze 39 soil samples collected from four industrial areas in Xuzhou City using inductively coupled plasma mass spectrometry and atomic fluorescence spectrometry. The descriptive statistics of heavy metals (HMs) in the soil profiles showed that the HM content at three depths was highly variable, and most coefficients of variation (CVs) showed moderate variability. The enrichment of Cd at all depths exceeded the risk screening value, and Cd pollution occurred in four plants. The enrichment of the other HMs at three depths was mainly concentrated in the pharmaceutical plant A and chemical plant C. It was found that the different HMs had different vertical distribution characteristics. For the different industrial plants, the raw materials and products not only made the spatial distribution characteristics of the HMs different, but also caused the HM types and contents to differ. The average single pollution indices of Cd in plant A, iron-steel plant B, and plant C indicated a slight pollution level. The other seven HMs in A, B, and C and all HMs in chemical plant D belonged to the safe category. The mean values of the Nemerow pollution index in the four industrial plants belonged to the warning category. The analysis showed that none of the HMs posed potential noncarcinogenic health risks, and only the carcinogenic health risks of Cr in plants A and C were unacceptable. The carcinogenic effect of Cr through the inhalation intake of resuspended soil particulates and that of Cd, Ni, and As via direct oral ingestion were the main exposure pathways.


Introduction
Generally, heavy metals (HMs) refer to copper (Cu), lead (Pb), zinc (Zn), nickel (Ni), chromium (Cr), cadmium (Cd), and mercury (Hg). Metalloid arsenic (As) also belongs to the HM category, due to its similarity of chemical properties and environmental behaviors [1]. With rapid urbanization and industrialization, the soil environment has been polluted, especially by HMs [2,3]. The soil pollution induced by HMs is considered a globally challenging issue [4].
Soil HM pollution has been widespread and serious across China since the economic reform and opening in the late 1970s [5,6]. Some statistics show that over 10 million hectares of land in China are threatened by soil HM pollution [7]. According to an extensive survey on the environmental quality of soils from 2005 to 2013 by the Ministry of Environmental Protection and the Ministry of Land and Resources [8], it was shown that 16.1% of sampling sites exceeded the environmental quality standards in agricultural soils, and HMs accounted for 82.4% of contaminated samples [9]. Statistically, 7.0%, 1.6%, 2.7%, 2.1%, 1.5%, plants. Studies on the spatial distribution characteristics of HMs in soils and the assessment of their environmental health risks are helpful for understanding the sources of soil HM pollution and making decisions regarding soil remediation in industrial plants.
Most studies related to industrial soil HM pollution were carried out in large cities, and less information is available for the medium-sized cities in China. At present, the mediumsized cities in China, which have relatively poor environmental management systems, are facing increased environmental problems, due to rapid urbanization and industrialization. The environmental problems induced in other medium-sized cities around the world are similar to those in China.
Xuzhou is a typical medium-sized city in North China, with a gross domestic product of 675.523 billion RMB in 2018 and a permanent population of 8.802 million at the end of 2018. To our knowledge, there has been little information published on the pollution level, spatial distribution, and environment health risk of soil HMs in the industrial plants of Xuzhou. Consequently, it is of great significance to investigate the pollution level, identify the fundamental sources, and assess the potential health risks of soil HMs in industrial plants. This study can be used as a reference by other similar cities around the world.
In this study, based on on-site soil samplings and lab experiments in four industrial plants in Xuzhou, the objectives were (1) to investigate the pollution level of HMs using the single factor contaminant index and the Nemerow pollution index, (2) to present the spatial distribution characteristics of HMs, and (3) to quantify the environmental health risk utilizing a USEPA recommended model. The presented results will be helpful for both environmental administration officers and factory managers to prevent soil HM pollution by industrial plants in similar medium-sized cities.

Study Area
The study area was located in Xuzhou (33 • 43 01"-34 • 58 34" N, 116 • 21 29"-118 • 40 01" E), Jiangsu Province, with alternating plains and mountains. The plain area accounts for 90%, with the elevation being 20-50 m, while the mountain area accounts for 10%, with the elevation generally being between 100 and 300 m. Xuzhou has a warm temperate monsoon climate, with an annual sunshine duration of 2284-2495 h, an annual average temperature of 14 • C, an annual frost-free period of 200-220 days, and annual average precipitation of 800-930 mm. The climate is characterized by four distinct seasons, abundant sunshine, moderate rainfall, and similar periods of rain and heat. The basic parameters of the soil in the study area are presented in Table 1. The digital elevation and geomorphological information of the study area are shown in Figure 1.  [41,42]. pH: Soil-determination of pH-potentiometry [43]. Moisture content was determined on-site using a three-parameter soil speed tester. Instrument model: Delta-T ML3-KIT [44].
One pharmaceutical plant (A), one iron-steel plant (B), and two chemical plants (C and D) were selected in the study area in Jiawang County, Xuzhou, representing different types of industry. The main soil type in the study area is cinnamon soil.
The annual production capacity of the pharmaceutical plant (A) is 199.48 tons. The main raw and auxiliary materials include ether, ethanol, toluene, dichloromethane, methanol, chloroform, hydrochloric acid, sodium hydroxide, and so on. The iron-steel plant (B) has an annual production capacity of 600,000 tons of continuous casting slab and 900,000 tons of hot-rolled wire rod. The chemical plant (C) has an annual production capacity of 80,000 tons of rutile titanium dioxide, 250,000 tons of by-product ferrous sulfate, and 300,000 tons of sulfur to produce sulfuric acid. The main raw and auxiliary materials include ilmenite, iron powder, caustic soda, hydrochloric acid, limestone powder, sulfur, and so on. The chemical plant (D) is one of the largest producers of vat dyes in China. It has 8 categories and 31 products. The main products are vat dyes and the designed capacity is 9550 tons per year. The main raw and auxiliary materials include anthraquinone, insurance powder, soda ash, sodium acetate, 3,9-dibromobenzoanthrone, benzanthrone, sulfuric acid, hydrochloric acid, chlorobenzene, triethylene glycol, toluene, 1,2-dichloroethane, and lignin. One pharmaceutical plant (A), one iron-steel plant (B), and two chemical and D) were selected in the study area in Jiawang County, Xuzhou, representin types of industry. The main soil type in the study area is cinnamon soil.
The annual production capacity of the pharmaceutical plant (A) is 199.48 main raw and auxiliary materials include ether, ethanol, toluene, dichlorometh anol, chloroform, hydrochloric acid, sodium hydroxide, and so on. The iron-stee has an annual production capacity of 600,000 tons of continuous casting slab an tons of hot-rolled wire rod. The chemical plant (C) has an annual production c 80,000 tons of rutile titanium dioxide, 250,000 tons of by-product ferrous su 300,000 tons of sulfur to produce sulfuric acid. The main raw and auxiliary ma clude ilmenite, iron powder, caustic soda, hydrochloric acid, limestone powd and so on. The chemical plant (D) is one of the largest producers of vat dyes in has 8 categories and 31 products. The main products are vat dyes and the design ity is 9550 tons per year. The main raw and auxiliary materials include anthr insurance powder, soda ash, sodium acetate, 3,9-dibromobenzoanthrone, benz sulfuric acid, hydrochloric acid, chlorobenzene, triethylene glycol, toluene, 1,2 ethane, and lignin.

The Sampling Layout Scheme
Based on previous data collection and an on-site survey, and according to the Technical Guidelines for Environmental Site Investigation, Technical Guidelines for Environmental Site Monitoring [10], and Technical Guidelines of Soil and Groundwater Self-Monitoring for Enterprises In Production [45], the sample points were arranged in key areas, such as the production area, sewage treatment area, hazardous chemical warehouse, calcination area, and warehouse. There were 10 soil sampling points in plants A, C, and D, and 9 soil sampling points in plant B. The sampling points at the four plants are shown in Figure 2.
Soil samples were collected at the target depth using a Geoprobe drilling rig (6610DT) into a PE sampling tube with a length of approximately 20-30 cm. The two ends of the sampling tube were sealed and labeled. Soil samples from the A, B, and C plants were

Laboratory Testing and Analysis
After sample collection, samples were immediately stored in a 0-4 • C refrigerator and sent to the laboratory within 48 h. Laboratory testing and analysis was carried out at the laboratory of the Jiangsu Fangzheng Environmental Protection Group Co., Ltd. To reduce the impacts of the process of sample preservation and transportation, 2-3 cm soil at the bottom of each sampling tube was removed and 10-60 g soil at the bottom of each sample was taken for testing and analysis. According to "Soil and Sediment-Determination of Aqua Regia Extracts of 12 Metal Elements-Inductively Coupled Plasma Mass Spectrometry" [46], the inductively coupled plasma mass spectrometry method with aqua regia treatment was used for the determination of Cu, Pb, Zn, Ni, Cr, Cd, and As. According to "Soil Quality-Analysis of Total Mercury, Arsenic and Lead Contents-Atomic Fluorescence Spectrometry" [47], the atomic fluorescence spectrometry method was used for the determination of Hg. All metal data mentioned in this study refer to total contents. At least two laboratory blank samples were prepared for each batch of samples, to ensure that the determination was accurate. The correlation coefficient of the standard curve established for each analysis was greater than 0.999. Each batch of samples was tested with parallel double samples and spiked recovery samples at a proportion of at least 10%. Under constant conditions, the relative deviation of the two independent determination results obtained was no more than 12%. Information regarding accuracy is presented in Table 2.

Pollution Level Evaluation Method
Both the single-factor contaminant index P i and the Nemerow pollution index P N were used to evaluate the HM pollution level in soil [48,49]. P i and P N were calculated using Equations (1) and (2), respectively.
where P i (unitless) is the single-factor contaminant index, C i (mg·kg −1 ) is the detected concentration of each HM i, and S i (mg·kg −1 ) is the evaluation criterion of each HM i.
The evaluation criterion of this study was the risk screening value (RSV) for soil contamination of agricultural land [11]. The P i value was classified as safety (P i ≤ 1), slight pollution (1 < P i ≤ 2), low pollution (2 < P i ≤ 3), moderate pollution (3 < P i ≤ 5), and high pollution (P i > 5).
where P N (unitless) is the Nemerow pollution index, P imax is the maximum value of P i , and P iave is the arithmetic mean value of P i . The P N value was classified as safety (P N ≤ 0.7), warning (0.7 < P N ≤ 1), slight pollution (1 < P N ≤ 2), moderate pollution (2 < P N ≤ 3), and heavy pollution (P N > 3) [50].

Exposure Assessment of HMs
Exposure assessment of HMs can be performed through direct oral ingestion of the soil (OISER j ), dermal absorption via soils adhered to exposed skin (DCSER j ), and the inhalation intake of resuspended soil particulates (PISER j ) [1]. Three types of exposure dose were calculated using Equations (4)-(6), respectively [51].
where the subscript j represents a noncarcinogenic (nc) or carcinogenic (ca) effect, the subscript a represents adults, OSIR a (mg·d −1 ) is the ingestion rate, ED a (year) is the exposure duration, EF a (d·year −1 ) is the exposure frequency, ABS o (unitless) is the oral absorption factor, BW a (kg) is the body weight, AT j (d) is the average exposure time for adults in their whole lifetime, SAE a (cm 2 ) is the exposed skin surface area, SSAR a (mg·cm −2 ) is the adherence rate of soil on skin, E v (unitless) is the daily exposure frequency of the dermal contact event, ABS d (unitless) is the dermal absorption factor, PM 10 (mg·m −3 ) is the content of inhalable particulates in ambient air, DAIR a (m 3 ·days −1 ) is the daily air inhalation rate of adults, PIAF (unitless) is the retention fraction of inhaled particulates in the body, f spo (unitless) is the fraction of soil-borne particulates in outdoor air, f spi (unitless) is the fraction of soil-borne particulates in indoor air, EFO a (d·a −1 ) is the outdoor exposure frequency, and EFI a (d·a −1 ) is the indoor exposure frequency.

Health Risk Assessment
The noncarcinogenic effect risk was assessed using the hazard quotient (HQ, including HQ ois , HQ dcs , and HQ pis ), describing three noncarcinogenic risk indices. They were calculated using Equations (7)-(11) [52].
where HQ ois , HQ dcs , and HQ pis are hazard quotients via the three pathways, respectively. C sur (mg·kg −1 ) is the mean concentration of each HM in the soil samples; SAF (unitless) is the soil allocation factor; RfD o , RfD d , and RfD i (mg·kg −1 ·d −1 ) are the corresponding reference doses via different pathways; ABS gi (unitless) is the digestive tract absorption efficiency factor, and RfC (mg·m −3 ) is the respiratory inhalation reference concentration. The hazard index (HI) is equal to the sum of HQ and was used to estimate the total noncarcinogenic effects created by all HMs and exposure pathways. If the HI value is more than 1, there is a chance that an adverse noncarcinogenic health effect will occur [53].
CR ois , CR dcs , and CR pis describing the carcinogenic risk indices via three pathways were calculated by Equations (12)- (16).
where SF o , SF d , and SF i (mg·kg −1 ·d −1 ) are the carcinogenic slope factor via the three pathways, and IUR (m 3 ·mg −1 ) is the breathing inhalation carcinogen unit. For a single metal, this is negligible when CR is less than 10 −6 , while it is unacceptable when CR is more than 10 −4 [54]. For all the carcinogenic HMs or exposure pathways, the acceptable level of accumulated carcinogenic risk is less than 10 −5 [55]. Since the health effect mechanism of Pb is different from the others, the health risk for Pb was not considered in this study. Table 3 shows the parameters of the health risk assessment model in this study. Table 3. Parameters of the human health risk assessment model in this study.

Statistical Analysis
Descriptive statistical indices including the mean, sample variance (SV), sample standard deviation (SSD), coefficient of variation (CV), kurtosis, and skewness were used. CV ≤ 20%, 20% < CV ≤ 50%, 50% < CV ≤ 100%, and CV > 100% indicated low, moderate, high, and exceptionally high variability, respectively [50]. Pearson coefficient correlation analysis was used to investigate the potential pollution sources (industrial activities) of HMs. If there is a significant positive correlation between HMs in soils, they have a high possibility of possessing similar sources. Table 4 shows the descriptive statistical characteristics of HMs in the soil profiles. The HM content at the three depths was highly variable and most CVs were more than 20% and less than 50%, showing moderate variability, except for Hg at a 1 m depth and Pb, Ni, and As at 5 m depth. This can probably be attributed to the different HM background contents and industrial activities.  [41,42]. RSV: Risk screening value for soil contamination of agricultural land [11]. RIV: Risk intervention value for soil contamination of agricultural land [11]. Soil environmental quality: Risk control standard for soil contamination of agricultural land [11]. AVSP: The average value of the soil profile.

The Descriptive Statistics of HMs in Soil Profiles
The mean values of HM contents at different depths were as follows: At 1 m depth, Cu, Pb, Zn, Ni, Cd, and As were 1. From the average value of the soil profile (AVSP), the content of HMs was highly variable as well. With the exception of Zn, Cr, Hg, and Cd, the AVSPs of the other four HMs were more than the BG and less than the RSV. The AVSPs of Cu, Pb, Ni, Cd, and As were 1.18, 1.42, 1.19, 3.80, and 1.13 times the BG. The maximum values of Pb, Ni, and Cd all exceeded the RSV, and the maximum value of Cd was 6.93 times the RSV. The AVSP of Cd was more than the RSV in the soil profile. This indicates that Cd exceeded its RSV value and may pose a risk to human health.
The order of the CVs of HMs was Hg (103.26%) > Pb (59.65%) > Cd (59.37%) > Zn (56.40%) > Cr (54.02%) > As (50.51%) > Ni (48.02%) > Cu (42.23%). While Hg showed exceptionally high variability, the other seven HMs showed moderate to high variability. This might be due to the influence of different degree of chemical processes on the soil HM pollution in the studied area. It can be seen that Hg, Pb, Cd, and Ni had a positive skewness, with a large value in Table 4. Similar studies showed that HMs have positive skewness if they are disturbed by human activities.

The Descriptive Statistics of HMs in the Different Industrial Plants
In Table 5, the mean values of Cu, Pb, Zn, Ni, Cr, Cd, and As for all depths are higher than the BV in A plant; As and Cd for all depths were higher than the BV in B plant; Cu, Pb, Zn, Ni, Cd, and As for all depths were higher than the BV in C plant; and Cd and Hg for all depths were higher than the BV in D plant. This indicates that many HMs were enriched to some extent in the four industrial plants. Cd in A, B, and C and at 4.5-6 m depth in D exceeded the RSV, indicating that Cd pollution occurred in all four industrial plants. Generally, the enrichment of soil HMs at different depths was mainly concentrated in the A and C industrial plants. This indicates that A and C probably cause greater HM pollution than B and D. The Hg at 5 m depth in the A plant, and Zn at 5 m depth, Ni at 1 m depth, and Cr at 1 m depth in the C plant, showed low variability. In comparison, Cd at 2.5-3 m depth, and Zn and As at 4.5-6 m depth, in the D plant, showed high variability. The rest showed moderate variability. Generally, most HMs showed a moderate variability in all the industrial plants in this study.   Figure 3 and Figure S1 show the vertical distribution of the content of the eight HMs in the different industrial plants. Comparisons were made between different layers for the same HM. It should be noted that the 5 m depth layer at the S08 site in industrial plant C was polluted by Pb and Ni through industrial activities. The content of Pb and Ni at a 5 m depth was higher than that of the upper two layers, indicating that Pb and Ni had migrated from the surface to the 5 m depth layer. HM pollution in deep soil may endanger the quality of the groundwater environment. The content of Hg at 1 m depth at S09 in B and at S03 in D was higher than the BV of Hg (0.04 mg·kg −1 ). Hg was introduced into certain sites through industrial activities in the B and D industrial plants.

Vertical Distribution Characteristics of Soil HMs in the Four Industrial Plants
same HM. It should be noted that the 5 m depth layer at the S08 site in industrial plant C was polluted by Pb and Ni through industrial activities. The content of Pb and Ni at a 5 m depth was higher than that of the upper two layers, indicating that Pb and Ni had migrated from the surface to the 5 m depth layer. HM pollution in deep soil may endanger the quality of the groundwater environment. The content of Hg at 1 m depth at S09 in B and at S03 in D was higher than the BV of Hg (0.04 mg·kg −1 ). Hg was introduced into certain sites through industrial activities in the B and D industrial plants.  The content of Cd at most sites of the four industrial plants exceeded not only the BV (0.1 mg·kg −1 ) but also the RSV (0.3 mg·kg −1 ), including S01, S02, S04, and S06-S08 in A, S01-S09 in B, S01-S04 and S06-S10 in C, and S07-S09 in D, which can be seen in Table 5. The content of Cd at the S02 site in the D region even exceeded 2.0 mg·kg −1 . Cd-induced pollution occurred in all industrial plants.

Spatial Distribution of HMs in Four Industrial Plants
It can be seen that Cu and Hg increased from north to south, and the content was the highest at the S04 site located between the No. 5 and No. 6 synthesis workshops, as shown in Figure 3. This high content of Cu and Hg might be related to industrial activities. Generally, with the exception of Ni and As, the No. 5 and No. 6 synthesis workshops and the dangerous goods warehouse had a strong effect on the content of the HMs.
With the exception of Hg and Pb, the content of the other six HMs in the west was relatively higher than that in the east, as shown in Figure S2(A). All HMs had a high content around the S07 and S08 sites in the southeast, whose sites were located near the calcination section in the C plant, as shown in Figure S2(B). High contents of Cu, Zn, and Ni were concentrated in the S01 site, surrounding the chemical laboratory building, while the high content of Hg was concentrated in S03 and S04, surrounding the dangerous goods warehouse and sewage treatment station, as seen in Figure S2(C).
The sequence of the average content of Cu and Ni was C > A > B > D. The sequence of the average content of Pb, Cr, and Zn was A > C > B > D. The sequence of the average content of Hg was D > A > B > C, while that of As was A > B > C > D. As for Cd, the proportion of sites exceeding the RSV was as high as 79.5% overall.
Three possible reasons may explain the different spatial distribution characteristics of the different HMs. The first reason is that the raw materials and products in different industries were different, causing the types and content of HMs to differ. The second reason might be that the contamination times of soils in the different industrial plants varied. Different HM elements have different migration and adsorption characteristics. The content of Cd at most sites of the four industrial plants exceeded not only the BV (0.1 mg·kg −1 ) but also the RSV (0.3 mg·kg −1 ), including S01, S02, S04, and S06-S08 in A, S01-S09 in B, S01-S04 and S06-S10 in C, and S07-S09 in D, which can be seen in Table 5. The content of Cd at the S02 site in the D region even exceeded 2.0 mg·kg −1 . Cd-induced pollution occurred in all industrial plants.

Spatial Distribution of HMs in Four Industrial Plants
It can be seen that Cu and Hg increased from north to south, and the content was the highest at the S04 site located between the No. 5 and No. 6 synthesis workshops, as shown in Figure 3. This high content of Cu and Hg might be related to industrial activities. Generally, with the exception of Ni and As, the No. 5 and No. 6 synthesis workshops and the dangerous goods warehouse had a strong effect on the content of the HMs.
With the exception of Hg and Pb, the content of the other six HMs in the west was relatively higher than that in the east, as shown in Figure S2A. All HMs had a high content around the S07 and S08 sites in the southeast, whose sites were located near the calcination section in the C plant, as shown in Figure S2B. High contents of Cu, Zn, and Ni were concentrated in the S01 site, surrounding the chemical laboratory building, while the high content of Hg was concentrated in S03 and S04, surrounding the dangerous goods warehouse and sewage treatment station, as seen in Figure S2C.
The sequence of the average content of Cu and Ni was C > A > B > D. The sequence of the average content of Pb, Cr, and Zn was A > C > B > D. The sequence of the average content of Hg was D > A > B > C, while that of As was A > B > C > D. As for Cd, the proportion of sites exceeding the RSV was as high as 79.5% overall.
Three possible reasons may explain the different spatial distribution characteristics of the different HMs. The first reason is that the raw materials and products in different industries were different, causing the types and content of HMs to differ. The second reason might be that the contamination times of soils in the different industrial plants varied. Different HM elements have different migration and adsorption characteristics. Table 6 shows the Pearson coefficient correlation analysis of HMs in the soils. The correlation between HMs could effectively indicate the sources of HMs. There were significant positive correlations between Cu and Pb, Ni, Cr, and As, with correlation coefficients of 0.364, 0.523, 0.577, and 0.418, respectively. The correlation coefficients between Pb and Ni, Cr, Cd, and As were 0.685, 0.422, 0.416, and 0.323, respectively. The correlation coefficients of Ni with Cr and As were 0.656 and 0.556, respectively. The correlation coefficient between Cr and As was 0.546. From the above, the correlation coefficients of Cu-Ni (0.523), Cu-Cr (0.577), Pb-Ni (0.685), Ni-Cr (0.656), and Ni-As (0.556) were relatively high, indicating that they might have the same sources. Cd was only positively correlated with Pb, and the correlation coefficient was only 0.416, indicating that the Cd might have a single pollution source.

Single Pollution Index of Surface Soil in the Different Industrial Plants
The average single pollution indices of Cd in the A, B, and C industrial plants were more than 1 and less than 2, which indicated a low pollution level, while the other seven HMs belonged to the safe category, as seen in Table S1. Most CVs of the single pollution index of HMs in the A and B industrial plants were more than 20% and less than 50%, which showed moderate variability. Only the CVs of the single pollution indices of Zn, Ni, and Cr in industrial plant C were less than 20%, showing low variability, and those of the other five HMs showed moderate variability. For industrial plant D, the average single pollution indices of all HMs were less than 1 and at a safe level, indicating that the soil in industrial plant D was generally good. The average values of the single pollution indices of Cu, Pb, Zn, Ni, Cr, Cd, Hg, and As in the four plants were 0.31, 0.21, 0.30, 0.34, 0.26, 1.11, 0.02, and 0.44, respectively (see Table S1). There were 28 slightly contaminated sites of Cd, accounting for 71.79%.

The Nemerow Pollution Index of Surface Soil in the Different Industrial Plants
In Table 7, the mean values of the Nemerow pollution index in the A, B, C, and D industrial plants were 0.86, 0.90, 0.94, and 0.70, respectively, which were all more than 0.7 and less than 1.0, corresponding to the warning level. All maximum values in the four industrial plants were more than 1.0, indicating that some sites in each plant had a slight pollution level. The CVs of the Nemerow pollution index in the A, B, C, and D industrial plants were 38.01%, 16.32%, 21.09%, and 48.17%, respectively, with low to moderate variability. The mean value of the Nemerow pollution index of all industrial plants was 0.85, showing a warning level, with the maximum Nemerow pollution index and CV being 1.40 and 32.22%, respectively. As seen in Table 8, there were 10, 19, and 10 sites with a safe level, warning level, and slight pollution level, accounting for 25.64%, 48.72%, and 25.64% of the 39 sites in the four industrial plants.  Table S2 shows the HQs and HIs of soil HMs under different exposure pathways at the different industrial plants. The HQs of Cr and As were higher than those of the other HMs at all industrial plants. As seen in Table 9, the health risk of As via the direct oral ingestion and inhalation intake of resuspended soil particulates reached 51.75% and 47.95%, respectively, and that of Cr through inhalation intake of resuspended soil particulates reached 53.08%. However, the average HQ ois , HQ dcs , and HQ pis were all less than 1, indicating that the Cr and As in all industrial plants did not pose potential noncarcinogenic health risks to humans. In addition, the noncarcinogenic risks of Cu, Zn, and Hg were mainly through direct oral ingestion, while those of Cd and Ni were mainly through the inhalation intake of resuspended soil particulates. The average HI s of soil HMs for the four industrial plants were all less than 1. The order of the HI values of soil HMs was As > Cr > Ni > Cd > Cu > Zn > Hg. Table 10 shows the health risk assessment of the carcinogenic effects of soil HMs under different exposure pathways at the different industrial plants. The CRs of Cr in plants A and C were unacceptable (Table 10). The average CRs of Cr, As, Cd, and Ni in the four industrial plants were 9.69 × 10 −5 , 9.48 × 10 −6 , 2.59 × 10 −6 , and 1.37 × 10 −5 , respectively, with an order of Cr > Ni > As > Cd. The carcinogenic effect of Cr through the inhalation intake of resuspended soil particulates reached 86.37%, and the carcinogenic effect of Cd, Ni, and As via direct oral ingestion reached 81.06%, 86.44%, and 87.83%, respectively (Table 7).

Conclusions
This study showed that the HMs at all depths were highly variable, and most CVs showed moderate variability. Cd at three depths even exceeded the RSV. Hg, Pb, Cd, and Ni might have been disturbed by industrial activities. HMs were enriched at all four industrial plants. Cd pollution occurred at all industrial plants. Different HMs had different vertical distribution characteristics. The raw materials and products not only made the soil HM spatial distribution characteristics different but caused the metal types and contents to be different. The correlation analysis of soil HMs suggested that Cu-Ni, Cu-Cr, Pb-Ni, Ni-Cr, and Ni-As might have the same sources. Cd might have a single pollution source. More attention should be paid to the raw materials during the production process to minimize HM pollution. The average single pollution indices of Cd at the A, B, and C industrial plants implied a slight pollution level, which indicated that Cd posed a possible risk to human health. The other seven HMs belonged to the safe category. The HQ s of Cr and As were higher than the those for other HMs and were all less than 1, indicating that Cr and As did not pose potential noncarcinogenic health risks to humans at all the industrial plants for the area investigated in this study.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/toxics11030217/s1, Figure S1: Vertical distribution of eight HM contents at different depths in different industrial plants; Figure S2: Spatial distribution characteristics of HMs in four industrial plants; Table S1: The single factor contaminant index of HMs in different industrial plants; Table S2 Conflicts of Interest: All authors certify that they have no affiliations with or involvement in any organization or entity with any financial interest or non-financial interest in the subject matter or materials discussed in this manuscript.