An Uncertainty Assessment of Human Health Risk for Toxic Trace Elements Using a Sequential Indicator Simulation in Farmland Soils

Toxic trace elements in farmland soils are potential threats to human health. In this study, we collected soil samples from the farmlands of southern Guangzhou. We used a sequential indicator simulation (SIS) to deal with the problem of skewed distribution in the sample data. We assessed the human health risks, as well as the uncertainties, of five toxic trace elements: arsenic (As), cadmium (Cd), chromium (Cr), lead (Pb), and mercury (Hg). The results were as follows: (1) The risk indices of two trace elements (Cd and Hg) were less than the standard threshold, which means that there was no human health risk due to Cd and Hg in the study area. However, the maximum risk indices of As, Cr, and Pb exceeded the standard threshold. In particular, the maximum risk index of Pb was twice the standard threshold; (2) The risk probabilities of As and Cr were less than 25% in most areas, and only a few parcels of farmland have a 100% risk probability. The risk map of Pb was used to identify contiguous areas of high-risk probability (i.e., 75%–100%) in the center of the study area. (3) E-type estimation by the SIS method overestimates the risk when the number of samples with concentrations above the threshold have a large proportion of total samples. Our conclusions are as follows: (1) The simulation results show that areas with high-risk indices were concentrated in the Panyu District, which is close to the Pearl River and the core urban area of Guangzhou; (2) Except for Pb, these trace elements are not likely to pose health risks in southern Guangzhou; (3) This study considers the risk probability found with the SIS method to be more reliable for visualizing regional risk.

farmland can easily enter the human body through breathing, skin contact, and food intake [3], which can result in lesions or cancerous changes in the human body [4]. For example, cadmium (Cd) intake can cause diseases such as lung cancer, kidney dysfunction, and high blood pressure [5]. Lead (Pb) intake can cause symptoms such as lead poisoning, anemia, kidney disease, and gastrointestinal cramps [6]. Hexavalent chromium (Cr) is a recognized carcinogen [7]. The body's central nervous system is affected negatively by arsenic (As) and mercury (Hg) intake [8,9]. The long-term accumulation of toxic trace elements in farmland comprises a major threat to human health. Regional health risk status of toxic trace elements are important for the protection of human health, food security, and farmland environments [10]. In order to obtain regional health risks accurately and effectively, it is important to quantify trace element toxicity [11]. Human health risk assessment is a key method for risk quantification [12,13], which can help in the management of health risks, pollution warnings, and soil restoration [14].
Human health risk assessment calculates the risk of toxic trace elements by using the elemental concentrations of the samples [11,15]. There are two types of methods for human health risk assessment: methods based on deterministic models and methods based on stochastic models.
The first type of method is based on deterministic models. These methods construct a risk assessment model using the mechanisms of the exposure pathway, dose-response relationship, and crop reflection relationship. This type of model obtains the required parameters by referring to relevant institutional guidelines or scholarly research. Then, the health risk of the study area is calculated by the average hazard index with the soil trace element concentrations at the sampling points [16,17]. For example, Gu et al. [18] referred to the manual exposure factors and parameters in the technical guidelines for contaminated site risk assessment [19], in order to obtain human health risk assessments in 28 urban parks in Guangzhou. They calculated the average hazard index of the carcinogenic and non-carcinogenic risk of each park and concluded that no human health risk existed in the urban parks in Guangzhou. Sawut et al. [7] used a series of existing studies to determine the intake rate, intake time, and other parameters of trace elements. They calculated average hazard indices for both children and adults, and obtained the human health risks of vegetables in Northwest China. However, due to the individual differences of soil environmental variables and a lack of prior knowledge, trace element concentrations may have a high degree of spatial variability. The parameters of this model, therefore, have considerable uncertainty. Most of these parameters are determined directly by referring to existing guidelines or research. Moreover, the researchers calculated the average value of the hazard index to obtain the degree of regional risk, which may lead to smoothing effects. This method might overestimate or underestimate the risk. Some researchers have estimated the regional risk without considering the uncertainties involved [20].
The second type of method is based on stochastic models. These methods use stochastic simulation to calculate the uncertainty of the model parameters for human health risk assessment [21,22]. The stochastic simulation takes values randomly based on the cumulative distribution of the variable data, and produces multiple simulation results to obtain the risk and its uncertainty [23]. For example, Chen et al. [20] used a Monte Carlo simulation to solve the parameter uncertainty problem caused by a lack of prior knowledge when evaluating human health risk in the surface soil of Beijing. They obtained the human health risk situation and its uncertainty in the Beijing sampling area. Albuquerque et al. [24] used a sequential Gaussian simulation to obtain the uncertainty of the high and low-risk areas while studying the risk levels of 14 trace elements in the Aviles Estuary in Spain, and visualized the risk using a mean value map. However, data conforming to a normal distribution is a requirement for these methods. For instance, a Monte Carlo simulation may lead to more errors when the number of samples is not sufficient or the probability distribution is biased. In addition, a sequential Gaussian simulation loses information when the data are processed and transformed.
Human health risk assessment based on the uncertainty model can objectively reflect the risk of toxic trace elements. However, the current uncertainty model research is based on the condition that the frequency distribution of the sample data follow a normal distribution. Meeting the normal Sustainability 2020, 12, 3852 3 of 17 distribution requirements may be difficult when using sample data of trace elements in farmland, especially when the trace element concentrations have abnormal values. Simple abnormal value processing methods are prone to smoothing effects, which results in the loss of information in key areas [14,25]. Moreover, due to the cost limitations of sampling, the number of samples is often limited. Abnormal values influence the frequency distribution, but often contain information about contamination by toxic trace elements. Therefore, we cannot simply eliminate the abnormal values. Therefore, we need to study the use of uncertainty models in obtaining human health risk when the frequency distribution of the sample data has a non-normal distribution.
To obtain the human health risk of toxic trace elements in farmland, this study uses a non-parametric statistical sequential indicator simulation method to address the problem of the sample data frequency distribution. This study focuses on (1) simulating the regional farmland health risk of each trace element; (2) assessing their uncertainty; and (3) mapping the spatial distribution of the human health risk probability.

