Spatial Distribution, Environmental Risk and Safe Utilization Zoning of Soil Heavy Metals in Farmland, Subtropical China

: Heavy metal (HM) accumulation in farmland soil can be transferred to the human body through the food chain, posing a serious threat to human health. Exploring the environmental risk and safe utilization zoning of soil HMs in farmland can provide the basis for the formulation of effective control strategies. Soil samples from typical subtropical farmland were collected in Jinhua City and analyzed for HMs (As, Cd, Cr, Cu, Hg, Ni, Pb, and Zn). The objective of this study was to explore the spatial distribution and environmental risk of soil HMs, and then divide the safe utilization area of soil HMs of farmland in Jinhua City. The results showed that the mean concentrations of soil HMs were, in descending order: Zn (76.05 mg kg − 1 ) > Cr (36.73 mg kg − 1 ) > Pb (32.48 mg kg − 1 ) > Cu (18.60 mg kg − 1 ) > Ni (11.95 mg kg − 1 ) > As (6.37 mg kg − 1 ) > Cd (0.18 mg kg − 1 ) > Hg (0.11 mg kg − 1 ), and all determined soil HMs did not exceed the risk screening values for soil contamination of agricultural land of China. The ﬁtted semi-variogram showed that the spatial autocorrelation of Cd, Hg, Pb, and Zn was weak, with island-shaped distribution, while As, Cr, Cu, and Ni had medium spatial autocorrelation, with strip-shaped and island-shaped distribution. The hot spot analysis and environmental risk probability showed that the environmental risks of As, Cd, Cu, Pb, Zn, and Cu were relatively high, whereas those of Cr, Hg, and Ni were relatively low. Safe utilization zones and basic safe utilization zones accounted for 89.35% and 8.58% of the total farmland area in Jinhua, respectively, and only a small part of the farmland soil was at risk of use.


Introduction
Soil is one of the most important natural resources for human survival and development, and it is also the source and sink of many pollutants [1,2]. With rapid industrialization and urbanization, heavy metals (HMs) have gradually accumulated in farmland soil [3,4], which have not only degraded the soil quality of farmland, but can also affect the growth of plants [5,6]. Furthermore, soil HMs can be transferred to the human body through the food chain, posing a serious threat to human health [7,8].
Recently, GIS-based geostatistical analysis methods have been widely used in the study of soil spatial distribution [9][10][11]. Among them, ordinary kriging quantifies the spatial distribution pattern by fitting a semi-variogram, which is one of the most commonly used methods to explore the spatial distribution of soil HMs [12]. Additionally, indicator kriging, which is a nonparametric approach, is more advanced due to its ability to take the Land 2021, 10, 569 2 of 13 data uncertainty into account and it is used to predict the conditional probability for the categorical data of unsampled sites [13,14]. Many scholars applied the indicator kriging method to predict the spatial distribution and pollution risk probability of soil HMs with the soil background value as the threshold [15][16][17]. The Getis-Ord spatial statistics index (G i *), also called hot spot analysis, is calculated to identify statistically significant spatial clusters of high values (hot spots) and low values (cold spots) [18]. The hot spots of soil HMs and their surroundings may have environmental risks and, consequently, the G i * index can be used for the identification of the environmental risk of soil HMs [12,19]. Currently, the safe utilization zoning of HMs in farmland has mainly been based on the assessment of soil HMs pollution, and a unified standard system has not yet been formed. Many assessment methods were proposed for measuring soil HMs pollution, such as the single pollution index [20], Nemerow index [21], geo-accumulation index [22], pollution load index [23], potential ecological risk index [24], enrichment factor index [25], etc. However, these methods are either not able to comprehensively reflect HMs pollution in soil, or may be subjective. The influence index of comprehensive quality of soil (IICQs) comprehensively considers the valence effect of elements, the environmental quality standard for soil, the background value of soil elements, and the specific load capacity of soil, which could scientifically and comprehensively evaluate HMs pollution in soil [26]. Therefore, this research adopts the IICQs to evaluate HMs pollution in farmland soil and carry out safe utilization zoning of soil HMs.
Jinhua has a long history of agricultural cultivation and is one of the most important granaries in Zhejiang Province. Therefore, exploring the pollution of HMs in farmland soil in Jinhua is essential to ensure food safety and human health. Previous studies on heavy metals in Jinhua mainly focused on urban soil, and little attention was paid to the characteristics and risk of heavy metals in farmland soil [27]. We hypothesized that the combined multiple methods would effectively assess the risk status of soil HMs and accurately delineate the safe utilization zoning. Specifically, the main objectives of this research were to (1) explore the spatial variability and distribution characteristics of HMs in farmland soils in the study area; (2) identify the environmental risks of HMs in farmland soils in the study area; (3) carry out safe utilization zoning of HMs in farmland soil based on the IICQs.

