1. Introduction
Soil is usually regarded as both a sink and source for heavy metals (HMs) and other pollutants [1]. An excess level of heavy metals in soil poses a potential threat to the health of plants, animals, and even humans via the food chain because of their persistence, high toxicity, non-biodegradable, and accumulative properties [2,3,4,5,6]. Consequently, soil heavy metal pollution has become a major concern in environmental management. Numerous studies have focused on HMs’ spatial distribution and pollution assessment during the last few decades [2,7,8,9,10,11]. In these studies, many pollution indices have been applied to evaluate soil contamination assessment. Among these indices, the geo-accumulation index [12], pollution load index [13], enrichment factor index [14], risk index [15], and Nemero index are the most commonly used methods. Additionally, the contour map is also a powerful tool for visualizing the contamination degree and grasping the spatial characteristics of HMs’ concentration.
Generally, researchers usually rely on interpolation methods, including kriging [16,17] and inverse distance [18,19] approaches, to estimate the HM concentrations and pollution indices of the unsampled area and therefore construct continuous surface distribution maps [20,21]. In most existing studies, SurferTM (Golden Software, Inc., Golden, CO, USA) and ArcGIS software (ESRI, Sacramento, CA, USA) can provide a way to realize interpolation algorithms in the form of a contour map [22,23,24,25]. Compared with ArcGIS, Surfer software is more efficient, being faster in data processing and plotting. However, when using this software, several time-consuming and labor-intensive problems may still arise due to a flood of geochemical data caused by multiple elements or too many sampling sites. For example, to assess soil contamination, researchers attempt to choose a proper method to calculate the pollution index several times in Microsoft ExcelTM, before generating contour maps via Surfer software. Then, they repeatedly plot contour maps based on the above calculation results or spatial distribution contour maps based on sampling data sets for each single heavy metal. Therefore, it is of great significance to eliminate the above routines and tedious manual operations. Gong developed the software EGAPI, which was built with the above pollution indices methods, representing one of the few attempts to improve the efficiency of pollution index calculation [26]. Nevertheless, this is far from enough, as there are still other problems that need to be solved, such as the need for a user-friendly interface, the lack of batch plotting HM contamination assessment or spatial distribution contour maps, and the requirement of the batch output for publication-quality contour maps.
To solve the problems outlined above and develop a more scientific, faster in data handling and contour map drawing software, a combination of Visual Basic (VB) programming language and Surfer automation technology would be a good solution. The Surfer automation technique allows scripter or other high-level programming languages, such as VB or C++, to call surfer’s drawing function and is widely used in soil moisture [27], exploration geochemistry [28], meteorology [29], and hydrology [30]. Similarly, an integrated software HMCA-Contour written in VB and based on Surfer automation was designed to eliminate routine and tedious operations involving datasets of HMs in soil topics. The main functions of HMCA-Contour include the following: (1) Batch calculating pollution indices for HMs; (2) batch plotting maps of the spatial distribution and soil contamination assessment of HMs; and (3) template setting for contour maps. This paper presents the details of HMCA-Contour through a case study.
2. Methodology
Several common contamination assessment indexes of soil HMs, including the geo-accumulation index (Igeo), pollution load index (PLI), enrichment factor (EF), potential ecological risk index (RI), single pollution index (SI), and Nemero index (NI), were selected and integrated into HMCA-Contour software.
2.1. Geo-Accumulation Index
The geo-accumulation index (Igeo ) was put forward by Müller [12]. Originally used in river bottom sediments, it can also be applied to the assessment of heavy metal contamination in soil [31].Igeois computed by the following formula:
Igeo=log2(Cn∕1.5Bn),
whereCn is the measured content of element n in soil,Bnis the geochemical background value of element n in soil, and the constant 1.5 is the background matrix correction factor due to lithogenic effects in different regions.
According to Müller [12], the geo-accumulation index consists of 7 grades or classes (Table 1), whereby the highest class 6 is an open class comprising all values of the index higher than class 5 and may reflect a 100-fold enrichment above the geochemical background value.
2.2. Pollution Load Index
The pollution load index (PLI), proposed by Tomlinson et al. [13], was employed to determine the magnitude of heavy metal contamination in the sediment. PLI was determined as the nth root of the product of the n contamination factor (CF). It is expressed as follows:
PLI=(CF1×CF2×CF3×…×CFn)1∕n,
where n is the number of metals and CF is the heavy metal concentration in the soil/background value of the metal. Interpretation PLI > 1 indicates polluted, whereas PLI < 1 indicates no pollution.
2.3. Enrichment Factor Index
The enrichment factor (EF) was first proposed by Zoller et al. [32] based on a study of the enrichment degree and source of elements in the atmosphere over Antarctica. Originally applied to the atmosphere and seawater, EF was gradually introduced to environment media in continental environments, such as soils [14], stream sediments [33], and river suspended solids [34] and it has been widely used in environmental geochemical heavy metal pollution assessment. The EF is defined mathematically as
EF=(Ci/Cref)sample(Ci/Cref)background,
where the concentration of the heavy element i isCiandCref refers to the concentration of the reference element. The most common reference elements are Fe, Al, Ti, Mn, and Sc [34,35,36,37,38,39]. According to Sutherland [35], EF is divided into five categories: Minimal pollution (EF < 2); moderate pollution (EF = 2–5); significant pollution (EF = 5–20); strong pollution (EF = 20–40); and extreme pollution (EF > 40).
2.4. Ecological Risk Index
The potential ecological risk index (RI) method is a sediment evaluation method proposed by Swedish scientist Hakanson [15] in 1980. It is a relatively fast, simple, and standard method for evaluating the sediment pollution degree and potential ecological risk index from the perspective of sedimentology.
The potential ecological risk index (RI) was suggested by Hakanson [15] and conducted to evaluate the impact of heavy metals in soil or sediment [40,41,42]. Three aspects were considered in this method: The toxicity of various heavy metals; the synergistic action of multiple elements; and the sensitivity of various heavy metal contamination in soils or sediments. Three Equations (4)–(6) for computing RI are as follows:
Efi=Eri×Efi,
Cfi=Csi∕Bni,
RI=∑i=1n Eri,
whereTrirefers to the toxic-response factor for a given element i, with values of 10, 2, 5, 30, 5, 1, and 40 for As, Cr, Cu, Cd, Pb, Zn, and Hg, respectively;Efirepresents an ecological risk factor for heavy metal i;Cfi,Csi, andBniare the contamination factor, the concentration in the sediment, and the background reference value for element i, respectively; and the potential ecological index (RI) is defined as the sum of risk factors.
According to Hakanson [15], the Er and RI can be divided into four grades and these are described in Table 2.
2.5. Single and Nemero Pollution Index
The single pollution index (Pi) quantifies the pollution of one single heavy metal [43]. It was defined as the ratio betweenCiandSi:
Pi=Ci∕Si,
whereCi identifies the concentration of each heavy metal in each sample andSi refers to the soil standard of each heavy metal.
The Nemerow pollution index (PN ) provided by Nemerow [44] was initially employed to water quality assessment and gradually extended to soil contamination [45,46,47]. This method takes both the mean and maximum of single pollution levels of all the evaluation factors into consideration. It was adopted and calculated as
PN=Piave2+Pimax22,
wherePiaveandPimax represent the maximum value and average value of the single factor pollution index in each sample, respectively. Pi and PN were respectively divided into four and five grades [48,49] and their criteria are listed in Table 3.
3. Software Description 3.1. Software Operation Environment Requirements
HMCA-Contour is temporarily only provided for mainstream Microsoft Windows (since 7), and it requires Excel 2007 or later and Surfer 13.0 to be installed on the user’s computer, which enables Visual Basic to invoke the automation method provided by Surfer software (version 13.0). The download link can provide access to the installation file from the website (https://github.com/HMCAHome/HMCA-Contour (accessed on 10 October 2020)). After downloading the installation file, the user can then run HMCA-Contour setup.exe and follow the installation directions on the screen. Once the installation procedure is completed, the user can run executable file HMCA-Contour. The flow chart of the operation of this software can be seen in Figure 1.
3.2. Data Processing Feature
3.2.1. Data Input
HMCA-Contour accepts two types of input files, i.e., .xlsx or .xls files. With the sheet named “Sheet1”, the first two columns of the template input file must be X and Y coordinates, respectively, and followed by a series of elements (Figure 2) that can be in any order. The user can enter the “Data Processing” interface by clicking the “Data processing” tab in the toolbar of the HMCA-Contour main interface (Figure 3a). In this interface, three subpanels are given, including “File Location”, “Data Extraction”, and “Data Processing”.
In the “File Location” subframe, the user loads the data file by clicking the “Open” button to choose the input file. After selection of the input file, this program automatically extracts the first row data in the left box of the “Data Extraction” frame, and the user is capable of choosing target items by clicking direction buttons (buttons “>”, “<”, “>>”, and “<<” represent add, reduce, select all, and deselect all, respectively) (Figure 3b).
3.2.2. Data Processing
In the “Data Processing” module, six HMCA methods, as mentioned above, are provided in this software (Figure 3b). The user can select the target method to work with data. It is noted that EF requires additional parameters, including the reference element (Cref element) and background value (Cref Background) (Figure 3b). Once the evaluation method is selected, the user can click the “Calculate” button to start calculating and is prompted to input the background value for each item according to the default format in a message box (Figure 3(C1)). After inputting the background value in the default form for each element, the user clicks the “OK” button, and the whole data handling process is then completed in seconds (Figure 3(C2)). The resulting filename is in the form of “source file name”_”selected method”.xls and saved in the “Output File” path directory, which is the same as the “Source File” path in default. It is noted that the units of heavy metal concentrations are required to be uniform (mg/kg).
3.3. Batch Plotting Function
3.3.1. Data Input
In order to enter the plotting interface (Figure 3d), the user simply clicks the “Next Step” button in the “Data Processing” interface or the plotting tab in the main interface of this software. In this function, MS Excel XLSX, XLS, and Comma Separated Value (CSV) files can be used as original data files for HMCA-Contour. Like “Data Extraction” in “Data Processing” capability, the program can extract the data automatically from the original file in the left textbox when the user clicks the “open” button. At the same time, the grid geometry parameters (Minimum, Maximum, Spacing, and Nodes) are calculated automatically. Additionally, the default options for plotting are loaded, for example, “File name”, “Output Path”, “Image Type”, and “Image Quality” in the output filename description frame, and the “Font Face”, “Font Size”, “Symbol Size”, “Symbol Type” in the “Postmap” frame. The program provides four kinds of image quality, from low to high. Moreover, the common image types (.jpg, .bmp, and .tif) are also provided, enabling the user to choose them by clicking on the drop-down menu. A Surfer project file (.srf) is generated at the same time, which provides the user with a way to easily modify the image.
3.3.2. Plotting
After reading data from the source file, the user may choose the target items from the left text box to the right text box in the “Choosing Items” subpanel. The buttons for selecting target elements are the same as those mentioned above (Figure 3d). It is noted that two items (x, y) must be left in the left text box. The user selects a gridding algorithm from the drop-down menu of “Gridding Method”. Additionally, two other kinds of data files are available in this function—the boundary data (.BLN) and the geographic data (.dat, .txt, .xls, and .xlsx)—with the former one being employed for boundary determination and the latter one being employed for the geographic location. Both “blank” and “post map” files are loaded by clicking the “open” button. In the post map frame, the symbol size and type, and the symbol label font face, can be set with the corresponding drop button. Plotting capabilities are based on Surfer automation. The user does not need to worry about the running process, whose flow chart can be seen in Figure 4.
3.3.3. Gridding Method Descriptions
Surfer software provides twelve algorithms built in this program [50], including Inverse Distance to a power, kriging, Minimum Curvature, Modified Shepard’s Method, Natural Neighbor, Nearest Neighbor, Polynomial Regression, Radial Basis Function, Triangulation with Linear Interpolation, Moving Average, Data Metrics, and Local Polynomial. According to the different characteristics of various types of data and different requirements, the user can scientifically choose the right interpolation through the HMCA-Contour interface to call Surfer software to grid discrete data without worrying about the gridding process. Since these interpolation algorithms are default in Surfer, the contour maps produced by HMCA-Contour and Surfer software can lack model accuracy in comparison to Arcgis under certain circumstances.
3.4. Template Setting
The user may enter the plotting setting interface by either clicking the “Next” button in the “Data processing” interface or clicking the function “Template Setting” tab in the main interface of HMCA-Contour. The Template Setting interface has five tabs: “Color scale”; “Axis Ticks”; “Axis Title”; “Axis Labels”; and “Map Title & Scale bar & North Arrow” (Figure 5). This function is designed to build a publication-ready template for a series of thematic maps. The details of the tabs are given below.
3.4.1. Color Scale
The program provides three group parameters to customize a color scale (Figure 5a,b). Firstly, the user can type a title for the color scale that consists of a prefix, element, and suffix. The “prefix” and “suffix” text boxes can be activated by their front check boxes. When activated, the “title” text can be formed by elements and activated text boxes. Secondly, two forms of color scale bars are offered: The simple one is a system-based preset, which provides several kinds of color scale bars in the “Preset” drop-down menu, whilst the advanced one requires the user to click the “more” button (Figure 5a). As soon as the “more” button is clicked, there is a pop window through which the user is capable of adjusting a desired color scale example based on “nums”, “zlevel”, and “lcolor” parameters (Figure 5b). The user then clicks the “ok” or “return” button. Note that the user must choose one method to show the color scale bar. Lastly, the color scale bar label is set as the axis label.
3.4.2. Axis Ticks
Four kinds of axis ticks are provided, including bottom, left, right, and top ticks. Each of them has the same parameters (Figure 5c). Take the bottom axis tick as an example, where the user can set parameters in terms of five aspects: (1) A major or minor tick type drop-down menu that enables the user to choose the desired tick type from a list containing “none”, “in”, “out”, and “both”; (2) a major or minor tick length spin box, which allows the user to use an up or down arrow button to respectively increase or decrease the number value of the tick length in the text field by a 0.01 increment; (3) a major or minor gird line drop-down menu that offers two choices (“True and False”) for the user to decide whether or not to display it; (4) a minor per major spin button that enables the user to set the number of minor ticks between each two major ticks; and (5) an axis line width spin button makes it possible for the user to adjust the axis line width by clicking the up or down arrow button. The same approach of adjusting parameters can be applied for the remaining kinds of axis ticks.
3.4.3. Axis Titles
The user is capable of customizing the axis title properties (Figure 5d). For example, if the user wants to add a bottom axis title, the user can type text in the text box, and then change the angle of the title by clicking the up or down spin button. Additionally, the user can set the title text font, size, vertical and horizontal distance to the axis, bold, and italic by clicking their corresponding drop-down menus to make a selection. Similarly, the other three axis title types can be set in the same way.
3.4.4. Axis Labels
Four check boxes allow the user to decide which axis labels need to be displayed on a map (Figure 5e). Once a check box is clicked, the user can further adjust the corresponding axis label parameters. The user clicks the up or down spin button to set a proper distance value between the axis line and labels. Moreover, by clicking the drop-down buttons, the font face, size, bold, and italic and angle of axis label can be adjusted via their drop list.
3.4.5. Map Title, Scale Bar, and North Arrow
Four sections are divided (Figure 5f). Among them, the method of setting parameters in the map title section occurs in the same way as in the color scale title section. With respect to the scale bar component, the user needs to adjust the scale bar label and appearance separately. A scale bar is a tool that helps the user with the reading of this scale when a scale rule is not available, whose style can be approached from three aspects: The spacing between cycles (Cyclespacing); the increment of the label value (Increment); and the number of cycles (Numcycles). All of these parameters can be set from their corresponding drop-down lists. Moreover, the user is able to change the angle of the scale bar. The north arrow section enables a north arrow to be set up on the map display. The user can use the “symbol type” drop-down list to choose the style, and then specify the location of its center position from the “X position” and “Y position” drop-down lists.
4. Case Study 4.1. Study Area and Data
To test the applicability of HMCA-Contour, we examined the soil geochemical data taken from the Gerqin Mine area [51]. This mining area (83°23′00′′E~83°27′00′′E, 32°47′00′′ N~32°50′00′′ N) is located in Wuma Township, Gaize County, Tibet. It is located about 100 km northwest of Gaize county, with a height of 4700–5500 m, relative height of 100–700 m, and area of 1700 km2. The climate of this area is cold and dry, with a long sunshine time, a large temperature difference between the day and night, and annual precipitation of about 190 mm. There is a large amount of evaporation in the mining area.
The soil in the study area is mainly composed of Quaternary alluvial diluvium and weathering products of the bedrock parent material. From east to west, there are eight mineral deposits (occurrences), including Bolong, Duobuza, Rongna, Naruo, Sena, Tiegelong, Saijiao, and Gaerqin. In this study, we take the southeast of the Duolong mining area, referred to as the Gaerqin mining area, as the study area. Currently, the Duolong area is in the stage of pre-mining, so it is of great significance to know the extent of HM concentration distribution characteristics and assess soil HM contamination.
In this study, a total of 835 surface topsoil samples were collected, with a density of about four samples per km2. In the vicinity of the pre-set sampling point, radiating 10–30 m to the surrounding area, 4–5 points were combined equally to form a mixed sample, and attempts were made to ensure the consistency of the sampling depth and soil type. Subsequently, the concentrations of Cu, Pb, Zn, Cd, Cr, Hg, and As in the samples were measured by using X-ray fluorescence (XRF) and Atomic absorption spectrometry (AAS). More detailed information on sample preparation and analysis can be found in [49]. In this study, Pb, Zn, Cr, and Hg were selected to demonstrate the advantage of HMCA-Contour, and a statistic description of their data is shown in Table 4.
4.2. Pollution Indices Calculation
As mentioned above, original data are arranged in the form of Figure 2, and the file format is .xlsx. Through the data processing function window, the user loads the test data and chooses the target elements to be processed. After selecting theIgeo method, the user clicks the “calculate” button and theIgeo calculation message box then immediately pops up, requiring the user to input the corresponding background values (500, 500, 300, and 1.5 mg/kg for Pb, Zn, Cr, and Hg, respectively) for each selected heavy metal in the default form of “{500, 500, 300, 1.5}”. After inputting the background values and clicking the “OK” button, the calculatedIgeo result with its statistic description file (Table 4) is then saved in the output path directory.
4.3. Spatial Distribution of Heavy Metal Concentrations
Spatial distribution maps of Pb, Zn, Cr, and Hg were created from the original geochemical data by means of the “Inverse Distance to a power” method in this software. In the plotting function window, the user loads the calculated Igeo file and chooses Pb, Zn, Cr, and Hg as the target elements. Additionally, the user sets the parameters of “Axis Ticks”, “Axis Labels”, and “Axis Titles”. The next step is to click the “OK” button to plot, and the bottom progress bar shows specific information about the batch plot. After clicking the prefix and suffix buttons of both the map title and color scale title to activate text boxes, the user enters text into the corresponding text boxes.
The geochemical maps of soil heavy metals produced by HMCA-Contour are shown in Figure 6. The spatial distribution patterns of Pb and Zn concentration hotspots are almost the same in the center of the maps, indicating that they probably have similar geochemical behaviors. In comparison, the Cr concentration is mainly concentrated in the northeast of the study area. Ore and rock weathering may be responsible for the spatial distribution patterns of soil heavy metals. According to the geological facts in the study area, there are ultramafic rock units in the northeast part of the study area, which can be easily eroded and weathered into soil, resulting in the accumulation of Cr in the topsoil.
4.4. Heavy Metal Contamination Assessment Maps
Heavy metal contamination assessment maps were generated viaIgeo (Figure 7). When theIgeo value is higher than 0, soils are polluted. Therefore, there is no soil pollution caused by Pb, Zn, and Cd in this study area. However, based on the result map of CrIgeo , the most contaminated areas are located in the northeast part of the Duolong mining area, where theIgeo ranged from 0 to 5, indicating that the soil quality ranges from a moderately contaminated level to extremely contaminated level. Therefore, further studies are required to investigate and analyze the contamination source in this area. Knowing the extent of pollution and pollution sources is an essential and important step towards the prevention and control of heavy metal contamination.
5. Conclusions Contour maps are vital for allowing environmental geochemists to visualize the spatial distribution of soil heavy metals and contamination degree. A user-friendly, handy, and high-efficiency software, named HMCA-Contour, was designed based on Visual Basic 6.0 language and Surfer automation technology and a case study was introduced in this paper. HMCA-Contour can automatically calculate the soil pollution indices with their statistic description and plot the spatial distribution of soil heavy metal maps and soil contamination assessment maps. Meanwhile, several details of generated maps can be set according to the user’s requirements in the template setting function. Previous repetitive calculation and contour map drawing routines can be shortened by HMCA-Contour software. The real case study tells us that it only takes within a second to calculate the pollution indices and 3~5 s to draw a contour map for a single heavy metal. After a very small manual intervention, the user can easily achieve the expected results. The program can be widely used in pollution indices calculation and plotting maps of soil HMs’ spatial distribution and soil pollution levels. It also has implications for the development and improvement of heavy metal contamination assessment software.
Enlarge this image.Figure 3. The user interfaces of HMCA-Contour. (a) The main interface of HMCA-Contour software; (b) The 'Data Processing' interface; (C1) The 'EF' method calculation interface; (C2) The prompt window of EF_method calculation; (d) The interface of Plotting.
Enlarge this image.Figure 5. Plot setting interfaces of HMCA-Contour. (a) Color Scale; (b) Color Scale setting; (c) Axis Ticks; (d) Axis Labels; (e) Axis Titles; (f) Map Title & Scale Bar & North Arrow.
Enlarge this image.Figure 6. Spatial distribution maps of Pb, Zn, Cr, and Hg concentrations in the topsoil (N = 835).
Enlarge this image.Figure 7. Soil contamination assessment maps of Pb, Zn, Cr, and Hg byIgeo in the study area (N = 835).
| Class | Value | Soil Quality |
|---|---|---|
| 0 | Igeo≤ 0 | Practically uncontaminated |
| 1 | 0 <Igeo< 1 | Uncontaminated to moderately contaminated |
| 2 | 1 <Igeo< 2 | Moderately contaminated |
| 3 | 2 <Igeo< 3 | Moderately to heavily contaminated |
| 4 | 3 <Igeo< 4 | Heavily contaminated |
| 5 | 4 <Igeo< 5 | Heavily to extremely contaminated |
| 6 | 5 <Igeo | Extremely contaminated |
| Er | Single-Potential Ecological Risk (Er) | RI | Comprehensive Potential Ecological Risk (RI) |
|---|---|---|---|
| <40 | Low potential ecological risk | <90 | Low potential ecological risk |
| 40 ≤ Er < 80 | Moderate potential risk | 90 ≤ RI < 180 | Moderate potential ecological risk |
| 80 ≤ Er < 160 | Considerable potential risk | 180 ≤ RI < 360 | Strong potential ecological risk |
| 160 ≤ Er < 320 | High potential risk | 360 ≤ RI < 720 | Very strong potential |
| ≥320 | Significantly very high | ≥720 | Highly strong potential |
| Class | Pi | Appraisal Result | PN | Appraisal Result |
| 1 | Pi< 1 | Non-pollution | PN≤ 0.7 | Safety domain |
| 2 | 1 ≤Pi< 2 | Mild pollution | 0.7 <PN≤ 1 | Precaution domain |
| 3 | 2 ≤Pi< 5 | Moderate pollution | 1 <PN≤ 2 | Slightly polluted domain |
| 4 | Pi> 5 | Heavy pollution | 2 <PN≤ 3 | Moderately polluted domain |
| 5 | - | - | PN> 3 | Seriously polluted domain |
| Original Data | Calculated Data | ||||||||
|---|---|---|---|---|---|---|---|---|---|
| Variables | Min | Mean | Max | SD | Variables | Min | Mean | Max | SD |
| Pb | 3.15 | 30.64 | 361.86 | 26.34 | Pb Igeo | −4.84 | −7.9 | −1.05 | 0.72 |
| Zn | 38.25 | 98.05 | 1631.15 | 84.27 | Zn Igeo | −3.12 | −4.29 | 1.12 | 0.61 |
| Cr | 45.22 | 683.57 | 25751.6 | 1888.75 | Cr Igeo | −1.05 | −3.31 | 5.84 | 1.9 |
| Hg | 0.01 | 0.04 | 0.52 | 0.03 | Hg Igeo | −5.94 | −7.49 | −2.12 | 0.62 |
Sample quantity: 835; the unit of original data: mg/kg.
Author Contributions
Q.L. designed the study, developed the software, and drafted the manuscript; G.L. and W.C. revised and refined the manuscript; G.C. contributed to the data calculation and plotting. All authors have read and agreed to the published version of the manuscript.
Funding
This study was supported by the China Geological Survey Project (No. DD20190167), National Key Research and Development Project, Deep Resources Exploration and Mining (No. 2018YFC1803103), and National Natural Science Foundation of China (No. 21377098, 41807134).
Data Availability Statement
The data presented in this study are available on request from the corresponding author.
Acknowledgments
The authors wish to thank Zhifeng Zhang, Zebin Zhang, and Ruiliang Wang for testing this software. We also thank the anonymous reviewers and the associate editor for their comments, which improved the manuscript.
Conflicts of Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
1. Khan, K.; Lu, Y.; Khan, H.; Ishtiaq, M.; Khan, S.; Waqas, M.; Wei, L.; Wang, T. Heavy metals in agricultural soils and crops and their health risks in Swat District, northern Pakistan. Food Chem. Toxicol. 2013, 58, 449-458.
2. Fei, X.; Christakos, G.; Xiao, R.; Ren, Z.; Liu, Y.; Lv, X. Improved heavy metal mapping and pollution source apportionment in Shanghai City soils using auxiliary information. Sci. Total Environ. 2019, 661, 168-177.
3. Gobran, G.R.; Huang, P. Biogeochemistry of Trace Elements in the Rhizosphere, 1st ed.; Elsevier: Amsterdam, The Netherlands, 2005; pp. 1-182.
4. Noll, M.R. Trace elements in terrestrial environments: Biogeochemistry, bioavailability, and risks of metals. J. Environ. Qual. 2003, 32, 374.
5. Liu, G.; Wang, J.; Liu, X.; Liu, X.; Li, X.; Ren, Y.; Wang, J.; Dong, L. Partitioning and geochemical fractions of heavy metals from geogenic and anthropogenic sources in various soil particle size fractions. Geoderma 2018, 312, 104-113.
6. Liu, G.; Xue, W.; Wang, J.; Liu, X. Transport behavior of variable charge soil particle size fractions and their influence on cadmium transport in saturated porous media. Geoderma 2019, 337, 945-955.
7. Bai, J.; Cui, B.; Chen, B.; Zhang, K.; Deng, W.; Gao, H.; Xiao, R. Spatial distribution and ecological risk assessment of heavy metals in surface sediments from a typical plateau lake wetland, China. Ecol. Model. 2011, 222, 301-306.
8. Xie, Y.; Chen, T.B.; Lei, M.; Yang, J.; Guo, Q.J.; Song, B.; Zhou, X.Y. Spatial distribution of soil heavy metal pollution estimated by different interpolation methods: Accuracy and uncertainty analysis. Chemosphere 2011, 82, 468-476.
9. Tóth, G.; Hermann, T.; Szatmári, G.; Pásztor, L. Maps of heavy metals in the soils of the European Union and proposed priority areas for detailed assessment. Sci. Total Environ. 2016, 565, 1054-1062.
10. Liu, G.; Yu, Y.; Hou, J.; Xue, W.; Liu, X.; Liu, Y.; Wang, W.; Alsaedi, A.; Hayat, T.; Liu, Z. An ecological risk assessment of heavy metal pollution of the agricultural ecosystem near a lead-acid battery factory. Ecol. Indic. 2014, 47, 210-218.
11. Sutkowska, K.; Teper, L.; Czech, T.; Hulok, T.; Olszak, M.; Zogala, J. Quality of peri-urban soil developed from ore-bearing carbonates: Heavy metal levels and source apportionment assessed using pollution indices. Minerals 2020, 10, 1140.
12. Müller, G. Index of geoaccumulation in sediments of the Rhine River. Geojournal 1969, 2, 108-118.
13. Tomlinson, D.L.; Wilson, J.G.; Harris, C.R.; Jeffrey, D.W. Problems in the assessment of heavy-metal levels in estuaries and the formation of a pollution index. Helgoländer Meeresunters. 1980, 33, 566-575.
14. Covelli, S.; Fontolan, G. Application of a normalization procedure in determining regional geochemical baselines. Environ. Geol. 1997, 30, 34-45.
15. Hakanson, L. An ecological risk index for aquatic pollution control.a sedimentological approach. Water Res. 1980, 14, 975-1001.
16. Wang, X.; Dan, Z.; Cui, X.; Zhang, R.; Zhou, S.; Wenga, T.; Yan, B.; Chen, G.; Zhang, Q.; Zhong, L. Contamination, ecological and health risks of trace elements in soil of landfill and geothermal sites in Tibet. Sci. Total Environ. 2020, 715, 136639.
17. Zhen, J.; Pei, T.; Xie, S. Kriging methods with auxiliary nighttime lights data to detect potentially toxic metals concentrations in soil. Sci. Total Environ. 2019, 659, 363-371.
18. Petrik, A.; Thiombane, M.; Lima, A.; Albanese, S.; Buscher, J.T.; De Vivo, B. Soil contamination compositional index: A new approach to quantify contamination demonstrated by assessing compositional source patterns of potentially toxic elements in the Campania Region (Italy). Appl. Geochem. 2018, 96, 264-276.
19. Yan, W.; Mahmood, Q.; Peng, D.; Fu, W.; Chen, T.; Wang, Y.; Li, S.; Chen, J.; Liu, D. The spatial distribution pattern of heavy metals and risk assessment of moso bamboo forest soil around lead-zinc mine in Southeastern China. Soil Tillage Res. 2015, 153, 120-130.
20. Jiang, X.; Zou, B.; Feng, H.; Tang, J.; Tu, Y.; Zhao, X. Spatial distribution mapping of Hg contamination in subclass agricultural soils using GIS enhanced multiple linear regression. J. Geochem. Explor. 2019, 196, 1-7.
21. Lu, A.; Wang, J.; Qin, X.; Wang, K.; Han, P.; Zhang, S. Multivariate and geostatistical analyses of the spatial distribution and origin of heavy metals in the agricultural soils in Shunyi, Beijing, China. Sci. Total Environ. 2012, 425, 66-74.
22. Wang, X.S. Assessment of heavy metal pollution in Xuzhou urban topsoils by magnetic susceptibility measurements. J. Appl. Geophys. 2013, 92, 76-83.
23. Zhang, Y.; Li, S.; Chen, Z.; Wang, F.; Chen, J.; Wang, L. A systemic ecological risk assessment based on spatial distribution and source apportionment in the abandoned lead acid battery plant zone, China. J. Hazard. Mater. 2018, 354, 170-179.
24. Huang, Y.; Li, T.; Wu, C.; He, Z.; Japenga, J.; Deng, M.; Yang, X. An integrated approach to assess heavy metal source apportionment in peri-urban agricultural soils. J. Hazard. Mater. 2015, 299, 540-549.
25. Wang, M.; Bai, Y.; Chen, W.; Markert, B.; Peng, C.; Ouyang, Z. A GIS technology based potential eco-risk assessment of metals in urban soils in Beijing, China. Environ. Pollut. 2012, 161, 235-242.
26. Gong, Q.J.; Deng, J.; Xiang, Y.C.; Wang, Q.F.; Yang, L.Q. Calculating pollution indices by heavy metals in ecological geochemistry assessment and a case study in parks of Beijing. J. China Uni. Geosci. 2008, 19, 230-241.
27. Lin, J.; Yang, K.; Chen, H.; Wang, J.Y.; Zhang, C.G.; Huang, C.R. The visualization of soil moisture based on VB and Surfer. Adv. Mater. Res. 2013, 774-776, 1859-1864.
28. Hou, C.Q.; Li, Z.Q.; Wu, X.B. Application of VB and Surfer Automation technique in geochemical exploration. Comput. Tech. Geophys.Geochem. Explor. 2010, 32, 565-570.
29. Zhang, T.J.; Wang, X.W.; Zhang, H.; Liu, X.W. Application of surfer automation technique in automatic mapping for meteorological data. Arid Meteorol. 2007, 2, 90-93.
30. Ni, J.J.; Zhang, X.H. Fast visualizing of hydrological data with VB 6.0 embedded with compiler of Surfer 7.0. Hydrogr. Surv. Charting 2005, 1, 64-66.
31. Loska, K.; Wiechuła, D.; Korus, I. Metal contamination of farming soils affected by industry. Environ. Int. 2004, 30, 159-165.
32. Zoller, W.H.; Gladney, E.; Duce, R.A. Atmospheric concentrations and sources of trace metals at the South Pole. Science 1974, 183, 198-200.
33. N'Guessan, Y.M.; Probst, J.L.; Bur, T.; Probst, A. Trace elements in stream bed sediments from agricultural catchments (Gascogne region, S-W France): Where do they come from? Sci. Total Environ. 2009, 407, 2939-2952.
34. Woitke, P.; Wellmitz, J.; Helm, D.; Kube, P.; Lepom, P.; Litheraty, P. Analysis and assessment of heavy metal pollution in suspended solids and sediments of the river Danube. Chemosphere 2003, 51, 633-642.
35. Sutherland, R.A. Bed sediment-associated trace metals in an urban stream, Oahu, Hawaii. Environ. Geol. 2000, 39, 611-627.
36. Lee, D.S.; Garland, J.A.; Fox, A.A. Atmospheric concentrations of trace elements in urban areas of the United Kingdom. Atmos. Environ. 1994, 28, 2691-2713.
37. Loska, K.; Cebula, J.; Pelczar, J.; Wiechuła, D.; Kwapuliński, J. Use of Enrichment, and Contamination Factors Together with Geoaccumulation Indexes to Evaluate the Content of Cd, Cu, and Ni in the Rybnik Water Reservoir in Poland. Water Air Soil Pollut. 1997, 93, 347-365.
38. Schiff, K.C.; Weisberg, S.B. Iron as a reference element for determining trace metal enrichment in Southern California coastal shelf sediments. Mar. Environ. Res. 1999, 48, 161-176.
39. Tam, N.F.Y.; Yao, M.W.Y. Normalisation and heavy metal contamination in mangrove sediments. Sci. Total Environ. 1998, 216, 33-39.
40. Cui, J.; Zang, S.; Zhai, D.; Wu, B. Potential ecological risk of heavy metals and metalloid in the sediments of Wuyuer River basin, Heilongjiang Province, China. Ecotoxicology 2014, 23, 589-600.
41. Jiang, H.H.; Cai, L.M.; Wen, H.H.; Hu, G.C.; Chen, L.G.; Luo, J. An integrated approach to quantifying ecological and human health risks from different sources of soil heavy metals. Sci. Total Environ. 2020, 701, 134466.
42. Li, L.; Wu, J.; Lu, J.; Min, X.; Xu, J.; Yang, L. Distribution, pollution, bioaccumulation, and ecological risks of trace elements in soils of the northeastern Qinghai-Tibet Plateau. Ecotox. Environ. Safe. 2018, 166, 345-353.
43. Lee, C.S.L.; Li, X.D.; Shi, W.Z.; Cheung, S.C.N.; Thornton, I. Metal contamination in urban, suburban, and country park soils of Hong Kong: A study based on GIS and multivariate statistics. Sci. Total Environ. 2006, 356, 45-61.
44. Nemerow, N.L. Scientific Stream Pollution Analysis, 1st ed.; Scripta Book Company: Washington, DC, USA, 1974; pp. 1-68.
45. 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.
46. Mazurek, R.; Kowalska, J.; Gąsiorek, M.; Zadrożny, P.; Józefowska, A.; Zaleski, T.; Kępka, W.; Tymczuk, M.; Orłowska, K. Assessment of heavy metals contamination in surface layers of Roztocze National Park forest soils (SE Poland) by indices of pollution. Chemosphere 2017, 168, 839-850.
47. Zhou, J.; Feng, K.; Pei, Z.; Lu, M. Pollution assessment and spatial variation of soil heavy metals in Lixia River Region of Eastern China. J. Soils Sediments 2016, 16, 748-755.
48. Wang, Z.; Liu, X.; Qin, H. Bioconcentration and translocation of heavy metals in the soil-plants system in Machangqing copper mine, Yunnan Province, China. J. Geochem. Explor. 2019, 200, 159-166.
49. EPAC. The Technical Specification for Soil Environmental Monitoring; China Environmental Protection Agency: Beijing, China, 2005.
50. Chen, H.H.; Li, X.; Ding, W.X. Twelve kinds of gridding methods of Surfer 8.0 in isoline drawing. Chin. J. Eng. Geophys. 2007, 4, 52-57.
51. Liu, Q.P.; Deng, S.Q.; Zhao, Y.Y.; Li, X.S. Sources Analysis and Element Geochemical Characteristics of Regional Soil Heavy Metals from the Ga'erqin Ore Deposit, Tibet. Acta Geosci. Sin. 2018, 39, 481-490.
Qingping Liu
1,2,
Guannan Liu
2,*,
Wei Chen
2,* and
Guoliang Chen
3
1School of the Earth Science and Resources, China University of Geosciences, Beijing 100081, China
2MNR Key Laboratory of Metallogeny and Mineral Assessment, Institute of Mineral Resources, Chinese Academy of Geological Sciences, Beijing 100037, China
3Tibet Huatailong Ming Co., Ltd., Lhasa 850212, China
*Authors to whom correspondence should be addressed.
© 2021. This work is licensed under http://creativecommons.org/licenses/by/3.0/ (the “License”). Notwithstanding the ProQuest Terms and Conditions, you may use this content in accordance with the terms of the License.