Study Area
The accumulation of toxic trace elements in farmland, especially around megacities, may result in a high probability of human health risks [18,20]. Guangzhou, a rapidly growing city, is one of the largest cities in South China. Guangzhou is located in the center of Guangdong Province, and is adjacent to Hong Kong and Macau. Guangzhou has a maritime subtropical monsoon climate, with an annual average temperature of 22.8 • C and average annual precipitation of 1800 mm. The Pearl River system passes through Guangzhou. It has a population of over 20 million. Human activity has a huge impact on this city. Due to the high degree of urbanization, few farmlands remain in the main urban areas of Guangzhou. Therefore, we selected the southern part of Guangzhou (i.e., Panyu and Nansha District) as the study area in this study, which has a lower urbanization level and contiguous farmland.

Soil Sampling and Laboratory Analysis
We designed the sampling plan according to land use/cover data of the study area. We sampled in the farmland as evenly as possible, with 195 samples (0-20 cm) throughout the farmland in 2015. Figure 1 shows the spatial distribution of the sample points.  In this study, the selection of trace elements was based on the existing toxic trace element contamination research of Guangzhou [26][27][28]; we also refer to the "12th Five-Year Plan for Heavy Metals Contamination Comprehensive Prevention and Control" of China, and the "13th Five-Year Plan for Heavy Metals Contamination Comprehensive Prevention and Control" of Guangdong Province. We selected five trace elements for study: lead (Pb), cadmium (Cd), chromium (Cr), mercury (Hg), and arsenic (As). We determined the total amounts of these five trace elements. We used the HF-HNO 3 -HClO 4 digestion method to process lead (Pb), cadmium (Cd), and chromium (Cr). We determined Cd and Pb by graphite furnace atomic absorption spectrometry and we determined Cr by flame atomic absorption spectrophotometry. The atomic absorption spectrometer (graphite furnace) used was the Z-2000 (HITACHI). The flame atomic absorption spectrometer used was the Z-5300 (HITACHI). Arsenic (As) and mercury (Hg) were processed with HCl-HNO 3 and were determined by reduction-gasification atomic fluorometry. The atomic fluorescence photometer used was the AFS-930 double-pipe atomic fluorescence spectrometer (Beijing Titan Instruments).