Study Area
Jinhua City (28 • 32 ~29 • 41 N, 119 • 14 ~120 • 47 E) is located in the middle of Zhejiang Province, eastern China, with an area of about 1.09 × 10 4 km 2 ( Figure 1). The study area has a subtropical monsoon climate, warm and humid, with abundant rainfall and four distinct seasons. The annual average temperature and precipitation are about 16.3~17.6 • C and 1150~1909 mm, respectively. The vegetation is subtropical evergreen broad-leaved forests. The area of farmland in Jinhua is about 2400 km 2 , and is dominated by paddy fields and dry land. Jinhua is also an important grain production area in Zhejiang Province, with annual grain output of 4.35 × 10 5 tons. In recent years, the concentrations of heavy metals in farmland soils in Jinhua have gradually increased with the massive use of pesticides and industrial emissions.

Soil sampling and Laboratory Analyses
According to the principle of typicality, representativeness, and consistency, soil samples were collected by using a 5 km × 5 km grid, then the non-farm sampling points were removed after field investigations, and a total of 209 farmland topsoil samples (0-20 cm) were finally collected ( Figure 1). For each sample, five subsamples within a 100 m radius from the sampling site were collected and mixed into a composite sample, stored in polyethylene bags, and taken back to the laboratory for further analysis.
All soil samples were air-dried and ground in a dust-free laboratory, and after that passed through a 2 mm nylon sieve to remove stones and other debris, and then passed through a 0.149 mm sieve for testing. Soil pH was measured by a pH meter with a ratio of 1:2.5 soil to deionized water [28]. Soil Cr, Cd, Cu, Pb, Ni, and Zn were digested with HCl-HNO3-HClO4 on a hot plate, and then measured by inductively coupled plasma-atomic emission spectrometry (ICP-AES); soil As and Hg were digested with aqua regia water bath digestion, and then measured by an inductively coupled plasma-atomic fluorescence spectrometer (ICP-AFS). All reagents used in the analyses were of premium grade, and all water was ultrapure water. The quality control was carried out with the national soil standard reference material (GBW07424, GSS-24) provided by the Center of National Standard Reference Material of China [29]. Additionally, 20% of the samples were randomly selected for testing repeatedly, and the error remained within 5%. The recovery rate ranged from 81% to 109%, and the measurement results were within the allowable range of error.

Ordinary Kriging
Ordinary kriging is an unbiased optimal estimation of soil HMs concentration by fitting a semi-variogram, which is used to quantify the structural characteristics and random changes of regionalized variables [11,24,30]. The calculation formula of the semi-variogram is as follows:

