<figure><img src="image_1.png" /></figure>

HAL
open science

# Improved Mapping of Potentially Toxic Elements in Soil via Integration of Multiple Data Sources and Various Geostatistical Methods

Fang Xia, Bifeng Hu, Youwei Zhu, Wenjun Ji, Songchao Chen, Dongyun Xu, Zhou Shi

<figure><img src="image_2.png" /></figure>

## To cite this version:

Fang Xia, Bifeng Hu, Youwei Zhu, Wenjun Ji, Songchao Chen, et al., Improved Mapping of Potentially Toxic Elements in Soil via Integration of Multiple Data Sources and Various Geostatistical Methods. Remote Sensing, 2020, 12 (22), 10.3390/rs12223775. [https://hal-03356743](https://hal.inrae.fr/hal-03356743)

HAL Id: hal-03356743

https://hal.inrae.fr/hal-03356743

Submitted on 28 Sep 2021

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_4.png" /></figure>

Article

# Improved Mapping of Potentially Toxic Elements in Soil via Integration of Multiple Data Sources and Various Geostatistical Methods

Fang Xia<sup>1,2,3</sup>, Bifeng Hu<sup>3,4,5,6</sup>①, Youwei Zhu<sup>7</sup>, Wenjun Ji<sup>8</sup>, Songchao Chen<sup>3,9</sup>①, Dongyun Xu<sup>3,4</sup> and Zhou Shi<sup>3,4,*</sup>

Received: 27 September 2020; Accepted: 4 November 2020; Published: 17 November 2020

<figure><img src="image_5.png" /></figure>

check for
updates

**Abstract:** Soil pollution by potentially toxic elements (PTEs) has become a core issue around the world. Knowledge of the spatial distribution of PTEs in soil is crucial for soil remediation. Portable X-ray fluorescence spectroscopy (p-XRF) provides a cost-saving alternative to the traditional laboratory analysis of soil PTEs. In this study, we collected 293 soil samples from Fuyang County in Southeast China. Subsequently, we used several geostatistical methods, such as inverse distance weighting (IDW), ordinary kriging (OK), and empirical Bayesian kriging (EBK), to estimate the spatial variability of soil PTEs measured by the laboratory and p-XRF methods. The final maps of soil PTEs were outputted by the model averaging method, which combines multiple maps previously created by IDW, OK, and EBK, using both lab and p-XRF data. The study results revealed that the mean PTE content measured by the laboratory methods was as follows: Zn (127.43 mg kg⁻<fnref n="1" />) > Cu (31.34 mg kg⁻<fnref n="1" />) > Ni (20.79 mg kg⁻<fnref n="1" />) > As (10.65 mg kg⁻<fnref n="1" />) > Cd (0.33 mg kg⁻<fnref n="1" />). p-XRF measurements showed a spatial prediction accuracy of soil PTEs similar to that of laboratory analysis measurements. The spatial prediction accuracy of different PTEs outputted by the model averaging method was as follows: Zn (R<fnref n="2" /> = 0.71) > Cd (R<fnref n="2" /> = 0.68) > Ni (R<fnref n="2" /> = 0.67) > Cu (R<fnref n="2" /> = 0.62) > As (R<fnref n="2" /> = 0.50). The prediction accuracy of the model averaging method for five PTEs studied herein was improved compared with that of the laboratory and p-XRF methods, which utilized individual geostatistical methods (e.g., IDW, OK, EBK). Our results proved that p-XRF was a reliable alternative to the traditional laboratory analysis methods for mapping soil PTEs. The model averaging approach improved the prediction accuracy of the soil PTE spatial distribution and reduced the time and cost of monitoring and mapping PTE soil contamination.

<footnote n="1">¹ College of Economics and Management, Zhejiang A&F University, Hangzhou 311300, China;
xiaf@zafu.edu.cn</footnote>

<footnote n="2">² Zhejiang Province Key Cultivating Think Tank—Research Academy for Rural Revitalization of Zhejiang Province, Zhejiang A&F University, Hangzhou 311300, China</footnote>

<footnote n="3">³ Institute of Agricultural Remote Sensing and Information Technology Application, College of Environmental and Resource Sciences, Zhejiang University, Hangzhou 310058, China; bifeng.hu@inra.fr (B.H.); Songchao.Chen@inrae.fr (S.C.); xudongyun@zju.edu.cn (D.X.)</footnote>

<footnote n="4">⁴ Key Laboratory of Environment Remediation and Ecological Health, Ministry of Education, College of Environmental and Resource Sciences, Zhejiang University, Hangzhou 310058, China</footnote>

<footnote n="5">⁵ Sciences de la Terre et de l'Univers, Orléans University, 45067 Orléans, France</footnote>

<footnote n="6">⁶ Unité de Recherche en Science du Sol, INRA, 45075 Orléans, France</footnote>

<footnote n="7">⁷ Protection and Monitoring Station of Agricultural Environment, Bureau of Agriculture, Hangzhou 310020, China; zhuyw@zjagri.gov.cn</footnote>

<footnote n="8">⁸ College of land Science and Technology, China Agricultural University, Beijing 100085, China;
wenjun.ji@cau.edu.cn</footnote>

<footnote n="9">⁹ Unité InfoSol, INRAE, 45075 Orléans, France</footnote>

<footnote>* Correspondence: shizhou@zju.edu.cn</footnote>
**Keywords:** proximal soil sensing; potentially toxic elements; p-XRF; empirical Bayesian kriging; model averaging method

# 1. Introduction

Soil pollution by potentially toxic elements (PTEs) is one of the most important environmental problems nowadays around the world and it has attracted much public attention [1-7], especially in countries like China [8-10], India [11-13], Iran [14-16], Pakistan [17-19], Brazil [20-22], and Bangladesh [23-25], which are undergoing industrialization. Among these PTEs, As, Cd, Cu, Ni, and Zn were among the list of most concerned pollutants in soil around the world since they could reside in soils for a long time and could pose a great threat to human health acting via different pathways (e.g., derma contact, inhalation, and ingestion of PTEs-polluted food) [26-30]. For example, long-time exposure to Cd could increase the risk of lung cancer, bone fractures, and kidney dysfunction [31]. The negative chronic effects of As include kidney cancer, lung cancer, and liver cancer [32].

PTEs in the soil are caused by various sources, such as waste emission, smelting, mining activity, traffic exhaust, atmospheric deposition, wastewater irrigation, fertilizers, and pesticide application [33-37]. Regulation and remediation of soil PTE pollution is a challenging and urgent task.

Information on the soil PTE spatial patterns forms the basis for effective management and soil remediation. Accurate identification of the extent of polluted areas could contribute to the more efficient management of soils contaminated by PTEs and reduce health risks and financial costs. Interpolation methods such as inverse distance weighting (IDW) and ordinary kriging (OK) have been widely employed to explore the spatial distribution of PTEs [38-43]. Recently, a new empirical Bayesian kriging (EBK) method developed by Krivoruchko (2012) [44] has been introduced to estimate the spatial pattern of target variables in the fields of public health [45], meteorology [46], and geology [47]. EBK has the advantage of representing the stochastic spatial process locally as a stationary or non-stationary random field, where the parameters of the locally defined random field vary across space. To the best of our knowledge, no application of EBK to soil PTE mapping has been published.

All geostatistical methods mentioned above require the input of information from soil observations. Typically, many soil samples are needed to ensure estimation accuracy, especially for surveys conducted at large spatial scales. However, this effort is usually hindered by financial and time constraints because soil properties significantly vary at both spatial and temporal scales, requiring extensive surveys and ample soil sample sizes [48-51]. Therefore, a trade-off has to be made between minimizing the sampling size and maximizing the estimation accuracy.

Portable X-ray fluorescence spectroscopy (p-XRF) is a rapid, low-cost, non-destructive, and easily applicable analysis method that can measure soil PTEs and has been previously used to monitor PTE pollution in soil [52-55]. p-XRF allows for significantly reducing the cost and time for the measurement of PTEs, as it does not require chemical analysis procedures (e.g., acid digestion). These advantages of p-XRF could potentially enhance the spatial prediction of soil PTEs, with some attempts previously reported. For example, in a study conducted in a former smelting area, Kim et al. (2019) [56] proved that using p-XRF measurements as an auxiliary for COK improved the estimation of the PTE spatial distribution even with a smaller sample size. In another study, Xia et al. (2019) [57] reported that using p-XRF measurements as covariates in COK improved the PTE mapping accuracy in agricultural soils.

Scientists have produced multiple maps based on both laboratory and p-XRF data, as well as by different interpolation methods (i.e., IDW, OK), but until now, no single method has shown better performance than others for the whole range of soil PTEs [38]. However, no method has shown much better performance for different PTEs. Therefore, we suggest combining the results of different methods to obtain better performance results. Herein, we assumed that a PTE distribution map produced using different methods and data sources was complementary and, thus, more accurate. We merged different maps using the model averaging method [58]. Previous researchers have improved estimation
accuracy for estimating the spatial distribution of different soil properties, such as soil pH [59] and soil carbon [60], by applying model averaging. However, no study has been reported yet on employing the model averaging method to improve the mapping of soil PTEs using different data sources and interpolation methods.

To fill these gaps, we aimed to: (1) compare the performance of different interpolation methods (i.e., IDW, OK, and EBK) to estimate the PTE spatial distribution based on the laboratory and p-XRF measurements; (2) test the feasibility of using p-XRF as an alternative to traditional laboratory analyses to map soil PTEs; and (3) evaluate the performance and gain of model averaging based on a map produced using different data sources and interpolation methods.

## 2. Materials and Methods

### 2.1. Study Area and Soil Sampling

This survey was conducted in Fuyang County, Hangzhou, South China. The study area covers an area of 1821.03 km² (latitude 29°44′45″–30°11′58.5″N and longitude 119°25′–120°19′30″E) (Figure 1). Fuyang County is characterized by a subtropical monsoon climate, with a mean annual precipitation of approximately 1487 mm and a mean annual temperature of 17.8 °C. We collected 293 surface soil (0–20 cm) samples from cropland in the survey region using a stainless steel shovel, where each soil sample was mixed with five other soil samples within a radius of 10 m. All the soil samples were packed into polyethylene bags and brought back to the lab. The soil samples were air-dried in the laboratory. Subsequently, after stones and roots were removed from the samples, we grinded and sieved the samples to less than 2 mm before sending them for chemical analysis. Fuyang County is marked by a long history of mining activities and PTEs smelting activities, as well as important agricultural grain production regions [53]. Previous studies have revealed that mining activities, paper mills, cement factory, and metallurgic activities, as well as agricultural activities, are the main sources of PTEs in soil in Fuyang County [53,61,62].

<figure><img src="image_6.png" /><figcaption>Figure 1. Study area and sampling locations.</figcaption></figure>

### 2.2. Chemical Analysis

In this study, the soil pH was measured by a pH meter using a water/soil mixture with a ratio of 2.5:1. The pre-treatments included weighing portions of 0.2 g, microwave digestion (4 mL HNO₃, 2 mL HCl, and 2 mL HF) at a temperature of 120 °C for 45 min, liquid transition, and filtration [63]. Thereafter, the total concentration of all PTEs was determined by liquid chromatography-inductively coupled plasma mass spectrometry (LC-ICP-MS) (NexION 300X, PerkinElmer Inc., 940 Winter Street.
Waltham, MA 02451, USA). The total concentration of PTEs was measured twice in approximately 15% of soil samples to ensure the accuracy and repeatability of the method [63].

## 2.3. p-XRF

We used the “Soils Mode” of the portable Niton XL2 GOLDD XRF Analyzer (Thermo Fisher Scientific Inc., Billerica, MA, USA) to obtain p-XRF spectra information. All the soil samples were tested three times in the same place in the laboratory, with a scanning time of 1.5 min [63]. The Cubist algorithm was used to predict the total content of As, Cd, Cr, Cu, Ni, and Zn, based on the p-XRF spectra data. Detailed information on the processes of PTE estimation using p-XRF spectra data is provided in our previous studies [63]. The R² of the PTEs predicted by p-XRF spectra data is 0.50, 0.69, 0.96, 0.87, and 0.97 for As, Cd, Cu, Ni, and Zn, respectively.

## 2.4. Mapping Methods

### 2.4.1. IDW Method

IDW is one of the most popularly used deterministic interpolation methods in environmental and soil sciences [64,65]. IDW estimates the target variable value using the value from the nearby sampled sites. The weights assigned to the nearby sites are determined by the inverse of their distance from the estimated location. The IDW formula used in this study is as follows:

$$Z(x_0) = \frac{\sum_{i=1}^{N} \frac{x_i}{h_{ij}^2}}{\sum_{i=1}^{N} \frac{1}{h_{ij}^2}}, \quad (1)$$

where $Z(x_0)$ expresses the estimated value of the interpolated site; $N$ means the number of nearby samples used for interpolation. In this study, when we conducted the IDW interpolation, neighbor samples within the distance of 14.5 km were used for the estimation of PTEs in unvisited locations and the number of $N$ varied between 10 and 15; $x_i$ represents the value of $i$th samples; and $h_{ij}$ is the distance between the estimated point and sampled locations.

### 2.4.2. OK Method

OK is another popular method for exploring spatial patterns of environmental variables [66–70]. The experimental semivariogram was used to reveal the spatial dependence of the variables. The corresponding formula is as follows:

$$\gamma^*(h) = \frac{1}{2N(h)} \sum_{i=1}^{N(h)} [Z(x_i) - Z(x_i + h)]^2, \quad (2)$$

where $\gamma^*(h)$ represents the semi-variance; $N(h)$ represents the pairs of sample numbers separated by distance $h$; $Z(x_i)$ means the observed value at location $i$; and $Z(x_i + h)$ expresses the observed value at location $i + h$. Based on the experimental semivariogram, a model was then fitted to employ interpolation. Detailed information on the OK method was reported by Webster and Oliver (2007) [66].

### 2.4.3. EBK Method

EBK was developed by Krivoruchko (2012) [44] and aimed to resolve the drawbacks of classical geostatistical models. In EBK, the stochastic spatial process is represented locally as a stationary or non-stationary random field, and the corresponding parameters vary across space. EBK comprises two geostatistical models: the intrinsic random function kriging and the linear mixed model. The result of EBK is a robust non-stationary algorithm for spatial interpolation of geophysical corrections. When the data coverage is sufficient, the EBK algorithm helps extend local trends. When the data coverage
is poor, EBK integrates an a priori background mean. EBK differs from other kriging methods by including the error introduced by estimating the underlying semivariogram [71]. Detailed information on EBK can be found in Krivoruchko (2012) [44].

### 2.4.4. Model Averaging Using the Granger-Ramanathan Algorithm

The Granger-Ramanathan (GR) method was developed by Granger and Ramanathan in 1984 [72]. This method assumes that better predictions are produced by a linear combination of different individual predictions with the weight of each prediction determined by the ordinary least squares (OLS) algorithm. In this study, we fitted a linear regression model between the measured PTE content of the validation dataset and the PTEs predicted by the map obtained using different methods and data sources. The PTE<sub>GR</sub> outcomes produced using the GR method were obtained as follows:

$$PTE_{GR} = \sum_{i=1}^{N} (\alpha_{i} \cdot PTE_{i}) + \beta, \quad (3)$$

where $\alpha_i$ represents the regression coefficient; $PTE_i$ is the PTE prediction of the i-th PTE map ($n = 3$ in this study); and $\beta$ is the intercept. The values of $\alpha$ and $\beta$ are calculated using the OLS method. The sum of weights $\alpha_i$ does not need to be equal to 1.

### 2.5. Assessment of Model Performance

We divided the collected soil samples into a training dataset ($n = 195$) and an independent validation dataset ($n = 98$) with a ratio of 2:1. The performance of the map produced by the single method or model averaging was assessed through independent validation. Four indices, including the coefficient of determination ($R^2$), root mean square error (RMSE), relative root mean square error (RRMSE), and bias, were estimated to evaluate the accuracy of different maps.

### 2.6. Data Analysis

The summary statistics were conducted using R Studio [73]. The variograms of all PTEs were fitted by GS+ software (version 9.0; Gamma Software Design, 2008). Then, the parameters of variograms were recorded and used for ordinary interpolation performed using ArcGIS 10.3 (ESRI, ArcGIS 10.3, Redlands, CA, USA). Other geostatistical interpolations were also employed in ArcGIS 10.3 (ESRI, ArcGIS 10.3, Redlands, CA, USA). Maps of the PTEs content in Fuyang County were produced using IDW, OK, EBK, and GR model averaging. The main steps of data analysis, as used in this study, are presented in Figure 2.

<figure><img src="image_7.png" /><figcaption>Figure 2. Flowchart of methods used in this study.</figcaption></figure>
# 3. Results

## 3.1. Summary Statistics

In Table 1, we summarize the descriptive statistics of the As, Cd, Cu, Ni, and Zn contents measured by laboratory analysis (LC-ICP-MS) and p-XRF. The mean content of Cd was higher than the risk screening value for soils (0.3 mg kg⁻¹) provided in National Standards of China (GB15618–2018), while the mean content of other PTEs was lower than the corresponding risk screening value. This indicates that soil in the study area underwent Cd pollution. Generally, the PTEs measurements by p-XRF were consistent with those by LC-ICP-MS. However, the coefficient of variation (CV) of the PTE content determined by p-XRF was lower than that of LC-ICP-MS.

<table><caption>**Table 1.** Statistics of the potentially toxic elements (PTEs) content measured by LC-ICP-MS and p-XRF (mg/kg) (n = 293).</caption><thead><tr><th>Variable</th><th>Min</th><th>Median</th><th>Mean</th><th>Max</th><th>SD</th><th>CV (%)</th></tr></thead><tbody><tr><td>LC-ICP-MS As</td><td>0.20</td><td>10.65</td><td>12.45</td><td>37.37</td><td>6.51</td><td>52.29</td></tr><tr><td>p-XRF As</td><td>1.69</td><td>12.44</td><td>12.49</td><td>31.43</td><td>4.82</td><td>38.59</td></tr><tr><td>LC-ICP-MS Cd</td><td>0.09</td><td>0.33</td><td>0.41</td><td>2.35</td><td>0.31</td><td>75.61</td></tr><tr><td>p-XRF Cd</td><td>0.08</td><td>0.37</td><td>0.43</td><td>2.53</td><td>0.26</td><td>60.47</td></tr><tr><td>LC-ICP-MS Cu</td><td>7.79</td><td>26.31</td><td>31.34</td><td>177.22</td><td>21.80</td><td>69.56</td></tr><tr><td>p-XRF Cu</td><td>3.39</td><td>26.60</td><td>31.49</td><td>177.27</td><td>21.74</td><td>69.04</td></tr><tr><td>LC-ICP-MS Ni</td><td>4.38</td><td>20.63</td><td>20.79</td><td>51.97</td><td>8.76</td><td>42.14</td></tr><tr><td>p-XRF Ni</td><td>0.85</td><td>21.12</td><td>21.08</td><td>52.04</td><td>8.59</td><td>40.75</td></tr><tr><td>LC-ICP-MS Zn</td><td>58.52</td><td>107.60</td><td>127.43</td><td>793.94</td><td>72.30</td><td>56.74</td></tr><tr><td>p-XRF Zn</td><td>54.06</td><td>109.11</td><td>127.95</td><td>893.77</td><td>72.48</td><td>56.65</td></tr></tbody></table>

Note: SD means standard deviation; CV means the coefficient of variation.

## 3.2. Mapping of PTEs Using Data from Individual Sources

We fitted the semivariogram of different PTEs to analyze the spatial structure of different PTEs. The nugget/sill ratio defined by Cambardella et al. (1994) [74] was calculated to analyze the spatial dependence of different PTEs. As presented in Figure 3 and Table 2, the nugget/sill ratio values were 47.27%, 7.60%, 20.85%, 27.16%, and 18.27% for As, Cd, Cu, Ni, and Zn measured by the LC-ICP-MS method, respectively. Meanwhile, the nugget/sill ratio values were 49.98%, 21.05%, 20.90%, 37.48%, and 18.53% for As, Cd, Cu, Ni, and Zn measured by the p-XRF method, respectively. These values indicate strong spatial dependence for Cd, Cu, and Zn and moderate spatial dependence for As and Ni. With regard to the range, Cu has the largest spatial dependence range. Considering the background information of the study area and parameters of the semivariogram, the spatial dependence of Cu and Ni in Fuyang County is mainly controlled by intrinsic variation in soil characteristics such as texture and mineralogy, while for As, Cd, and Zn, it is controlled by both intrinsic variation and extrinsic variations such as industrial activities and mining activities, as well as agricultural actions [74].
<figure><img src="image_8.png" /><figcaption>Figure 3. Cont.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 3.** Semivariogram of PTEs content measured by LC-ICP-MS (left column) and p-XRF (right column).</figcaption></figure>

<table><caption>**Table 2.** Parameters of the semivariogram of different PTEs.</caption><thead><tr><th>Variable</th><th>Model Type</th><th>Nugget</th><th>Sill</th><th>Nugget/Sill (%)</th><th>Range (m)</th><th>R<sup>2</sup></th></tr></thead><tbody><tr><td>LC-ICP-MS As</td><td>Spherical</td><td>22.99</td><td>48.64</td><td>47.27</td><td>20200</td><td>0.897</td></tr><tr><td>LC-ICP-MS Cd</td><td>Spherical</td><td>0.01</td><td>0.172</td><td>0.76</td><td>20930</td><td>0.589</td></tr><tr><td>LC-ICP-MS Cu</td><td>Spherical</td><td>142.0</td><td>681.1</td><td>20.85</td><td>23830</td><td>0.640</td></tr><tr><td>LC-ICP-MS Ni</td><td>Spherical</td><td>20.80</td><td>76.59</td><td>27.16</td><td>16120</td><td>0.823</td></tr><tr><td>LC-ICP-MS Zn</td><td>Gaussian</td><td>1730.0</td><td>9470.0</td><td>18.27</td><td>13050</td><td>0.518</td></tr><tr><td>p-XRF As</td><td>Gaussian</td><td>12.91</td><td>25.83</td><td>49.98</td><td>11190</td><td>0.805</td></tr><tr><td>p-XRF Cd</td><td>Gaussian</td><td>0.024</td><td>0.124</td><td>21.05</td><td>11800</td><td>0.533</td></tr><tr><td>p-XRF Cu</td><td>Spherical</td><td>144.0</td><td>688.9</td><td>20.90</td><td>23910</td><td>0.622</td></tr><tr><td>p-XRF Ni</td><td>Spherical</td><td>28.70</td><td>76.58</td><td>37.48</td><td>18740</td><td>0.869</td></tr><tr><td>p-XRF Zn</td><td>Gaussian</td><td>1810.0</td><td>9770.0</td><td>18.53</td><td>13090</td><td>0.498</td></tr></tbody></table>

The map of the content of As, Cd, Cu, Ni, and Zn determined by LC-ICP-MS and p-XRF was outputted using the IDW, OK, and EBK methods (Figures 4–8). The results showed that the As values were lower than the risk screening values regulated by the Chinese National Standards (GB15618-2018) in most of the study area (Figure 4). The elevated As values in the soil were predominantly detected in the northwestern and central parts of Fuyang County. The low As values were concentrated in the southeast of Fuyang County. Many mines are located in the northwest of the study area, whereas its central part is represented by the flood plain—an important agricultural production area. Long-term application of chemical fertilizers and active mining may have contributed to a relatively high As
content in the soil in the northwestern and central parts of the survey region. Although the spatial pattern of As content was well illustrated in Figure 4b,d,f compared with Figure 4a,c,e, the As content would be overestimated in areas with a relatively low value of As, while underestimated in regions with a relatively high value of As in Figure 4b,d,f because the As measured by p-XRF has a narrower range compared with LC-ICP-MS with a larger minimum and smaller maximum of As content (Table 1).

<figure><img src="image_10.png" /><figcaption>**Figure 4.** Maps of As content measured by LC-ICP-MS (left column) and p-XRF (right column) using IDW (a,b), OK (c,d), and EBK (e,f) geostatistical methods.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 5.** Maps of Cd content measured by LC-ICP-MS (left column) and p-XRF (right column) using IDW (a,b), OK (c,d), and EBK (e,f) geostatistical methods.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 6.** Maps of Cu content measured by LC-ICP-MS (left column) and p-XRF (right column) using IDW (a,b), OK (c,d), and EBK (e,f) geostatistical methods.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 7.** Maps of Ni content measured by LC-ICP-MS (left column) and p-XRF (right column) using IDW (a,b), OK (c,d), and EBK (e,f) geostatistical methods.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 8.** Maps of Zn content measured by LC-ICP-MS (left column) and p-XRF (right column) using IDW (a,b), OK (c,d), and EBK (e,f) geostatistical methods.</figcaption></figure>
The content of Cd showed a spatial pattern similar to As (Figure 5). Some regions revealed pollution by Cd, where the Cd content in the soil was higher than the risk screening value regulated by the Chinese National Standards (GB15618-2018). High Cd values in the soil were primarily located in the northwestern, north, and central parts of Fuyang County. We attribute Cd pollution to the long-term application of chemical fertilizers, mining activity, and industrial activity in these regions, especially considering that agricultural and mining activities have been previously identified as important sources of Cd and As in Fuyang County [75]. The relative smaller spatial dependence range of Cd measured by p-XRF could explain the appearance of artificial lines apparently crossing the study area presented in Figure 5b,d,f.

The results showed that the Cu content in most areas of the survey region was lower than the risk screening value regulated by the Chinese National Standards (GB15618-2018) (Figure 6). The elevated Cu content was only found in the southeast of Fuyang County, downstream of the Fuchun River. We argue that Cu contained in chemical fertilizers and pesticides applied in other regions of the study area could have been transported by the Fuchun River and then accumulated in the soil of this region by irrigation, which could also be one of the main sources of As and Cd accumulated in soil in the southeast of Fuyang County (Figures 4 and 5).

The Ni content in most parts of the survey region was lower than the risk screening value regulated by the Chinese National Standards (GB15618-2018) (Figure 7). A high Ni value was located only in the northwest and north of Fuyang County. Waste discharge from mining activity was likely the primary source of Ni accumulated in the soil of this area.

The Zn content in most parts of the survey region was lower than the risk screening value regulated by the Chinese National Standards (GB15618-2018) (Figure 8). High values of Zn were mainly found in the north, northwest, and southeast of Fuyang County. Zn in these regions was likely mainly sourced from mining activities.

As presented in Figures 4–8, the maps of the PTEs generated by different interpolation methods and measurement methods showed very similar spatial patterns. However, the maps generated by the IDW method revealed more local variations of the PTEs, while the maps produced by OK provided the smoothest spatial patterns of PTEs. The local high values of PTEs in the center of the study area are mainly rooted in industrial activities, waste emission of factories, and mining activities with a small scale.

It is worth noting that there are some artificial lines in Figures 4d, 6d and 7d which were produced by the OK method, and that dramatic changes of PTE content occur in the areas located on the different sides of these artificial lines. Adjusting the parameters such as the search window and number of neighbor samples may resolve this issue and could obtain smoother maps of soil PTEs. In addition, improving the sample size is also expected to avoid the appearance of these artificial lines.

### 3.3. Spatial Modeling of PTEs Using Model Averaging

In this study, model averaging was performed to produce the final PTE map using the six maps constructed by a combination of LC-ICP-MS measurements with IDW, OK, and EBK and p-XRF measurements with IDW, OK, and EBK. As presented in Figure 9, in general, the maps produced by model averaging method show similar spatial patterns to the maps produced by IDW, OK and EBK.
<figure><img src="image_15.png" /><figcaption>Figure 9. PTE maps produced by model averaging (right column) algorithms.</figcaption></figure>

### 3.4. Comparison of Performance of Different Mapping Methods

In this study, we validated all the models by comparing their estimations with lab measurements (LC-ICP-MS measurements) using the independent validation dataset (Table 3, Figures S1–S4). We estimated the accuracy of the measuring and interpolation methods using their different combinations (Table 3). The absolute values of the RMSE of the different interpolation methods follow a decreasing order: Zn > Cu > Ni > As > Cd (Table 3). With regard to the RRMSE, the interpolation models clearly have a lower RRMSE value for Ni and Zn compared with As, Cd, and Cu. In addition, the LC-ICP-MS spatial estimation outperformed p-XRF for As, Ni, and Zn with a smaller value of the RMSE and RRMSE (Table 3) [76,77].
<table><caption>**Table 3.** Performance of different models using combined data evaluated by an independent validation dataset.</caption><thead><tr><th>Algorithms</th><th>Method</th><th>Element</th><th>R<sup>2</sup></th><th>RMSE (mg kg<sup>-1</sup>)</th><th>RRMSE (%) <sup>a</sup></th><th>Bias (mg kg<sup>-1</sup>)</th></tr></thead><tbody><tr><td>IDW</td><td>LC-ICP-MS</td><td>As</td><td>0.38</td><td>4.71</td><td>37.83</td><td>0.21</td></tr><tr><td>IDW</td><td>p-XRF</td><td>As</td><td>0.28</td><td>4.96</td><td>39.71</td><td>0.71</td></tr><tr><td>OK</td><td>LC-ICP-MS</td><td>As</td><td>0.34</td><td>4.90</td><td>39.36</td><td>0.32</td></tr><tr><td>OK</td><td>p-XRF</td><td>As</td><td>0.17</td><td>5.44</td><td>43.55</td><td>0.79</td></tr><tr><td>EBK</td><td>LC-ICP-MS</td><td>As</td><td>0.44</td><td>4.37</td><td>35.10</td><td>0.34</td></tr><tr><td>EBK</td><td>p-XRF</td><td>As</td><td>0.29</td><td>4.93</td><td>39.47</td><td>0.57</td></tr><tr><td>IDW</td><td>LC-ICP-MS</td><td>Cd</td><td>0.65</td><td>0.18</td><td>43.90</td><td>-0.01</td></tr><tr><td>IDW</td><td>p-XRF</td><td>Cd</td><td>0.59</td><td>0.19</td><td>44.19</td><td>0.01</td></tr><tr><td>OK</td><td>LC-ICP-MS</td><td>Cd</td><td>0.65</td><td>0.19</td><td>46.34</td><td>-0.04</td></tr><tr><td>OK</td><td>p-XRF</td><td>Cd</td><td>0.61</td><td>0.19</td><td>44.19</td><td>0.002</td></tr><tr><td>EBK</td><td>LC-ICP-MS</td><td>Cd</td><td>0.58</td><td>0.19</td><td>46.34</td><td>-0.005</td></tr><tr><td>EBK</td><td>p-XRF</td><td>Cd</td><td>0.58</td><td>0.19</td><td>44.19</td><td>0.01</td></tr><tr><td>IDW</td><td>LC-ICP-MS</td><td>Cu</td><td>0.60</td><td>13.15</td><td>41.96</td><td>1.08</td></tr><tr><td>IDW</td><td>p-XRF</td><td>Cu</td><td>0.60</td><td>13.18</td><td>41.85</td><td>1.27</td></tr><tr><td>OK</td><td>LC-ICP-MS</td><td>Cu</td><td>0.60</td><td>13.53</td><td>43.17</td><td>1.00</td></tr><tr><td>OK</td><td>p-XRF</td><td>Cu</td><td>0.59</td><td>13.40</td><td>42.55</td><td>1.22</td></tr><tr><td>EBK</td><td>LC-ICP-MS</td><td>Cu</td><td>0.58</td><td>13.47</td><td>42.98</td><td>1.67</td></tr><tr><td>EBK</td><td>p-XRF</td><td>Cu</td><td>0.58</td><td>13.45</td><td>42.71</td><td>1.68</td></tr><tr><td>IDW</td><td>LC-ICP-MS</td><td>Ni</td><td>0.61</td><td>5.54</td><td>26.65</td><td>0.58</td></tr><tr><td>IDW</td><td>p-XRF</td><td>Ni</td><td>0.53</td><td>6.22</td><td>29.51</td><td>1.20</td></tr><tr><td>OK</td><td>LC-ICP-MS</td><td>Ni</td><td>0.63</td><td>5.52</td><td>26.55</td><td>0.84</td></tr><tr><td>OK</td><td>p-XRF</td><td>Ni</td><td>0.52</td><td>6.46</td><td>30.65</td><td>1.82</td></tr><tr><td>EBK</td><td>LC-ICP-MS</td><td>Ni</td><td>0.64</td><td>5.32</td><td>25.59</td><td>0.53</td></tr><tr><td>EBK</td><td>p-XRF</td><td>Ni</td><td>0.55</td><td>6.04</td><td>28.65</td><td>0.99</td></tr><tr><td>IDW</td><td>LC-ICP-MS</td><td>Zn</td><td>0.67</td><td>35.39</td><td>27.77</td><td>10.51</td></tr><tr><td>IDW</td><td>p-XRF</td><td>Zn</td><td>0.64</td><td>36.60</td><td>28.60</td><td>10.36</td></tr><tr><td>OK</td><td>LC-ICP-MS</td><td>Zn</td><td>0.59</td><td>35.80</td><td>28.09</td><td>2.73</td></tr><tr><td>OK</td><td>p-XRF</td><td>Zn</td><td>0.59</td><td>35.71</td><td>27.91</td><td>2.70</td></tr><tr><td>EBK</td><td>LC-ICP-MS</td><td>Zn</td><td>0.64</td><td>35.79</td><td>28.09</td><td>10.13</td></tr><tr><td>EBK</td><td>p-XRF</td><td>Zn</td><td>0.63</td><td>36.27</td><td>28.35</td><td>10.41</td></tr></tbody></table>

Note: a: the RRMSE means relative root mean square error. It is equal to the ratio of the mean of the square root of residuals squared with the mean content of PTEs [76,77].

As listed in Table 3, generally for As, Cu, Ni, and Zn, all the methods overestimated PTEs content when mapping the spatial distribution of PTEs with a positive value of bias. The results presented in Figures S1–S4 indicate a general trend of overestimation of low PTEs concentration and underestimation of high levels of PTEs concentration when mapping all the elements using different methods. In terms of Cd, the results produced by different methods using LC-ICP-MS measurements slightly underestimated the value of Cd, while using p-XRF measurements slightly overestimated the value of Cd.

Overall, our results revealed that the accuracy of the LC-ICP-MS spatial estimation was higher than that of p-XRF for all PTEs, as also reported by others related studies [78,79]. However, the accuracy of the p-XRF spatial prediction of Cd, Cu, Zn, and Ni contents was very similar to that of the laboratory measurements. The best spatial prediction of different PTEs was obtained from various combinations of measuring methods and algorithms, highlighting the necessity of using model averaging to improve the prediction accuracy even further.

As for As in the soil, the most accurate prediction results were produced by a combination of the LC-ICP-MS measurements and the EBK algorithm, with an R² of 0.44, RMSE of 4.37 mg kg⁻¹, RRMSE of 35.10%, and a bias of 0.34 mg kg⁻¹. The best prediction for Cd was obtained by the combination of the LC-ICP-MS measurements and the IDW method, with an R² of 0.65, RMSE of 0.18 mg kg⁻¹,
RRMSE of 43.90%, and a bias of -0.01 mg kg⁻¹. With regard to Cu, the most accurate estimation was acquired by the combination of the LC-ICP-MS measurements and the OK method, with an R² of 0.60, RMSE of 13.53 mg kg⁻¹, RRMSE of 43.17%, and a bias of 1.00 mg kg⁻¹. For Ni, the combination of the LC-ICP-MS measurements and the EBK method provided the best prediction results, with an R² of 0.64, RMSE of 5.32 mg kg⁻¹, RRMSE of 25.59%, and a bias of 0.53 mg kg⁻¹. As for Zn, using the LC-ICP-MS measurements and the IDW method produced the best estimation accuracy with an R² of 0.67, RMSE of 35.39 mg kg⁻¹, RRMSE of 27.77%, and a bias of 10.51 mg kg⁻¹.

### 3.5. PTE Mapping Using Model Averaging

Herein, we applied the GR algorithm to map the spatial distribution of PTEs based on the individual maps. These maps were produced by the IDW, OK, and EBK methods using the LC-ICP-MS and p-XRF measurements. Then, we compared the results produced by model averaging with optimized results produced by individual interpolation methods and measurement methods.

As shown in Table 4, the PTEs map produced by model averaging showed minor improvement compared with the map with the highest accuracy produced based on individual data sources (i.e., LC-ICP-MS or p-XRF measurements) using the IDW, OK, or EBK algorithms. Compared with the IDW, OK, and EBK algorithms, the maps of all PTEs produced by model averaging had higher accuracy and were characterized by a higher R², and a lower value of the RMSE and RRMSE. Especially, compared with individual interpolation methods and measurement methods, model averaging essentially reduced the bias value of the prediction results. Thus, model averaging can be used to improve the estimation of the spatial distribution of soil PTEs, especially of As.

<table><caption>**Table 4.** Accuracy parameters of PTEs maps obtained by an optimized individual method and model averaging.</caption><thead><tr><th>Algorithms</th><th>Method</th><th>R<sup>2</sup></th><th>RMSE (mg kg<sup>-1</sup>)</th><th>RRMSE (%) <sup>a</sup></th><th>Bias (mg kg<sup>-1</sup>)</th></tr></thead><tbody><tr><td rowspan="3">Individual method</td><td>As</td><td>0.44</td><td>4.37</td><td>35.10</td><td>0.34</td></tr><tr><td>Model averaging</td><td>As</td><td>0.50</td><td>4.08</td><td>32.77</td><td>-0.005</td></tr><tr><td>Improvement (%)</td><td>As</td><td>13.64</td><td>-6.64</td><td>-6.64</td><td>-98.53</td></tr><tr><td rowspan="3">Individual method</td><td>Cd</td><td>0.65</td><td>0.18</td><td>43.90</td><td>-0.01</td></tr><tr><td>Model averaging</td><td>Cd</td><td>0.68</td><td>0.17</td><td>41.46</td><td>-0.0002</td></tr><tr><td>Improvement (%)</td><td>Cd</td><td>4.62</td><td>-5.56</td><td>-5.56</td><td>-98.00</td></tr><tr><td rowspan="3">Individual method</td><td>Cu</td><td>0.60</td><td>13.53</td><td>43.17</td><td>1.00</td></tr><tr><td>Model averaging</td><td>Cu</td><td>0.62</td><td>12.64</td><td>40.33</td><td>0.02</td></tr><tr><td>Improvement (%)</td><td>Cu</td><td>3.33</td><td>-6.58</td><td>-6.58</td><td>-98.00</td></tr><tr><td rowspan="3">Individual method</td><td>Ni</td><td>0.64</td><td>5.32</td><td>25.59</td><td>0.53</td></tr><tr><td>Model averaging</td><td>Ni</td><td>0.67</td><td>5.09</td><td>24.48</td><td>0.002</td></tr><tr><td>Improvement (%)</td><td>Ni</td><td>4.69</td><td>-4.32</td><td>-4.34</td><td>-99.62</td></tr><tr><td rowspan="3">Individual method</td><td>Zn</td><td>0.67</td><td>35.39</td><td>27.77</td><td>10.51</td></tr><tr><td>Model averaging</td><td>Zn</td><td>0.71</td><td>29.19</td><td>22.91</td><td>-0.09</td></tr><tr><td>Improvement (%)</td><td>Zn</td><td>5.63</td><td>-17.52</td><td>-17.50</td><td>-99.14</td></tr></tbody></table>

Note: <sup>a</sup> the RRMSE means relative root mean square error. It is equal to the ratio of the mean of the square root of residuals squared with the mean content of PTEs [76,77].

## 4. Discussion

### 4.1. Spatial Prediction Accuracy Using Different Interpolation Methods

In this study, the prediction accuracy of different interpolation methods for estimating the PTE content measured by the LC-ICP-MS and p-XRF methods is summarized in Tables 3 and 4. The interpolation methods used herein had a good performance for all PTEs, except for As. The bad performance of the interpolation method for As may mainly be attributed to two sources. On the one hand, as we have described in Section 2.3, the content of soil As measured by p-XRF clearly
has lower accuracy compared with others elements, which then also leads to lower accuracy for the interpolation results. On the other hand, as presented in Figure 10, the histograms of the content of As measured by both LC-ICP-MS and p-XRF clearly deviated from a normal distribution even when a log10() transformation was employed. However, a normal distribution of data was usually expected to employ kriging interpolation. This could also explain the relative bad performance of the interpolation results for As.

<figure><img src="image_16.png" /><figcaption>**Figure 10.** Histogram of the original data of content of As measured by LC-ICP-MS (a) and p-XRF (b), and log10(-)-transformed content of As measured by LC-ICP-MS (c) and p-XRF (d).</figcaption></figure>

EBK provides the highest spatial prediction accuracy for As distribution, IDW provides the highest spatial prediction accuracy for Cd distribution, OK provides the highest spatial prediction accuracy for Cu distribution, IDW provides the highest spatial prediction accuracy for Ni distribution, and ENK provides the highest spatial prediction accuracy for Zn distribution. The EBK method did not perform better than IDW and OK, according to the study results. These findings confirmed that spatial interpolation algorithms are data-specific and variable-specific [64], and there exists no interpolation method that always shows the best performance for all the variables. However, different methods used herein did not output significantly different estimations of PTEs content in the study area. We attribute this fact to the relatively high sampling density in this study. As stated elsewhere, when the sampling density is high, most interpolation methods tend to output similar results [80]. Dirks et al. (1998) [81] found that when using a high-density dataset, kriging methods did not show any major improvement in prediction accuracy compared with simple methods (e.g., IDW, nearest neighbors).

With regard to the PTE spatial distribution, different interpolation methods produced very similar results. However, local variations in PTEs became apparent in the map produced by the IDW algorithm. The map produced by the OK method smoothed out most of the local details. These findings were consistent with those of previous studies [82–84]. In the PTE estimation results produced by the EBK algorithm, we detected a more dramatic PTE variation in the study area. This is because, unlike in OK, in EBK, the stochastic spatial process is represented locally as a stationary or non-stationary random field, so the prediction varies more strongly [71,85]. In addition, the sample size is another important factor to affect the interpolation accuracy [86–88]. In this study, we collected a total of 293 soil samples. Although this sample size has proved to be able to capture the spatial variation of soil heavy metals at
the county level [89–93], it is still far from an extensive survey and more soil samples are expected to improve the interpolation accuracy. Therefore, in further study, more soil samples need to be collected to analyze the effects of the sampling density on the interpolation accuracy.

## 4.2. p-XRF Measurements as an Alternative Means for PTE Mapping

The summary statistics extracted from the p-XRF measurements of As, Cd, Cu, Ni, and Zn in the soil were only slightly worse compared to the laboratory analysis (LC-ICP-MS) results, proving the potential of p-XRF to be used as an alternative method for measuring As, Cd, Cu, Ni, and Zn contents in soil, which could reduce the time and economy cost for the measurement of PTEs content in the laboratory. With regard to the spatial prediction, the prediction accuracy produced by different methods using p-XRF was close to that of the traditional laboratory analysis, except for As. The application of the p-XRF method reduces the time, labor, and financial costs associated with soil sampling, soil sample pre-processing, and laboratory analysis [94]. As such, it enables us to collect and analyze more soil samples with a higher sampling density, which is crucial for conducting soil surveys at large spatial scales. Another potential application of p-XRF is to help understand the general status of soil PTE pollution and provide prior knowledge for designing more efficient soil sampling schemes.

## 4.3. Potential of Model Averaging for PTE Mapping

Our study results demonstrated the superiority of using model averaging (GR) to improve the mapping of all soil PTEs (Figure 8, Table 3). As shown in Table 4, using the model averaging method reduced the RRMSE value of all PTEs. Especially, the model averaging method obviously reduced the bias value which confirms the potential of model averaging for improved mapping of soil PTEs. The degree of improvement varied for different elements. The most significant improvement when employing model averaging was shown for As, where the $R^2$ of the prediction of model averaging was significantly improved, while the RMSE, RRMSE, and bias decreased compared with other interpolation methods. In addition, the prediction results of model averaging retained more local variations in soil PTEs compared with the OK method, while some of the extremely abrupt variations in soil PTEs were smoothed out compared with EBK, making the results seem closer to the actual situation. In this study, we only tested the feasibility of GR for mapping soil PTEs using model averaging. In our future work, some other model averaging approaches, such as variance weighted, Bayesian model averaging, piecewise linear decision tree, and weighted model averaging, will be tested to map soil PTEs with larger expected improvement [95–97].

## 4.4. Method Limitations and Potential Improvements

Herein, we proved the potential of p-XRF to serve as an alternative method for laboratory analysis and the potential of model averaging to map soil PTEs. However, some limitations of the application of p-XRF and model averaging still need to be overcome.

First, PTEs from the obtained soil samples were measured using p-XRF after soil sample pre-processing (i.e., we air-dried the samples in the laboratory, removed stones and roots, and, finally, grounded and sieved the samples to be less than 2 mm in size). There are still some limitations to realize in situ and on-the-go measurements of soil PTEs because the soil spectra are extremely susceptible to interference from the external environment (e.g., soil moisture and particle size) and reduce the prediction accuracy [98,99]. Therefore, p-XRF should have a better possibility of application if we could reduce the negative effects of the environmental factors and then realize in situ and even on-the-go measurements [100]. Advanced pre-processing algorithms, such as external parameter orthogonalization and direct standardization, can be used to deal with this issue [101,102].

Second, as reported by Chen et al. (2020) [96], the sample size substantially affects the accuracy of model averaging methods. Therefore, it is necessary to compare the accuracy of spatial interpolation algorithms with that of the model averaging method using datasets with different sizes or sampling densities. This could help achieve a better trade-off between prediction accuracy and cost.
After that, the sampling design and spatial distribution of the samples also have essential effects on the prediction results. The sampling design could affect the reliability of the variogram [64,103]. The sampling design also has a significant influence on the prediction accuracy with irregular designs preferred to regular ones [88]. Sample spacing must relate to the scale or scales of the variation of the primary variable in a region, otherwise samples might be too sparsely spaced to identify a correlation and could result in a pure nugget [66]. Spatial clustering of samples affects the accuracy of the estimations, but the effects vary among different interpolation methods [64,66]. Data collected at a suitable range of separations to capture changes in the scales of the variation deserve higher prediction accuracy [104]. Therefore, more attention should be paid to the soil sampling design and spatial distribution of soil samples in the future.

Finally, the performance of the model averaging method may vary among different survey regions and soil properties. Thus, extended applications are needed to confirm the feasibility of the model averaging approach for mapping soil and other PTEs in different areas and at variable spatial scales. Moreover, other types of model averaging methods, such as variance weighted, Bayesian model averaging, piecewise linear decision tree, and weighted model averaging, could also be tested to improve the accuracy of PTE mapping.

## 5. Conclusions

This study is the first to compare the mapping of soil PTEs using individual geostatistical methods (i.e., IDW, OK, and EBK) and a model averaging approach based on the laboratory and p-XRF measurements. The feasibility of p-XRF as an alternative method of laboratory analysis for mapping soil PTEs was confirmed. The main study conclusions are as follows:

1. The average content of all PTEs was lower than the corresponding risk screening value regulated by the Chinese National Standards (GB15618-2018), except for Cd.

2. The EBK method showed no clear advantage over IDW and OK, according to the study results. Moreover, no method regarded herein consistently outperformed other methods when interpolating different PTEs.

3. The p-XRF spatial prediction accuracy of Cd, Cu, Ni, and Zn was similar to that of the laboratory measurements, confirming that p-XRF can be used as a reliable alternative to traditional laboratory analysis for mapping Cd, Cu, Ni, and Zn in the soil.

4. The model averaging method improved the mapping accuracy of all soil PTEs regarded herein compared with individual geostatistical algorithms (i.e., IDW, OK, and EBK). The improvement in mapping accuracy was the clearest for As.

Our study proved that p-XRF and model averaging provide a cost-effective, reliable, and highly accurate alternative to map soil PTEs, which could essentially reduce the time and economy cost for soil sampling and monitoring. This advantage is more obvious when soil PTEs pollution monitoring and remediation are carried out at a large spatial scale. This study provides an implication for more efficient monitoring and remediation of soil PTE pollution.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2072-4292/12/22/3775/s1, Figure S1: Scatter plot of spatial prediction values using IDW Vs PTEs content measured by LC-ICP-MS method, Figure S2: Scatter plot of spatial prediction values using OK Vs PTEs content measured by LC-ICP-MS method, Figure S3: Scatter plot of spatial prediction values using EBK Vs PTEs content measured by LC-ICP-MS method, Figure S4: Scatter plot of spatial prediction values using model averaging Vs PTEs content measured by LC-ICP-MS method.

**Author Contributions:** F.X.: conceptualization, methodology, software, writing—review and editing, resources. B.H.: conceptualization, methodology, data curation, writing—review and editing. Y.Z.: supervision, resources. W.J.: writing—review and editing, software. S.C.: writing—review and editing. D.X.: data curation. Z.S.: conceptualization, supervision, writing—review and editing, resources. All authors have read and agreed to the published version of the manuscript.
**Funding:** This research was funded by the National Key Research and Development Program of China (2018YFC1800105) and the Scientific Research and Development Fund of Zhejiang A&F University (Talent Startup Project)(2020FR060).

**Acknowledgments:** This research was funded by the National Key Research and Development Program of China (2018YFC1800105). We also acknowledge the support from the Talent Introduction Project of Zhejiang A&F University. We acknowledge the support received by Bifeng Hu from the China Scholarship Council (201706320317) for three years of Ph.D. study in the French National Institute for Agriculture, Food, and Environment (INRAE) and Orléans University in France. We also thank the editor and reviewers for their valuable and constructive comments on our manuscript which helped us essentially improve the manuscript.

**Conflicts of Interest:** The authors have declared that no competing interests exist. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript or in the decision to publish the results.

## References

1. Lacarce, E.; Saby, N.P.A.; Martin, M.P.; Marchant, B.P.; Boulonne, L.; Meersmans, J.; Jolivet, C.; Bispo, A.; Arrouays, D. Mapping soil Pb stocks and availability in mainland France combining regression trees with robust geostatistics. *Geoderma* **2012**, *170*, 359–368. [CrossRef]

2. Hu, B.F.; Xue, J.; Zhou, Y.; Shao, S.; Fu, Z.; Li, Y.; Chen, S.C.; Qi, L.; Shi, Z. Modelling bioaccumulation of heavy metals in soil-crop ecosystems and identifying its controlling factors using machine learning. *Environ. Pollut.* **2020**, 114308. [CrossRef] [PubMed]

3. Jia, X.L.; Hu, B.F.; Marchant, B.P.; Zhou, L.Q.; Shi, Z.; Zhu, Y.W. A methodological framework for identifying potential sources of soil heavy metal pollution based on machine learning: A case study in the Yangtze Delta, China. *Environ. Pollut.* **2019**, *250*, 601–609. [CrossRef] [PubMed]

4. Jiang, Y.F.; Chen, S.C.; Hu, B.F.; Zhou, Y.; Liang, Z.Z.; Jia, X.L.; Huang, M.X.; Wei, J.; Shi, Z. A comprehensive framework for assessing the impact of potential agricultural pollution on grain security and human health in economically developed areas. *Environ. Pollut.* **2020**, 114653. [CrossRef]

5. Lequy, E.; Saby, N.P.A.; Ilyin, I.; Bourin, A.; Sauvage, S.; Leblond, S. Spatial analysis of trace elements in a moss bio-monitoring data over France by accounting for source, protocol and environmental parameters. *Sci. Total Environ.* **2017**, *590*, 602–610. [CrossRef]

6. Fu, T.T.; Zhao, R.Y.; Hu, B.F.; Jia, X.L.; Wang, Z.G.; Zhou, L.Q.; Huang, M.X.; Li, Y.; Shi, Z. Novel framework for modelling the cadmium balance and accumulation in farmland soil in Zhejiang Province, East China: Sensitivity analysis, parameter optimisation, and forecast for 2050. *J. Clean Prod.* **2020**, 123674. [CrossRef]

7. Nriagu, J.O. A history of global metal pollution. *Science* **2017**, *272*, 223. [CrossRef]

8. Zhao, F.J.; Ma, Y.; Zhu, Y.G.; Tang, Z.; McGrath, S.P. Soil contamination in China: Current status and mitigation strategies. *Environ. Sci. Technol.* **2015**, *49*, 750–759. [CrossRef]

9. Yang, Q.; Li, Z.; Lu, X.; Duan, Q.; Huang, L.; Bi, J. A review of soil heavy metal pollution from industrial and agricultural regions in China: Pollution and risk assessment. *Sci. Total Environ.* **2018**, *642*, 690–700. [CrossRef]

10. Shi, A.; Shao, Y.; Zhao, K.; Fu, W. Long-term effect of E-waste dismantling activities on the heavy metals pollution in paddy soils of southeastern China. *Sci. Total Environ.* **2020**, *705*, 135971. [CrossRef]

11. Kumar, V.; Sharma, A.; Kaur, P.; Sidhu, G.P.S.; Bali, A.S.; Bhardwaj, R.; Thukral, A.K.; Cerda, A. Pollution assessment of heavy metals in soils of India and ecological risk assessment: A state-of-the-art. *Chemosphere* **2019**, *216*, 449–462. [CrossRef] [PubMed]

12. Dogra, N.; Sharma, M.; Sharma, A.; Keshavarzi, A.; Minakshi; Bhardwaj, R.; Thukral, A.K.; Kumar, V. Pollution assessment and spatial distribution of roadside agricultural soils: A case study from India. *Int. J. Environ. Health Res.* **2020**, *30*, 146–159. [CrossRef] [PubMed]

13. Adimalla, N. Heavy metals contamination in urban surface soils of Medak province, India, and its risk assessment and spatial distribution. *Environ. Geochem. Health* **2020**, *42*, 59–75. [CrossRef] [PubMed]

14. Doabi, S.A.; Karami, M.; Afyuni, M. Heavy metal pollution assessment in agricultural soils of Kermanshah province, Iran. *Environ. Earth Sci.* **2019**, *78*, 70. [CrossRef]

15. Solgi, E.; Abbasitabr, H.; Shayesteh, K.; Mortazavi, S. Mapping and risk assessment of heavy metals in agricultural soils of the Siakh Darengoun Region, Shiraz, Iran. *Environ. Resour. Res.* **2020**, *8*, 83–95.
16. Taati, A.; Salehi, M.H.; Mohammadi, J.; Mohajer, R.; Díez, S. Pollution assessment and spatial distribution of trace elements in soils of Arak industrial area, Iran: Implications for human health. *Environ. Res.* **2020**, *187*, 109577. [CrossRef]

17. Hafeez, F.; Zafar, N.; Nazir, R.; Javeed, H.M.R.; Rizwan, M.; Asad, S.A.; Iqbal, A. Assessment of flood-induced changes in soil heavy metal and nutrient status in Rajanpur, Pakistan. *Environ. Monit. Assess.* **2019**, *191*, 234. [CrossRef]

18. Rehman, I.U.; Ishaq, M.; Ali, L.; Muhammad, S.; Din, I.U.; Yaseen, M.; Ullah, H. Potentially toxic elements' occurrence and risk assessment through water and soil of Chitral urban environment, Pakistan: A case study. *Environ. Geochem. Health* **2020**, *1–14*. [CrossRef]

19. Hashmi, M.Z.; Kanwal, A.; Pongpiachan, S.; Su, X.; Nasim, W.; Saand, M.A.; Mubeen, M.; Akram, R.; Wang, S.; Ahmed, Z. Arsenic distribution and metabolism genes abundance in Paddy soils from Punjab and Sindh provinces, Pakistan. *Arab. J. Geosci.* **2020**, *13*, 1–10. [CrossRef]

20. Nascimento, V.X.; Barros, A.; Azevedo, J.A.; Miranda, P.R.; da Costa, J.G. Bioavailability of heavy metals in mangrove soil in Alagoas, Brazil. *Biosci. J.* **2019**, *818–825*. [CrossRef]

21. Davila, R.B.; Fontes, M.P.F.; Pacheco, A.A.; da Silva Ferreira, M. Heavy metals in iron ore tailings and floodplain soils affected by the Samarco dam collapse in Brazil. *Sci. Total Environ.* **2020**, *709*, 136151. [CrossRef] [PubMed]

22. de Menezes, M.D.; Bispo, F.H.A.; Faria, W.M.; Gonçalves, M.G.M.; Curi, N.; Guilherme, L.R.G. Modeling arsenic content in Brazilian soils: What is relevant? *Sci. Total Environ.* **2020**, *712*, 136511. [CrossRef] [PubMed]

23. Nahar, K.; Ali, M.M.; Khanom, A.; Alam, M.K.; Azad, M.A.K.; Rahman, M.M. Levels of heavy metal concentrations and their effect on net nitrification rates and nitrifying archaea/bacteria in paddy soils of Bangladesh. *Appl. Soil Ecol.* **2020**, *156*, 103697. [CrossRef]

24. Akber, M.A.; Rahman, M.A.; Islam, M.A.; Islam, M.A. Potential ecological risk of metal pollution in lead smelter-contaminated agricultural soils in Khulna, Bangladesh. *Environ. Monit. Assess.* **2019**, *191*, 351. [CrossRef] [PubMed]

25. Islam, M.S.; Ahmed, M.K.; Al-Mamun, M.H.; Eaton, D.W. Human and ecological risks of metals in soils under different land-use types in an urban environment of Bangladesh. *Pedosphere* **2020**, *30*, 201–213. [CrossRef]

26. Hu, B.F.; Jia, X.L.; Hu, J.; Xu, D.Y.; Xia, F.; Li, Y. Assessment of heavy metal pollution and health risks in the soil-plant-human system in the Yangtze River Delta, China. *Int. J. Environ. Res. Public Health* **2017**, *14*, 1042. [CrossRef]

27. Jia, X.L.; Fu, T.T.; Hu, B.F.; Shi, Z.; Zhou, L.Q.; Zhu, Y.W. Identification of the potential risk areas for soil heavy metal pollution based on the source-sink theory. *J. Hazard. Mater.* **2020**, *393*, 122424. [CrossRef]

28. Hu, B.F.; Shao, S.; Ni, H.; Fu, Z.Y.; Hu, L.S.; Zhou, Y.; Min, X.X.; She, S.F.; Chen, S.C.; Huang, M.X.; et al. Current status, spatial features, health risks, and potential driving factors of soil heavy metal pollution in China at province level. *Environ. Pollut.* **2020**, 114961. [CrossRef]

29. Saby, N.P.A.; Thioulouse, J.; Jolivet, C.C.; Ratié, C.; Boulonne, L.; Bispo, A.; Arrouays, D. Multivariate analysis of the spatial patterns of 8 trace elements using the French soil monitoring network data. *Sci. Total Environ.* **2009**, *407*, 5644–5652. [CrossRef]

30. Liu, Y.; Fei, X.F.; Zhang, Z.; Li, Y.; Tang, J.; Xiao, R. Identifying the sources and spatial patterns of potentially toxic trace elements (PTEs) in Shanghai suburb soils using global and local regression models. *Environ. Pollut.* **2020**, *264*, 114171. [CrossRef]

31. Satarug, S.; Baker, J.R.; Urbenjapol, S.; Haswell-Elkins, M.; Reilly, P.E.; Williams, D.J.; Moore, M.R. A global perspective on cadmium pollution and toxicity in non-occupationally exposed population. *Toxicol. Lett.* **2003**, *137*, 65–83. [CrossRef]

32. Smith, A.H.; Goycolea, M.; Haque, R.; Biggs, M.L. Marked increase in bladder and lung cancer mortality in a region of northern Chile due to arsenic in water. *Am. J. Epidemiol.* **1998**, *147*, 660–669. [CrossRef] [PubMed]

33. Fei, X.F.; Christakos, G.; Xiao, R.; Ren, Z.; Liu, Y.; Lv, X.N. Improved heavy metal mapping and pollution source apportionment in Shanghai City soils using auxiliary information. *Sci. Total Environ.* **2019**, *661*, 168–177. [CrossRef] [PubMed]

34. Cheng, J.L.; Shi, Z.; Zhu, Y.W. Assessment and mapping of environmental quality in agricultural soils of Zhejiang Province, China. *J. Environ. Sci.* **2007**, *19*, 50–54. [CrossRef]
35. Shao, S.; Hu, B.F.; Fu, Z.Y.; Wang, J.Y.; Lou, G.; Zhou, Y.; Jin, B.; Li, Y.; Shi, Z. Source identification and apportionment of trace elements in soils in the Yangtze River Delta, China. *Int. J. Environ. Res. Public Health* **2018**, *15*, 1240. [CrossRef]

36. Hu, B.F.; Zhou, Y.; Jiang, Y.F.; Ji, W.J.; Fu, Z.Y.; Shao, S.; Li, S.; Huang, M.X.; Zhou, L.Q.; Shi, Z. Spatio-temporal variation and source changes of potentially toxic elements in soil on a typical plain of the Yangtze River Delta, China (2002–2012). *J. Environ. Manag.* **2020**, *271*, 110943. [CrossRef]

37. Yang, Y.; Yang, X.; He, M.; Christakos, G. Beyond mere pollution source identification: Determination of land covers emitting soil heavy metals by combining PCA/APCS, GeoDetector and GIS analysis. *Catena* **2020**, *185*, 104297. [CrossRef]

38. Xiang, M.T.; Li, Y.; Yang, J.; Li, Y.; Li, F.; Hu, B.F.; Cao, Y. Assessment of Heavy Metal Pollution in Soil and Classification of Pollution Risk Management and Control Zones in the Industrial Developed City. *Environ. Manag.* **2020**. [CrossRef]

39. Liu, S.; Zhao, H.; Wu, K.; Zhang, Z.; Hou, Y.; Chen, T.; Jin, Q. Evaluation of heavy metal distribution characteristics of agricultural soil–rice system in a high geological background area according to the influence index of comprehensive quality (IICQ). *Environ. Sci. Pollut. Res.* **2020**, *27*, 20920–20933. [CrossRef]

40. Hu, B.F.; Shao, S.; Fu, Z.Y.; Li, Y.; Ni, H.; Chen, S.C.; Zhou, Y.; Jin, B.; Shi, Z. Identifying heavy metal pollution hot spots in soil-rice systems: A case study in South of Yangtze River Delta, China. *Sci. Total Environ.* **2019**, *658*, 614–625. [CrossRef]

41. Mukhopadhyay, S.; Chakraborty, S.; Bhadoria, P.B.S.; Li, B.; Weindorf, D.C. Assessment of heavy metal and soil organic carbon by portable X-ray fluorescence spectrometry and NixPro™ sensor in landfill soils of India. *Geoderma. Reg.* **2020**, *20*, e00249. [CrossRef]

42. Wang, Y.; Duan, X.; Wang, L. Spatial distribution and source analysis of heavy metals in soils influenced by industrial enterprise distribution: Case study in Jiangsu Province. *Sci. Total Environ.* **2020**, *710*, 134953. [CrossRef] [PubMed]

43. Hu, B.F.; Wang, J.Y.; Jin, B.; Li, Y.; Shi, Z. Assessment of the potential health risks of heavy metals in soils in a coastal industrial region of the Yangtze River Delta. *Environ. Sci. Pollut. Res.* **2017**, *24*, 19816–19826. [CrossRef] [PubMed]

44. Krivoruchko, K. *Empirical Bayesian Kriging*; Esri: Redlands, CA, USA, 2012.

45. Onovo, A.; Atobatele, A.; Kalaiwo, A.; Obanubi, C.; James, E.; Gado, P.; Odezugo, G.; Ogundehin, D.; Magaji, D.; Russell, M. Using Supervised Machine Learning and Empirical Bayesian Kriging to reveal Correlates and Patterns of COVID-19 Disease outbreak in sub-Saharan Africa: Exploratory Data Analysis. *SSRN Electr. J.* **2020**. [CrossRef]

46. Gupta, A.; Kamble, T.; Machiwal, D. Comparison of ordinary and Bayesian kriging techniques in depicting rainfall variability in arid and semi-arid regions of north-west India. *Environ. Earth Sci.* **2017**, *76*, 512. [CrossRef]

47. Giustini, F.; Ciotoli, G.; Rinaldini, A.; Ruggiero, L.; Voltaggio, M. Mapping the geogenic radon potential and radon risk by using Empirical Bayesian Kriging regression: A case study from a volcanic area of central Italy. *Sci. Total Environ.* **2019**, *661*, 449–464. [CrossRef]

48. Ji, W.J.; Li, S.; Chen, S.C.; Shi, Z.; Rossel, R.A.V.; Mouazen, A.M. Prediction of soil attributes using the Chinese soil spectral library and standardized spectra recorded at field conditions. *Soil Till. Res.* **2016**, *155*, 492–500. [CrossRef]

49. Li, S.; Shi, Z.; Chen, S.C.; Ji, W.J.; Zhou, L.Q.; Yu, W.; Webster, R. In situ measurements of organic carbon in soil profiles using vis–NIR spectroscopy on the Qinghai–Tibet Plateau. *Environ. Sci. Technol.* **2015**, *49*, 4980–4987. [CrossRef]

50. Caporale, A.G.; Adamo, P.; Capozzi, F.; Langella, G.; Terribile, F.; Vingiani, S. Monitoring metal pollution in soils using portable-XRF and conventional laboratory-based techniques: Evaluation of the performance and limitations according to metal properties and sources. *Sci. Total Environ.* **2018**, *643*, 516–526. [CrossRef]

51. Liu, G.; Liu, B.; Yang, L.; Hu, W.; Qu, M.; Lu, F.; Huang, B. Using pXRF to assess the accumulation, sources, and potential ecological risk of potentially toxic elements in soil under two greenhouse vegetable production systems in North China. *Environ. Sci. Pollut. Res.* **2020**, *27*, 11105–11115. [CrossRef]

52. Peinado, F.M.; Ruano, S.M.; González, M.B.; Molina, C.E. A rapid field procedure for screening trace elements in polluted soil using portable X-ray fluorescence (XRF). *Geoderma* **2010**, *159*, 76–82.
53. Hu, B.F.; Chen, S.C.; Hu, J.; Xia, F.; Xu, J.F.; Li, Y.; Shi, Z. Application of portable XRF and VNIR sensors for rapid assessment of soil heavy metal pollution. *PLoS ONE* **2017**, *12*, e0172438.

54. Wan, M.X.; Hu, W.Y.; Qu, M.K.; Tian, K.; Zhang, H.D.; Wang, Y.; Huang, B. Application of arc emission spectrometry and portable X-ray fluorescence spectrometry to rapid risk assessment of heavy metals in agricultural soils. *Ecol. Indic.* **2019**, *101*, 583–594.

55. Declercq, Y.; Delbecque, N.; De Grave, J.; De Smedt, P.; Finke, P.; Mouazen, A.M.; Nawar, S.; Vandenberghe, D.; Van Meirvenne, M.; Verdoodt, A. A Comprehensive Study of Three Different Portable XRF Scanners to Assess the Soil Geochemistry of An Extensive Sample Dataset. *Remote Sens.* **2019**, *11*, 2490.

56. Kim, H.R.; Kim, K.H.; Yu, S.; Moniruzzaman, M.; Hwang, S.I.; Lee, G.T.; Yun, S.T. Better assessment of the distribution of As and Pb in soils in a former smelting area, using ordinary co-kriging and sequential Gaussian co-simulation of portable X-ray fluorescence (XRF) and ICP-AES data. *Geoderma* **2019**, *341*, 26–38.

57. Xia, F.; Hu, B.F.; Shao, S.; Xu, D.; Zhou, Y.; Zhou, Y.; Zhou, Y.; Huang, M.X.; Li, Y.; Chen, S.C.; et al. Improvement of Spatial Modeling of Cr, Pb, Cd, As and Ni in Soil Based on Portable X-ray Fluorescence (PXRF) and Geostatistics: A Case Study in East China. *Int. J. Environ. Res. Public Health* **2019**, *16*, 2694.

58. Arsenault, R.; Gatien, P.; Renaud, B.; Brissette, F.; Martel, J.L. A comparative analysis of 9 multi-model averaging approaches in hydrological continuous streamflow simulation. *J. Hydrol.* **2015**, *529*, 754–767.

59. Malone, B.P.; Minasny, B.; Odgers, N.P.; McBratney, A.B. Using model averaging to combine soil property rasters from legacy soil maps and from point data. *Geoderma* **2014**, *232*, 34–44.

60. O'Rourke, S.M.; Stockmann, U.; Holden, N.M.; McBratney, A.B.; Minasny, B. An assessment of model averaging to improve predictive power of portable vis-NIR and XRF for the determination of agronomic soil properties. *Geoderma* **2016**, *279*, 31–44.

61. Yin, X.; Yao, C.; Song, J.; Li, Z.; Zhang, C.; Qian, W.; Bi, D.; Li, C.X.; Teng, Y.; Wan, H.D.; et al. Mercury contamination in vicinity of secondary copper smelters in Fuyang, Zhejiang Province, China: Levels and contamination in topsoils. *Environ. Pollut.* **2009**, *157*, 1787–1793.

62. Zhang, X.Y.; Lin, F.F.; Wong, M.T.; Feng, X.L.; Wang, K. Identification of soil heavy metal sources from anthropogenic activities and pollution assessment of Fuyang County, China. *Environ. Monit. Assess.* **2009**, *154*, 439. [CrossRef] [PubMed]

63. Xu, D.Y.; Chen, S.C.; Xu, H.Y.; Wang, N.; Zhou, Y.; Shi, Z. Data fusion for the measurement of potentially toxic elements in soil using portable spectrometers. *Environ. Pollut.* **2020**, *263*, 114649.

64. Li, J.; Heap, A.D. Spatial interpolation methods applied in the environmental sciences: A review. *Environ. Modell. Softw.* **2014**, *53*, 173–189. [CrossRef]

65. Hu, B.F.; Wang, J.Y.; Fu, T.T.; Li, Y.; Shi, Z. Application of Spatial Analysis on Soil Heavy Metal Contamination: A Review. *Chin. J. Soil. Sci.* **2017**, *48*, 1014–1024. (In Chinese)

66. Webster, R.; Oliver, M.A. Geostatistics for Environmental Scientists, 2nd ed.; Wiley: Hoboken, NJ, USA, 2008.

67. Hu, B.F.; Zhao, R.Y.; Chen, S.C.; Zhou, Y.; Jin, B.; Li, Y.; Shi, Z. Heavy metal pollution delineation based on uncertainty in a coastal industrial city in the Yangtze River Delta, China. *Int. J. Environ. Res. Public Health* **2018**, *15*, 710. [CrossRef]

68. Goovaerts, P. Geostatistics for Natural Resources Evaluation; Oxford University Press: Oxford, UK, 1997.

69. Zhou, Y.; Jia, Z.; Wang, J.; Chen, L.; Zou, M.; Li, Y.; Zhou, S. Heavy metal distribution, relationship and prediction in a wheat-rice rotation system. *Geoderma* **2019**, *354*, 113886.

70. Yang, L.; Xu, H.; Yu, S. Estimating PM2.5 concentrations in Yangtze River Delta region of China using random forest model and the Top-of-Atmosphere reflectance. *J. Environ. Manag.* **2020**, *272*, 111061.

71. Gribov, A.; Krivoruchko, K. Empirical Bayesian kriging implementation and usage. *Sci. Total Environ.* **2020**, *722*, 137290. [CrossRef]

72. Granger, C.W.; Ramanathan, R. Improved methods of combining forecasts. *J. Forecast.* **1984**, *3*, 197–204.

73. Team, R. RStudio: Integrated development for R; RStudio, Inc.: Boston, MA, USA, 2015; Volume 42, p. 14. Available online: http://www.rstudio.com (accessed on 6 November 2020).

74. Cambardella, C.A.; Moorman, T.B.; Novak, J.M.; Parkin, T.B.; Karlen, D.L.; Turco, R.F.; Konopka, A.E. Field-scale variability of soil properties in central Iowa soils. *Soil Sci. Soc. Am. J.* **1994**, *58*, 1501–1511. [CrossRef]

75. Zhang, Y.Y.; Ma, J.; Wei, H.Y.; Shi, T.R. Heavy metals in typical farmland soils of Zhejiang Province: Levels, sources and ecological risks. *Ecol. Environ. Sci.* **2019**, *28*, 1233–1241. (In Chinese)
76. Despotovic, M.; Nedic, V.; Despotovic, D.; Cvetanovic, S. Evaluation of empirical models for predicting monthly mean horizontal diffuse solar radiation. *Renew. Sustain. Energ. Rev.* **2016**, *56*, 246–260.

77. Eyduran, E.; Zaborski, D.; Waheed, A.; Celik, S.; Karadas, K.; Grzesiak, W. Comparison of the predictive capabilities of several data mining algorithms and multiple linear regression in the prediction of body weight by means of body measurements in the indigenous Beetal goat of Pakistan. *Pak. J. Zool.* **2017**, *49*, 57–265.
[CrossRef]

78. Lee, H.; Choi, Y.; Suh, J.; Lee, S.H. Mapping copper and lead concentrations at abandoned mine areas using element analysis data from ICP-AES and portable XRF instruments: A comparative study. *Int. J. Environ. Res. Public Health* **2016**, *13*, 384. [CrossRef]

79. Horta, A.; Azevedo, L.; Neves, J.; Soares, A.; Pozza, L. Integrating portable X-ray fluorescence (pXRF) measurement uncertainty for accurate soil contamination mapping. *Geoderma* **2021**, *382*, 114712. [CrossRef]

80. Burrough, P.A.; McDonnell, R.A. Creating continuous surfaces from point data. In *Principles of Geographic Information Systems*; Oxford University Press: Oxford, UK, 1998.

81. Dirks, K.N.; Hay, J.E.; Stow, C.D.; Harris, D. High-resolution studies of rainfall on Norfolk Island Part II: Interpolation of rainfall data. *J. Hydrol.* **1998**, *208*, 187–193.

82. Xie, B.; Jia, X.; Qin, Z.; Zhao, C.; Shao, M.A. Comparison of interpolation methods for soil moisture prediction on China's Loess Plateau. *Vadose Zone J.* **2020**, *19*, e20025.

83. Metahni, S.; Coudert, L.; Gloaguen, E.; Guemiza, K.; Mercier, G.; Blais, J.F. Comparison of different interpolation methods and sequential Gaussian simulation to estimate volumes of soil contaminated by As, Cr, Cu, PCP and dioxins/furans. *Environ. Pollut.* **2019**, *252*, 409–419. [CrossRef]

84. Shen, Q.; Wang, Y.; Wang, X.; Liu, X.; Zhang, X.; Zhang, S. Comparing interpolation methods to predict soil total phosphorus in the Mollisol area of Northeast China. *Catena* **2019**, *174*, 59–72. [CrossRef]

85. Krivoruchko, K.; Gribov, A. Evaluation of empirical Bayesian kriging. *Spat. Stat.-NETH* **2019**, *32*, 100368.
[CrossRef]

86. Yan, B.; Pan, Y.; Zhao, C. Spatial variability and reasonable sampling number of regional soil heavy metals. *Trans. Chin. Soc. Agric. Eng.* **2008**, *24*, 260–264.

87. Zhengping, X.; Xingwei, Y.; Xiangsuo, D.; Xian, L. Spatial Variability of Soil Nutrient and Reasonable Sampling Number. *Trans. Chin. Soc. Agric. Eng.* **2002**, *18*, 6–9.

88. Li, J.; Heap, A.D. A review of comparative studies of spatial interpolation methods in environmental sciences: Performance and impact factors. *Ecol. Inform.* **2011**, *6*, 228–241. [CrossRef]

89. Song, D.; Zhuang, D.; Jiang, D.; Fu, J.; Wang, Q. Integrated health risk assessment of heavy metals in Suxian County, South China. *Int. J. Environ. Res. Public Health* **2015**, *12*, 7100–7117. [CrossRef] [PubMed]

90. Pan, L.; Ma, J.; Hu, Y.; Su, B.; Fang, G.; Wang, Y.; Wang, L.; Xiang, B. Assessments of levels, potential ecological risk, and human health risk of heavy metals in the soils from a typical county in Shanxi Province, China. *Environ. Sci. Pollut. Res.* **2016**, *23*, 19330–19340. [CrossRef] [PubMed]

91. Ungureanu, T.; Iancu, G.O.; Pintilei, M.; Chicos, M.M. Spatial distribution and geochemistry of heavy metals in soils: A case study from the NE area of Vaslui county, Romania. *J. Geochem. Explor.* **2017**, *176*, 20–32.
[CrossRef]

92. Lü, J.; Jiao, W.B.; Qiu, H.Y.; Chen, B.; Huang, X.X.; Kang, B. Origin and spatial distribution of heavy metals and carcinogenic risk assessment in mining areas at You'xi County southeast China. *Geoderma* **2018**, *310*, 99–106. [CrossRef]

93. Wang, S.; Cai, L.M.; Wen, H.H.; Luo, J.; Wang, Q.S.; Liu, X. Spatial distribution and source apportionment of heavy metals in soil from a typical county-level city of Guangdong Province, China. *Sci. Total Environ.* **2019**, *655*, 92–101. [CrossRef]

94. Lemière, B. A review of pXRF (field portable X-ray fluorescence) applications for applied geochemistry. *J. Geochem. Explor.* **2018**, *188*, 350–363. [CrossRef]

95. Liu, C.; Maheu, J.M. Forecasting realized volatility: A Bayesian model-averaging approach. *J. Appl. Economet.* **2009**, *24*, 709–733. [CrossRef]

96. Chen, S.C.; Mulder, V.L.; Heuvelink, G.B.M.; Poggio, L.; Caubet, M.; Dobarco, M.R.; Walter, C.; Arrouays, D. Model averaging for mapping topsoil organic carbon in France. *Geoderma* **2020**, *366*, 114237. [CrossRef]

97. Zhou, Y.; Xue, J.; Chen, S.C.; Zhou, Y.; Liang, Z.Z.; Wang, N.; Shi, Z. Fine-Resolution Mapping of Soil Total Nitrogen across China Based on Weighted Model Averaging. *Remote Sens.* **2020**, *12*, 85. [CrossRef]
98. Ji, W.J.; Rossel, R.A.V.; Shi, Z. Improved estimates of organic carbon using proximally sensed vis-NIR spectra corrected by piecewise direct standardization. *Eur. J. Soil. Sci.* **2015**, *66*, 670–678. [CrossRef]

99. Minasny, B.; McBratney, A.B.; Bellon-Maurel, V.; Roger, J.M.; Gobrecht, A.; Ferrand, L.; Joalland, S. Removing the effect of soil moisture from NIR diffuse reflectance spectra for the prediction of soil organic carbon. *Geoderma* **2011**, *167*, 118–124. [CrossRef]

100. Weindorf, D.C.; Paulette, L.; Man, T. In-situ assessment of metal contamination via portable X-ray fluorescence spectroscopy: Zlatna, Romania. *Environ. Pollut* **2013**, *182*, 92–100. [CrossRef]

101. Ji, W.J.; Viscarra Rossel, R.; Shi, Z. Accounting for the effects of water and the environment on proximally sensed vis-NIR spectra and their calibrations. *Eur. J. Soil. Sci.* **2015**, *66*, 555–565. [CrossRef]

102. Nawar, S.; Cipullo, S.; Douglas, R.K.; Coulon, F.; Mouazen, A.M. The applicability of spectroscopy methods for estimating potentially toxic elements in soils: State-of-the-art and future trends. *Appl. Spectrosc. Rev.* **2020**, *55*, 525–557. [CrossRef]

103. Brus, D.J.; de Gruijter, J.J. Estimating of non-ergodic variograms and their sampling variance by design-based sampling strategies. *Math. Geol.* **1994**, *26*, 437–454. [CrossRef]

104. Laslett, G.M. Kriging and splines: An empirical comparison of their predictive performance in some applications. *J. Am. Stat. Assoc.* **1994**, *89*, 391–400. [CrossRef]

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_17.png" /></figure>

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).