Sequential Indicator Simulation of Soil Data
Sequential indicator simulation is a stochastic simulation method based on non-parametric statistics. Its indicator model has good flexibility and does not require the data used to satisfy any certain distribution type [29][30][31]. Moreover, the simulation has a small smoothing effect. This method does not ignore abnormal values in areas, which makes the simulation result closer to the real situation [14]. Information about the extreme values in the abnormal region can be utilized [32].According to the existing research, SIS can be divided into the following four steps [29]: (1) The original sampled data are indicator transformed. The user obtains a priori conditional cumulative distribution function (CCDF) according to the Kriging method. Z(x i ), i = 1, 2, . . . , n is a set of sampled data under the K-group threshold (z k ). The more quantiles (thresholds) that are selected, the more reliable the CCDF is [33]. The indicator variables are encoded as 0 or 1 by the indicator function I(x i ) (Equation (1)) [29]: (2) Using the known raw data as conditional data, the user obtains the conditional prior probability F(x) = Z(x i ) ≤ z k C(x i ) . This probability uses the indicator Kriging estimate to extrapolate its conditional cumulative distribution function F z(x) .
(3) The user divides the study area into grids of uniform resolution and defines a random path through all grid nodes. This method extracts a value from F z(x) randomly as the simulation value at the first grid position. Then, the user obtains the estimated valueÎ of the CCDF at the target grid (Equations (2)) [29]: (4) The above new value is added to the conditional data set. The grid's prior condition cumulative distribution function F z(x) for the next simulation position is modeled by this new data set. This process is repeated until all grid nodes are simulated.
This method obtains n realizations by repeating the above-mentioned steps. In n realizations, the number of grid values exceeding the threshold is n 1 and the uncertainty (above threshold probability) of this grid is n 1 /n. Each realization of the sequential indicator simulation has equal probability [25,32]. Some studies have also used conditional means (E-type estimate) of n realizations to obtain the spatial distribution of trace element concentrations [29][30][31]. In the calculation process, we used the GS + 9.0 software to obtain the experimental semivariogram and variogram models of each threshold. We realized the stochastic simulation with the SGeMS v2.5b software.

Human Health Risk Assessment
Human health risks contain carcinogenic and non-carcinogenic risks. Toxic trace elements enter into the human body in different ways, such as oral intake, skin contact, and inhalation [34]. Existing research has demonstrated that As and Cd have carcinogenic effects [35]. Hexavalent chromium (Cr 6+ ) is carcinogenic to humans, but trivalent chromium (Cr 3+ ) has not been defined as being carcinogenic yet [35]. However, Cr 6+ is more unstable than Cr 3+ in nature [36]. Chromium in the soil is more inclined towards the stable state of Cr 3+ [37]. Therefore, the reference parameter is the parameter of Cr 3+ in this study, and we only consider the non-carcinogenic risk of Cr. The non-carcinogenic assessment of trace elements include Hg, Pb, and Cr. We used an exposure model based on the US EPA model [38]. The equations (Table A1) and parameters (Tables A2 and A3) refer to the Technical Guidelines For Risk Assessment of Contaminated Sites and Regional Screening Levels [19,38]. Due to the different dietary and living habits of people in different countries, the exposure parameters of health risk are different. The Technical Guidelines For Risk Assessment of Contaminated Sites was used to correct some parameters, according to the actual situation in China [19].