Soil Sampling and Laboratory Analyses
According to the principle of typicality, representativeness, and consistency, soil samples were collected by using a 5 km × 5 km grid, then the non-farm sampling points were removed after field investigations, and a total of 209 farmland topsoil samples (0-20 cm) were finally collected ( Figure 1). For each sample, five subsamples within a 100 m radius from the sampling site were collected and mixed into a composite sample, stored in polyethylene bags, and taken back to the laboratory for further analysis.
All soil samples were air-dried and ground in a dust-free laboratory, and after that passed through a 2 mm nylon sieve to remove stones and other debris, and then passed through a 0.149 mm sieve for testing. Soil pH was measured by a pH meter with a ratio of 1:2.5 soil to deionized water [28]. Soil Cr, Cd, Cu, Pb, Ni, and Zn were digested with HCl-HNO 3 -HClO 4 on a hot plate, and then measured by inductively coupled plasma-atomic emission spectrometry (ICP-AES); soil As and Hg were digested with aqua regia water bath digestion, and then measured by an inductively coupled plasma-atomic fluorescence spectrometer (ICP-AFS). All reagents used in the analyses were of premium grade, and all water was ultrapure water. The quality control was carried out with the national soil standard reference material (GBW07424, GSS-24) provided by the Center of National Standard Reference Material of China [29]. Additionally, 20% of the samples were randomly selected for testing repeatedly, and the error remained within 5%. The recovery rate ranged from 81% to 109%, and the measurement results were within the allowable range of error.

Ordinary Kriging
Ordinary kriging is an unbiased optimal estimation of soil HMs concentration by fitting a semi-variogram, which is used to quantify the structural characteristics and random changes of regionalized variables [11,24,30]. The calculation formula of the semivariogram is as follows: 2 (1) where γ(h) is the semi-variogram; h is the space separation distance between the two samples; Z(x i ) is the sample value of Z(x) at the spatial position x i ; Z(x i + h) is the sample value of Z(x) at the spatial position x i + h; N(h) is the number of sample point pairs.

Indicator Kriging
Indicator kriging is a non-linear kriging method for estimating the probability of exceeding a specific threshold value at a given location, thereby indicating environmental risks [17]. The key point of indicator kriging is to perform binary encoding of the original data according to the threshold (Z k ) [31], and the relationship is measured as follows: then, the semi-variogram was used to quantify the spatial correlation of the indicator codes, and calculate the probability of the indicator variable at the unsampled points through the following linear equation: where I (x α ; Z k ) is the indicated transformation value of the known sample point x α ; I*(x 0 ; Z k ) represents the estimated value of x 0 to be predicted; n is the number of samples, α = 1, 2, . . . , n; λ α is the weight coefficient, which can be solved by the following equations.
where η is the Lagrange multiplier; γ I is the indicated semi-variance value at the respective lag distances when Z k is the threshold, α = 1, 2, . . . , n.

Hot Spot Analysis
The Getis-Ord spatial statistics index (G * i ) is one of the local spatial autocorrelation analysis methods, which can identify statistically significant spatial clusters of hot spots and cold spots by calculating G * i statistics [18]. The hot spots of HMs represent areas with statistically significant spatial clustering of samples with high HMs concentrations, and the cold spots are areas in which the lowest HMs concentrations are spatially clustered [32]. Hot spots of HMs in the soil may have potential pollution risks. The formula of the G * i index is as follows: where x j represents the attribute value of spatial unit j; w i,j represents the spatial weight between spatial units i and j; n is the number of spatial units; X is the mean value; S is the standard deviation. A significantly positive G * i statistical value indicates hot spots, and the higher the G * i statistical value, the more densely the hot spots are clustered. A significantly negative G * i statistical value indicates cold spots, and the lower the G * i statistical value, the more tightly the cold spots are clustered [33,34]. The IICQs comprehensively consider factors such as the background value of soil HMs, soil quality standards, and valence effects, and can comprehensively reflect the degree of soil pollution [27]. The formula of the IICQs is as follows: where N indicates the number of HMs; C i indicates the measured concentration of HM i ; C Si indicates the soil environmental quality standard value of HM i , which is the soil risk screening value; C Bi indicates the background value of HM i ; n represents the oxidation number of HM i ; RIE represents relative impact equivalent; DDDB represents deviation degree of the determination concentration from the background value; DDSB represents the deviation degree of soil standard from the background value; X is the number of soil-measured values exceeding the soil risk screening values; and Y is the number of soil-measured values exceeding the background values. The pollution levels of soil HMs were classified according to the IICQs and "The National Soil Pollution Condition Investigation Report" [35], and then the safe utilization zoning of soil HMs in farmland was carried out based on the pollution levels (Table 1). Table 1. Classification of soil environmental quality and safe utilization.