Uncertainty Model Accuracy Evaluation
The accuracy of the uncertainty model was evaluated by the symmetric p-probability intervals (p-PI) [29]. The conditional cumulative probability distribution function F(u; ζ) at any location u allows for the computation of p-PI bounded by (1-p)/2 and (1+p)/2. For example, the 0.5-PI representation is In this case, the correct uncertainty model would show that there is a 0.5 probability that the actual observations at u falls into this interval [29]. We calculated the probability that an observation falls into p-PIs as follows (Equations (11) and (12)) [29]: Sustainability 2020, 12, 3852 6 of 17 A smaller root mean square error (RMSE) means higher accuracy [39]. The RMSE was calculated as follows (Equation (13)) [39]: Table 1 shows the results of the descriptive statistics of five trace elements in southern Guangzhou farmland. The main soil types of farmland were paddy soils and red soils in study area. The main soil texture was loam. The clay content was in the range of 22%-43%. The maximum, average, and minimum pH values were 8.2, 5.7, and 3.9, respectively. By comparing the background values of risk screening values for soil trace elements of the Pearl River Delta area [40], the maximum values of the five trace elements were seen to be higher than the soil background values. This indicated that trace elements have external inputs in the study area. Furthermore, the average values for Cd and Hg were larger than the soil background values, indicating that the Cd and Hg accumulation was high. Previous studies have found that urban land and mining areas had contamination of toxic trace elements in Guangzhou [18]. Soil trace elements of farmland in Guangzhou are likely to pose risks to crops, humans, and the environment. Due to the high urbanization development in Guangzhou, as seen in Figure 1, farmland in the Panyu and Nansha Districts is very close to urban areas. Therefore, studies of the human health risks of toxic trace elements in this farmland has been deemed necessary.

Preliminary Data Description
The coefficient of variation (CV) is the ratio of the standard deviation to the mean [41]. The CV reflects the degree of dispersion of variables. The coefficient of variation of the five trace elements was in the range of 32%-81%. The order of variation was As > Hg > Cd > Cr = Pb. The coefficients of variation of As and Hg were relatively large, indicating that these two trace elements may be greatly affected by human activity. The asymmetry coefficients (skewness) of the five trace elements were all greater than 0, as can be seen in Table 1; therefore, their data frequency distributions had a certain positive bias. The order of deviation was As > Hg > Pb > Cd > Cr. The Kolmogorov-Smirnov (K-S) test showed that As and Hg had a significance of less than 0.05, which means that their data distribution was not normal. We could not use those data to provide a good result with some uncertainty models, such as Monte Carlo or a sequential Gaussian simulation [42].

Human health Risk Threshold of Each Trace Element
We used the concentrations of trace elements to calculate human health risks at each sample point. Then, we calculated the average values of sample points to obtain the human health risk indices (CR or HI) for the five trace elements. According to the risk standard value, a cancer risk index (CR) below 10 −4 is acceptable, and the non-carcinogenic hazard index (HI) indicates no hazard when it has a value Sustainability 2020, 12, 3852 7 of 17 less than 1.0. As shown in Table 2, the CR of As and Cd were less than 10 −4 , and the HI of Hg, Pb, and Cr were less than 1.0. This result meant that the carcinogenic risk of the study area was acceptable and non-carcinogenic hazards did not exist in the study area. However, this method may ignore the risk information of abnormal value of some samples. According to these judgments and the US EPA's health risk model, we calculated the CR and HI risk thresholds for each trace element. The threshold of concentration of a carcinogenic trace element was calculated as C sur = 10 −4 / (Oral risk factor + Dermal risk factor + Inhalational risk factor). The threshold of concentration of a non-carcinogenic trace element was calculated as C sur = 1/ (Oral risk factor + Dermal risk factor + Inhalational risk factor). Table 2 shows those thresholds. As can be seen in Table 1, the concentrations of As, Pb, and Cr exceeded their risk thresholds. Therefore, certain risks existed in some areas of the study area.