Statistical Analysis
The descriptive statistics were produced and the Kolmogorov-Smirnov (K-S) test was conducted using SPSS 19.0 for Windows (SPSS Inc., Chicago, IL, USA). The fitting of the semi-variogram was performed in GS+ 9.0 for Windows (Gamma Inc., Ann Arbor, MI, USA). The kriging interpolation, hot spot analysis, and safe utilization zoning were estimated in ArcGIS 10.2 for Windows (ESRI Inc., Redlands, CA, USA).

Concentration Characteristics of Soil HMs
The descriptive statistics of HMs in farmland soils in Jinhua are summarized in Table 2. Soil pH ranged from 3.99 to 8.15, with an average pH of 5.34, which was mainly acidic. The study area has a subtropical monsoon climate with abundant rain and heat resources, and the strong leaching makes the soil acidic [36]. In addition, the application of chemical fertilizers intensifies the acidification of farmland soils [37]. The mean concentrations of soil HMs were, in descending order: Zn (76.05 mg kg −1 ) > Cr (36.73 mg kg −1 ) > Pb (32.48 mg kg −1 ) > Cu (18.60 mg kg −1 ) > Ni (11.95 mg kg −1 ) > As (6.37 mg kg −1 ) > Cd (0.18 mg kg −1 ) > Hg (0.11 mg kg −1 ), and all determined soil HMs did not exceed the risk screening values for soil contamination of agricultural land of China [38]. However, the Land 2021, 10, 569 6 of 13 mean concentrations of Cu, Hg, and Zn were higher than the corresponding background values of Jinqu Basin and Zhejiang Province, indicating anthropic inputs of these elements were relatively high. The mean concentrations of Cu, Cr, Ni, and Pb in farmland soil in Jinhua City were lower than those in Dehui except for Zn [39]. Similarly, the mean concentrations of As, Cd, Cr, Cu, Hg, Ni, Pb, and Zn in farmland soil in Jinhua City were all lower than those in the agricultural soil south of the Yangtze River Delta [19]. According to the coefficient of variation (CV) of HMs, Pb and Zn had moderate variability (20% < CV < 50%), and As, Cd, Cr, Cu, and Ni had high variability (50% < CV < 100%), whereas Hg had extremely high variability (CV > 100%) [4]. It could be suggested that the distribution of Hg in the study area was extremely uneven and might be dominated by anthropogenic sources [40]. Kurtosis and skewness indicated that the distribution of soil HMs in the study area was a skewed distribution.

Spatial Distribution Characteristics of Soil HM Concentration
The parameters of fitted semi-variogram models are shown in Table 3. The exponential model was fitted for all HMs. The R 2 of both As and Ni was higher than 0.8, and the RSS values were both close to 0, which indicated the fitted models could better reflect the spatial structure characteristics of soil HMs. The R 2 values of Cd, Cr, Cu, Pb, and Zn were higher than 0.5, and RSS values were also close to 0, which indicated the fitted models basically met the requirements. The ratio of nugget (C 0 ) and sill (C 0 + C) was usually used to reflect the spatial heterogeneity in the nugget variance [40]. The C/(C 0 + C) ratios of As, Cr, Cu, and Ni were between 0.25 and 0.75, which indicated a moderate spatial correlation, and the variable was dominated by the combined effects of natural and human factors. The C/(C 0 + C) ratios of Cd, Hg, Pb, and Zn were higher than 0.75, indicating that there was strong spatial randomness and weak autocorrelation [40,43]. The range of spatial autocorrelation of soil HM concentrations was between 8400 and 77,100 m, which was much larger than the sampling interval (5000 m), indicating that the sampling design in this paper was reasonable and could reflect the true variation characteristics of soil HMs [44]. Ni had the largest range (77,100 m), indicating that it was spatially correlated at a larger scale, and the intrinsic factors (such as soil parent material, climate, and topography) had a greater impact on its spatial variability in the soil [45]. However, the ranges of Cd, Hg, Pb, and Zn were smaller, which implied that the spatial variability of these elements was greatly affected by extrinsic factors, such as agricultural activities, industrial emissions, vehicle emissions, etc. [40,45]. Table 3. Parameters of the semi-variogram models for heavy metals in farmland soil.