Human Health Risk Assessment Based on SIS
We used the thresholds (Table 2) to perform a sequential indicator simulation. This method needs the a priori information of the corresponding indicator Kriging. As the distribution of sample points was non-uniform, we used the average nearest neighbor method to determine the step size. According to the rule of thumb, the number of steps is equal to half the maximum distance between pairs of sample points divided by the step size. We took the quartile of each trace element concentrations as the threshold value, and obtained experimental semivariogram and variogram models of each threshold condition. Table A4 lists the parameters of the experimental semivariograms and the fitted models for each trace element. According to the size of the study area and the operation time of grid division, we divided the study area into of 620 × 510 grid, where the size of each grid was 100 m × 100 m. We carried out 1000 realizations in each grid by SIS. Then, we obtained the risk index of CR or HI for each trace element by calculating the E-type estimates. The E-type estimates showed the spatial distribution of the risk index ( Figure 2).
As shown by the risk index distribution map (Figure 2), the clusters with high values of each trace element were mostly concentrated in the north and central regions of the study area. These areas were mainly close to the Pearl River and the urban areas of the Panyu District (Figure 1). This showed that urban areas have an impact on trace elements of farmland around cities. The cancer risk index (CR) of As exceeded the threshold (10 −4 ) in some areas. The cancer risk index (CR) of Cd did not exceed the threshold (10 −4 ), which meant there was no risk of Cd in the study area. The maximum non-carcinogenic hazard index (HI) of Cr and Pb exceeded the threshold (1.0). In particular, the maximum non-carcinogenic HI of Pb was twice the standard threshold. The non-carcinogenic hazard index (HI) of Hg did not exceed the threshold (1.0), which meant there was no risk of Hg in the study area. According to the spatial distribution of risk index, the risk areas comprised only a small part of the study area. As shown by the risk index distribution map (Figure 2), the clusters with high values of each trace element were mostly concentrated in the north and central regions of the study area. These areas were mainly close to the Pearl River and the urban areas of the Panyu District (Figure 1). This showed that urban areas have an impact on trace elements of farmland around cities. The cancer risk index (CR) of As exceeded the threshold (10 −4 ) in some areas. The cancer risk index (CR) of Cd did not exceed the threshold (10 −4 ), which meant there was no risk of Cd in the study area. The maximum non-carcinogenic hazard index (HI) of Cr and Pb exceeded the threshold (1.0). In particular, the maximum non-carcinogenic HI of Pb was twice the standard threshold. The non-carcinogenic hazard index (HI) of Hg did not exceed the threshold (1.0), which meant there was no risk of Hg in the study area. According to the spatial distribution of risk index, the risk areas comprised only a small part of the study area.

Uncertainty Assessment
From the above-mentioned CR or HI of trace elements, Cd and Hg were unlikely to pose a risk to human health in the study area. However, the maximum risk values of As, Pb, and Cr exceeded the threshold. We calculated the probability of each grid exceeding the risk for As, Pb, and Cr. The process was as follows: the risk threshold of each trace element was set as the criterion. Then, we counted the number (n1) of each grid when the trace element concentration exceeded the risk threshold in n realizations; the risk probability of the grid was calculated with p=n1/n. Figure 3 shows the spatial distribution of the risk probability.

Uncertainty Assessment
From the above-mentioned CR or HI of trace elements, Cd and Hg were unlikely to pose a risk to human health in the study area. However, the maximum risk values of As, Pb, and Cr exceeded the threshold. We calculated the probability of each grid exceeding the risk for As, Pb, and Cr. The process was as follows: the risk threshold of each trace element was set as the criterion. Then, we counted the number (n 1 ) of each grid when the trace element concentration exceeded the risk threshold in n realizations; the risk probability of the grid was calculated with p = n 1 /n. Figure 3 shows the spatial distribution of the risk probability.   Figure 3 shows that most of the risk probabilities of As and Cr ranged from 0% to 25%, which indicates a low-risk level. However, areas of high-risk probability (75% to 100%) also existed. The 100% risk probabilities for As and Cr only existed in a few grids. According to the size of the grid, a few parcels of farmland had certain health risks. Therefore, the overall risks of As and Cr were acceptable in the study area. We should give much more attention to the at-risk farmlands. Figure 3 Figure 3 shows that most of the risk probabilities of As and Cr ranged from 0% to 25%, which indicates a low-risk level. However, areas of high-risk probability (75% to 100%) also existed. The 100% risk probabilities for As and Cr only existed in a few grids. According to the size of the grid, a few parcels of farmland had certain health risks. Therefore, the overall risks of As and Cr were acceptable in the study area. We should give much more attention to the at-risk farmlands. Figure 3 also shows the risk probability of Pb. Most of the risk probabilities of Pb ranged from 25% to 50%, which indicates a moderate-risk level. The high-risk probability (75% to 100%) areas of Pb were concentrated in the central part of the study area. Furthermore, areas of 100% risk probability of Pb existed in some grids. The management should give priority to and carry out monitoring in these farmlands.
According to the grid statistics (Figure 4), most risk possibilities of As and Cr were less than 25%. Therefore, most areas were low-risk situations, which indicates that these trace elements were less likely to pose health risks in southern Guangzhou. The proportion of moderate-risk areas (25% < probability < 75%) for Pb was 80.06% in the study area. The proportion of high-risk areas (probability > 75%) for Pb was 1.2%, with an area of 2.8 km 2 . Therefore, the study area had a moderate-risk of Pb. The risk probability can help to determine the key risk control areas of farmland. The management should set a certain monitoring threshold, and consider more sampling and analysis of toxic trace elements in these farmlands. It is necessary to take more protection measures in these farmlands.  Figure 3 shows that most of the risk probabilities of As and Cr ranged from 0% to 25%, which indicates a low-risk level. However, areas of high-risk probability (75% to 100%) also existed. The 100% risk probabilities for As and Cr only existed in a few grids. According to the size of the grid, a few parcels of farmland had certain health risks. Therefore, the overall risks of As and Cr were acceptable in the study area. We should give much more attention to the at-risk farmlands. Figure 3 also shows the risk probability of Pb. Most of the risk probabilities of Pb ranged from 25% to 50%, which indicates a moderate-risk level. The high-risk probability (75% to 100%) areas of Pb were concentrated in the central part of the study area. Furthermore, areas of 100% risk probability of Pb existed in some grids. The management should give priority to and carry out monitoring in these farmlands. According to the grid statistics (Figure 4), most risk possibilities of As and Cr were less than 25%. Therefore, most areas were low-risk situations, which indicates that these trace elements were less likely to pose health risks in southern Guangzhou. The proportion of moderate-risk areas (25% < probability < 75%) for Pb was 80.06% in the study area. The proportion of high-risk areas (probability

Accuracy Evaluation of the Uncertainty Model
We evaluated the accuracy of the risk uncertainty. Nine p-PI intervals ("10%", "20%", "30%", "40%", "50%", "60%", "70%", "80%", and "90%") were calculated to obtain the probability that the sample point fell into the p interval. The closer the scatter is to the standard line (i.e., y = x), the more accurate the result. The results in Figure 5 show the accuracy of the evaluation results for the risk probability (CR or HI). The scattered lines are all close to the standard line. The RMSE level was also low: the RMSE of As, Pb, and Cr were 0.0218, 0.0208, and 0.0211, respectively. The RMSE indicated that the overall accuracy of the uncertainty evaluation was effective [39].
"40%", "50%", "60%", "70%", "80%", and "90%") were calculated to obtain the probability that the sample point fell into the p interval. The closer the scatter is to the standard line (i.e., y = x), the more accurate the result. The results in Figure 5 show the accuracy of the evaluation results for the risk probability (CR or HI). The scattered lines are all close to the standard line. The RMSE level was also low: the RMSE of As, Pb, and Cr were 0.0218, 0.0208, and 0.0211, respectively. The RMSE indicated that the overall accuracy of the uncertainty evaluation was effective [39].