Element
Model C 0 C 0 + C C/(C 0 + C) Range (m) R 2 RSS The spatial distribution of soil HMs in farmland of Jinhua City is shown in Figure 2. The As, Cr, Cu, and Ni presented strip-shaped and island-shaped distributions. The high concentrations of As were mainly distributed in the central and western parts of the study area, where light industry and agriculture are located, and the population and traffic is dense. In addition, the high concentrations of As were distributed in the east of Pujiang (IV), which was basically consistent with the distribution of farmland. Previous studies have shown that long-term application of phosphate fertilizers and pesticides can increase the As concentration in farmland soils [46,47]. The concentrations of Cr and Ni had similar spatial distributions, which were mainly concentrated in the central and western parts of the study area, indicating that they may have common sources [48]. Moreover, the concentrations of Cr and Ni in the southwest of Wuyi (III) were also relatively high, which might be related to nearby mines. The spatial distribution of Cu was similar to that of Cr. Furthermore, Cu had obvious high concentrations in the east of Yongkang (IX) and the southwest of Pan'an (V), which might be related to the fact that there are many hardware factories and plastic manufacturing factories there. The spatial distributions of Cd, Hg, Pb, and Zn concentrations were mainly island-shaped distributions, and showed obvious point source pollution, implying that they may be mainly affected by human activities [39]. High concentrations of Cd were mainly distributed in the western part of the study area, and were scattered in other counties. It has been demonstrated that long-term application of phosphate fertilizer, sewage irrigation, and mining and smelting of cadmium bearing ores can cause soil Cd pollution [9,49,50]. The high concentration of Hg mainly was distributed in a small area of the middle and western part of the study area, which might be caused by sewage irrigation and atmospheric deposition [49,51]. The high concentration of Pb was mainly distributed near the urban areas and main roads of the study area. Many studies indicated that the common sources of Pb in agricultural soils are industrial fumes, vehicle exhausts, and pesticides [39,52]. The high concentration of Zn showed a scattered distribution in the southern and central parts of the study area, which may have derived from both lithogenic sources and human activities [6,49,53].
part of the study area, which might be caused by sewage irrigation and atmospheric deposition [49,51]. The high concentration of Pb was mainly distributed near the urban areas and main roads of the study area. Many studies indicated that the common sources of Pb in agricultural soils are industrial fumes, vehicle exhausts, and pesticides [39,52]. The high concentration of Zn showed a scattered distribution in the southern and central parts of the study area, which may have derived from both lithogenic sources and human activities [6,49,53].