Comparison of Human Health Risk Results
The coefficient variation (CV) of As was large, which means that the accumulation of As was mainly caused by human factors. The areas with high concentrations of As were close to the Pearl River basin in this study. The Panyu and Nansha District are delta areas where the sediments of the Pearl River accumulate, and the irrigation of the farmland in these areas depends on the Pearl River and its tributaries. The discharge of sewage and wastewater results in the increase of As concentrations in the Pearl River [43,44], which may cause the accumulation of As. At the same time, the productive process of Guangzhou's non-ferrous metal production is the main source of As in the atmosphere. Dust deposition has a negative impact on outdoors environments, and the CR of As has been found to exceed the threshold in urban areas of Guangzhou [45], which is consistent with the results of this study. However, we also found that most of the risk of As did not exceed the unacceptable level in farmland of Guangzhou, and the risk management can be strengthened according to the risk probability.
From the results of the risk assessment for Cr, the risk probabilities of most farmlands were at a low-risk level (0%-25%). The health risk of Cr was found to be acceptable in Guangzhou City in other studies [18,45,46], in accordance with the results of this study. Very few farmlands exhibited a highrisk situation, and the spatial distribution did not show any tendency for aggregation. Therefore, we

Comparison of Human Health Risk Results
The coefficient variation (CV) of As was large, which means that the accumulation of As was mainly caused by human factors. The areas with high concentrations of As were close to the Pearl River basin in this study. The Panyu and Nansha District are delta areas where the sediments of the Pearl River accumulate, and the irrigation of the farmland in these areas depends on the Pearl River and its tributaries. The discharge of sewage and wastewater results in the increase of As concentrations in the Pearl River [43,44], which may cause the accumulation of As. At the same time, the productive process of Guangzhou's non-ferrous metal production is the main source of As in the atmosphere. Dust deposition has a negative impact on outdoors environments, and the CR of As has been found to exceed the threshold in urban areas of Guangzhou [45], which is consistent with the results of this study. However, we also found that most of the risk of As did not exceed the unacceptable level in farmland of Guangzhou, and the risk management can be strengthened according to the risk probability.
From the results of the risk assessment for Cr, the risk probabilities of most farmlands were at a low-risk level (0%-25%). The health risk of Cr was found to be acceptable in Guangzhou City in other studies [18,45,46], in accordance with the results of this study. Very few farmlands exhibited a high-risk situation, and the spatial distribution did not show any tendency for aggregation. Therefore, we should analyze the high-risk fields, in order to determine if they are associated with point-source pollution.
Traffic activities have been shown to have a serious impact on the Pb pollution in the urban surface soil of Guangzhou [28,47]. The concentration of Pb in the soil has been positively correlated with the traffic load of the road [48]. The spatial distribution of Pb's human health risk in farmland was scattered in this study. Most grids of risk probabilities (50%-100%) were near roads. In addition, most of these roads are expressways. Furthermore, Pb-containing substances enter the ground through industrial emission, coal combustion, and waste incineration [49,50]. Those problems also affect the farmland soil of this city. In this study, Pb is a toxic trace element which can pose health risks in some farmlands. Management should monitor the tendency of Pb accumulation to control the risk.
We found that there has no human health risk of Cd and Hg in the farmland soil of southern Guangzhou under the current conditions. In other studies on toxic trace elements in Guangzhou, the areas with high concentrations of Cd and Hg appeared in the core urban area of the city [51]. Those core urban areas are close to the north of the study area. At the same time, the human health risks of Cd and Hg in other studies were found to be acceptable for Guangzhou's urban environment and atmospheric particles [18,52,53]. These conclusions were similar to our results.

Selection of Risk Mapping Based on the SIS Method
The SIS method can obtain the probability (above threshold) and conditional means (E-type estimate) for each variable grid [29]. Some studies have obtained the spatial distribution of regional trace element concentrations by using E-type estimates [3,14,31]. According to the risk threshold, the spatial distribution of regional trace element concentrations also can provide a reference for regional risk management [3,31].
We considered three kinds of numbers: the number of sample points whose sample concentration exceeded the risk threshold, the number of grids whose E-type estimated value exceeded the risk threshold, and the number of grids with 100% risk probability (Table 3). Combined with the spatial distribution of risk probability (Figure 3), locations of the sample concentration above the threshold were consistent with 100% risk probability locations. The numbers were equal or close for sample points and 100% risk probability grids. At the same time, the number and distribution ( Figure 6) of over-risk grids (E-type estimate) of As and Cr were consistent with the over-risk sample point and the 100% risk probability grid. However, the number of over-risk grids (E-type estimate) of Pb were much greater than the 100% risk probability grid. This shows that these grids have a low risk probability but have a high E-type estimate value: which means that the high values in the grid have a small quantity, but a large contribution of the E-type estimate.  The risk threshold may have an impact on the accuracy of risk judgement made by the E-type estimate. Before the calculation of the E-type estimate, the sample data is divided into sub-intervals by quantiles. We then multiply the subinterval probability by the average of the subinterval to obtain the E-type estimate. The average value of the sub-interval is calculated by the head and tail values of these sub-intervals [54][55][56]. The number of sampling points whose sample concentration (Pb) exceeded the risk threshold was 42% of the total. We divided the Pb data into four sub-intervals (by quartile). We can see that the risk threshold was lower than the average of the third sub-interval (i.e., 2 nd quartile to 3 rd quartile). This is the reason why the E-type estimate is easily greater than the risk The risk threshold may have an impact on the accuracy of risk judgement made by the E-type estimate. Before the calculation of the E-type estimate, the sample data is divided into sub-intervals by quantiles. We then multiply the subinterval probability by the average of the subinterval to obtain the E-type estimate. The average value of the sub-interval is calculated by the head and tail values of these sub-intervals [54][55][56]. The number of sampling points whose sample concentration (Pb) exceeded the risk threshold was 42% of the total. We divided the Pb data into four sub-intervals (by quartile). We can see that the risk threshold was lower than the average of the third sub-interval (i.e., 2nd quartile to 3rd quartile). This is the reason why the E-type estimate is easily greater than the risk threshold. Therefore, when the number of samples with concentration above threshold has a large proportion of total sample amount, the E-type estimate overestimates the risk. Therefore, we tend to use the risk probability by using the SIS for risk judgement.

Conclusions
This study assessed the human health risks of toxic trace elements in the farmland of southern Guangzhou, a rapidly developing area of China. We also analyzed the risk uncertainty and visualized the regional risk probability. We drew the following conclusions.
The high-value areas of the human health risk index were concentrated in the north-central area (Panyun District) of the study area, which is close to the main urban area of Guangzhou, and has frequent human activities. We found that high urbanization poses a hazard to human health in farmland. We calculated the risk index of each trace element by using a risk assessment model of human health. The maximum risk indices of Cd and Hg were less than the standard threshold, while the maximum risk indices of As, Cr, and Pb exceeded the standard threshold. In particular, the maximum risk index of Pb was twice the standard threshold.
We visualized the risk uncertainty to obtain the risk probability distribution map. Most of the risk probabilities of As and Cr were between 0% and 25%, which indicates a low-risk level. The 100% human health risk probability areas of As and Cr were only a few parcels of farmland. Most of the risk probabilities of Pb were between 25% and 50%, which indicates a moderate-risk level. Some areas (2.8 km 2 ) had a 100% probability health risk of Pb; we should give more attention and take more protection measures for farmland in these areas. The RMSE of the risk uncertainty model was low. The order of RMSE was As > Pb > Cr and the accuracy of the uncertainty model was acceptable.
We found that the mean estimate using traditional methods may underestimate the human health risks. Moreover, the E-type estimate by the SIS method overestimated the human health risk when the number of samples with concentration above the threshold had a large proportion of total sample amount. Therefore, we suggest using the risk probability by the SIS method to visualize regional human health risks. Using the risk probability makes interpreting the risk situation more intuitive. It is relatively reliable and realistic, and can provide guidance for the prevention of toxic trace elements influencing human health.   A   Table A1. Equations of each parameter.

Parameter Type Equation
Carcinogenic exposure value Non-carcinogenic exposure value Note: OISER means oral ingestion of soil exposure. DCSER means soil exposure through skin contact routes. PISER means soil exposure to inhaled particulate matter.  Note: SFo means oral intake carcinogenic slope factor. SFd means skin contact carcinogenic slope factor. SFi means inhalation carcinogenic slope factor. RfDo means oral intake reference dose. RfDd means skin contact reference dose. RfDi means inhalation reference dose.