Environmental Risk of Soil HM Concentration
The distribution of soil HM hot spots and environmental risk probabilities are shown in Figure 3. Almost all the hot spots of HMs were present in areas with environmental risk, indicating that both the G * i index and indicator kriging could reflect the soil environmental risks of soil HMs well. However, As and Cu elements did not have hot spots in the highrisk area (VI), which might be because the G * i index is a commonly used local spatial autocorrelation method. An earlier study noted that the spatial autocorrelation analysis could easily compare and test the strength and significance of spatial dependence, and the local Moran's I is helpful to find the spatial outliers or hot spots of soil HMs, which are both considered effective exploratory spatial analysis methods for regional variables [54]. Hu et al. [19] reported that most interpolation methods have smoothing effects or prediction bias when predicting soil HM pollution, and the G * i index could not only define hot spots, but also provide credible probabilities for the results, which is of great significance for practical applications. However, the G * i index could only evaluate the known sample data, while the indicator kriging could perform risk assessment on unsampled points, which makes up for the deficiency of the G * i index [12,13,54]. Lin et al. [55] also demonstrated that the combination of hot spot analysis and geostatistical analysis was effective for exploring soil HM pollution.  (2) Basic safe utilization zone. The basic safe utilization zone accounted for 8.58% of the total farmland area and was mainly distributed in counties except Pujiang (IV). The soil HM concentrations in this zone slightly exceeded the screening values, but without constituting a threat to crop growth and human health. Therefore, priority protection strategies should be adopted, such as controlling the input of new pollutants and optimizing agricultural farming measures to reduce the risks of agricultural production [17].
(3) Low-risk monitoring zone. The area of the low-risk monitoring zone only accounted for 1.41% of the total farmland area, and was mainly distributed in the western and southern parts of Lanxi (VI), southern parts of Yiwu (VII), and southern parts of Wuyi (III), etc. The soil HM concentrations of farmland soil in this zone slightly exceeded the screening values, and pose a threat to crop growth and human health. Therefore, it is necessary to monitor the soil environment dynamically to prevent the input of pollution. Additionally, it is suggested that soil remediation measures such as adjusting the tillage system should be taken to improve the soil environment [4].
(4) Medium-risk early warning zone. The area of this zone was about 10.84 km 2 , and was mainly distributed in the western part of Lanxi (VI), Yiwu (VII) and the southern part of Wuyi (III). The concentration of soil HMs in these areas obviously exceeded the risk screening values, posing a greater threat to crop growth and human health. Therefore, strict early warning and control strategies should be implemented, and the quality of crops, soil and irrigation water should be strictly monitored. Additionally, economic

Safe Utilization Zoning of Soil HMs in Farmland
The IICQs was calculated according to Equations (6)- (9), and the grades of the farmland soil safe utilization zones were classified with reference to Table 1. The safe utilization zoning of soil HMs in farmland of Jinhua City were interpolated by ordinary kriging (Figure 4), and the areas of each safe utilization zone are shown in Table 4. HMs in the human body through the food chain [4]. (5) High-risk restricted zone. The high-risk restricted zone had the smallest area, and only accounted for 0.05% of the total farmland area in the study area, and was concentrated in the southern part of Wuyi (III). The high-risk restricted zone had extremely high levels of soil HM concentrations, which seriously exceeded the risk screening values and were not suitable for planting crops. Therefore, the area should have restricted agricultural use and the soil should be actively repaired [57].

Conclusions
In this paper, we explored the spatial distribution, environmental risk, and safe utilization of soil HMs in farmland of Jinhua City, subtropical China. The results showed that the mean concentrations of As (6.37 mg kg −1 ), Cd (0.18 mg kg −1 ), Cu (18.60 mg kg −1 ), Cr (36.73 mg kg −1 ), Hg (0.11 mg kg -1 ), Ni (11.95 mg kg -1 ), Pb (32.48 mg kg -1 ), and Zn (76.05 mg kg −1 ) did not exceed the risk screening values for soil contamination of agricultural land of China. However, the mean concentrations of Cu, Hg, and Zn were higher than the corresponding background values of the study area, indicating anthropic inputs of these elements were high. The C/(C0+C) ratios of As (0.576), Cr (0.502), Cu (0.502), and Ni (0.501) indicated a moderate spatial correlation, whereas the C/(C0+C) ratios of Cd (0.874), Hg (0.919), Pb (0.926), and Zn (0.904) indicated a strong spatial randomness and weak autocorrelation. As, Cr, Cu, and Ni showed strip-shaped and island-shaped distributions, while Cd, Hg, Pb, and Zn presented island-shaped distributions. The environmental risks of As, Cd, Cu, Pb, and Zn were relatively high, which need to be monitored. (1) Safe utilization zone. The safe utilization zone covered an area of 2144.47 km 2 , accounting for 89.35% of the total farmland area of Jinhua, and was evenly distributed in the northern part of Wucheng (I) and Wuyi (III) and other countries. The soil HM concentrations in the safe utilization area might exceed the soil background value, but did not exceed the risk screening values for soil contamination of agricultural land of China, so the soil was clean and safe to use, thus a strategy of key protection and prioritized use should be adopted. It is necessary to give priority to the development of agriculture and strictly prevent the input of pollutants [56].
(2) Basic safe utilization zone. The basic safe utilization zone accounted for 8.58% of the total farmland area and was mainly distributed in counties except Pujiang (IV). The soil HM concentrations in this zone slightly exceeded the screening values, but without constituting a threat to crop growth and human health. Therefore, priority protection strategies should be adopted, such as controlling the input of new pollutants and optimizing agricultural farming measures to reduce the risks of agricultural production [17].
(3) Low-risk monitoring zone. The area of the low-risk monitoring zone only accounted for 1.41% of the total farmland area, and was mainly distributed in the western and southern parts of Lanxi (VI), southern parts of Yiwu (VII), and southern parts of Wuyi (III), etc. The soil HM concentrations of farmland soil in this zone slightly exceeded the screening values, and pose a threat to crop growth and human health. Therefore, it is necessary to monitor the soil environment dynamically to prevent the input of pollution. Additionally, it is suggested that soil remediation measures such as adjusting the tillage system should be taken to improve the soil environment [4].
(4) Medium-risk early warning zone. The area of this zone was about 10.84 km 2 , and was mainly distributed in the western part of Lanxi (VI), Yiwu (VII) and the southern part of Wuyi (III). The concentration of soil HMs in these areas obviously exceeded the risk screening values, posing a greater threat to crop growth and human health. Therefore, strict early warning and control strategies should be implemented, and the quality of crops, soil and irrigation water should be strictly monitored. Additionally, economic crops could be planted instead of grain crops, which could reduce the accumulation of HMs in the human body through the food chain [4].
(5) High-risk restricted zone. The high-risk restricted zone had the smallest area, and only accounted for 0.05% of the total farmland area in the study area, and was concentrated in the southern part of Wuyi (III). The high-risk restricted zone had extremely high levels of soil HM concentrations, which seriously exceeded the risk screening values and were not suitable for planting crops. Therefore, the area should have restricted agricultural use and the soil should be actively repaired [57].

Conclusions
In this paper, we explored the spatial distribution, environmental risk, and safe utilization of soil HMs in farmland of Jinhua City, subtropical China. The results showed that the mean concentrations of As (6.37 mg kg −1 ), Cd (0.18 mg kg −1 ), Cu (18.60 mg kg −1 ), Cr (36.73 mg kg −1 ), Hg (0.11 mg kg −1 ), Ni (11.95 mg kg −1 ), Pb (32.48 mg kg −1 ), and Zn (76.05 mg kg −1 ) did not exceed the risk screening values for soil contamination of agricultural land of China. However, the mean concentrations of Cu, Hg, and Zn were higher than the corresponding background values of the study area, indicating anthropic inputs of these elements were high. The C/(C 0 + C) ratios of As (0.576), Cr (0.502), Cu (0.502), and Ni (0.501) indicated a moderate spatial correlation, whereas the C/(C 0 + C) ratios of Cd (0.874), Hg (0.919), Pb (0.926), and Zn (0.904) indicated a strong spatial randomness and weak autocorrelation. As, Cr, Cu, and Ni showed strip-shaped and island-shaped distributions, while Cd, Hg, Pb, and Zn presented island-shaped distributions. The environmental risks of As, Cd, Cu, Pb, and Zn were relatively high, which need to be monitored. The total area of safe and basic safe utilization areas accounts for 98% of farmland in the study area. For the safe and basic safe utilization area, a protection strategy should be implemented. The low-risk and medium-risk utilization areas should adopt dynamic monitoring and early warning prevention strategies, respectively, while the high-risk utilization areas should adopt strategies to restrict agricultural use. This study will provide a reference for land managers to formulate effective control and remediation strategies for protecting farmland soil. From the methodological perspective, this paper will provide a methodological reference for exploring soil HM pollution in areas with similar characteristics to the study area.

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.