Sources and Risk Characteristics of Heavy Metals in Plateau Soils Predicted by Geo-Detectors

Abstract

The ecological environment of the remote plateaus has been a challenge plaguing many countries due to global warming, increased human activities, and frequent occurrence of various heavy metal (HM) pollutants. This paper analyzes the primary factors and potential susceptible regions’ characteristics related to soil HMs (As, Cd, Pb, Cr, Ni, and Zn) in the southern Tibet Plateau (TP) using Geo-detectors and a set of multi-source data from remote sensing and geographic and geological information. The geoaccumulation index showed that the As pollution was the most serious. The precipitation, pH, geological age (GA), and rock type (RT) were the most critical factors affecting HMs. Additionally, soil HMs were relatively unaffected by land use and clay. Based on the risk detector, the most key impact types (range) of the predominant factors of HM’s sources were identified. For example, precipitation (100–200 mm), pH (7–7.5), DEM (>5500 m), RT (ophiolite), and GA (Ordovician) had the highest average concentration of As in each type (ranges) of the predominant factors. This work provides new data on the extent of soil HM contamination in the southern TP and predicts vulnerable areas of HM contamination, providing an important scientific basis for monitoring and managing HM in remote areas at high altitudes.

1. Introduction

The soil is an indispensable resource for human survival and agricultural production [1,2]. However, heavy metal (HM) pollution of the soil is a worldwide environmental concern, as it threatens the future of global food production and human health [3,4,5]. For example, arsenic (As) is harmful to the skin and affects blood pressure and the cardiovascular system, and cadmium (Cd) and lead (Pb) affect the nervous system, kidneys, and bones [6,7,8]. In recent decades, with urbanization and industrialization, HM pollution in soil has been developing continuously. Soil investigation is challenging and costly due to the complexity of soil properties and the concealment of HM pollution [8]. Thus, the identification and risk assessment of soil contamination sources are of great practical importance to ensure food security and the rational use of land resources.

Many assessments and source identification methods have been applied in recent years to decipher soil HMs [9,10]. The potential ecological risk index (PERI) and the geoaccumulation index (I_geo) are the most prevalent methods for evaluating HM pollution and its risk characteristics [11]. The main forms of traditional source apportionment include cluster analysis (CA), positive matrix factorization (PMF) and principal component analysis (PCA) [12,13,14,15]. Although it is possible to infer factors that might influence the HM distribution, the interpretation of the evaluation theory is based on the subjective judgment of the evaluator. In addition, these methods cannot quantitatively analyze the degree of influence of each factor and the interaction among the driving factors. Thus, it is difficult to ensure the accuracy of risk management [16]. A Geo-detector is a statistical method for spatial data which can quantitatively reveal the relationship between geographical phenomena (e.g., the spatial distribution of HMs) and various geographical factors (e.g., pH or soil type) [17]. Compared to PCA and PMF, Geo-detectors provide more robust statistics for revealing causal relationships between independent and dependent variables [18]. Recently, an increasing number of studies have applied geographical detectors to quantify the intrinsic associations of spatial heterogeneity among different influencing factors and variables, such as HM pollution [3,4], vegetation coverage [19], and urban green space [20]. Zhao, et al. [4] explored the relationship between high background levels of Cd in soil and soil type, landform type, rock type, and geological age using Geo-detectors and predicted the risk zone of Cd in soil. According to Qiao, et al. [3], the main influencing factors of soil pollution were identified by Geo-detectors as mining activities and the pH value, providing a theoretical basis for soil pollution control and restoration. Therefore, this study uses Geo-detectors to explore the relationship between HM distribution and influencing factors to determine the potential susceptible area characteristics affecting soil HMs.

Previous works have revealed that human activities (transportation, industry, and mining activities), topography factors, climate factors, and geological processes (e.g., volcanic activity and weathering of parent rocks) contribute greatly to soil HM accumulation [4,21,22,23]. The above results are of great significance for understanding the spatial heterogeneity of soil HMs. However, few studies have attempted to comprehensively analyze the link between HM distribution and geographic factors, regional geological conditions, and anthropogenic activities, and the extent of the dominant factors.

The Tibetan Plateau (TP) is the most unique geological–geographical–ecological unit on earth, characterized by a sensitive climate, complex geology, and fragile ecology [24]. It used to be considered less affected by industrial and agricultural activities and less contaminated with HMs. However, due to the increasing global warming, exploitation of natural resources, and urbanization, unexpected HM contamination has occurred in remote areas of land previously considered clean [25]. Yu, et al. [26] suggested that climate warming has exacerbated the melting of permafrost on the TP, accelerating the introduction of trace elements (e.g., As) into the environment. In addition, anthropogenic industrial sources have exacerbated the accumulation of HMs in soils [27]. Notably, the TP is an area with a high prevalence of endemic diseases, and the enrichment and distribution of HMs in the soil are of concern [28]. Currently, out of the relevant studies on soil HM pollution on the TP, most of them focus on a specific region or a single influencing factor (e.g., major roads or major mining areas) [29,30,31], lacking an overall analysis of large-scale watersheds [32]. Additionally, the data of soil HM sampling was recorded some time ago, and the distribution of HM in soil may have changed under the influence of human activities and climate change. Therefore, a comprehensive analysis of the driving factors and multi-factor interaction of the spatial distribution of HMs in the soil is critical to guide the response to the subsequent ecological environment issues. The Yarlung Tsangpo River basin (YTRB), located in the southern part of the TP of China, is a tectonic suture zone formed by the Himalayan orogeny. It is an area of concentrated population and agricultural development, and is known as the most essential “food store” of the TP [33]. Moreover, the YTRB is an ideal area to study the variations of HMs because of its longitudinal geography (>2000 km in length from east to west), wide range of elevation coverage (132–7168 m), climatic differences, and complex geologic age and rock types [34]. To the best of our knowledge, studies on the prediction of HMs in the highland environment have not been conducted. Therefore, in such a climate-sensitive and ecologically fragile region, it is crucial to explore the driving role of geological conditions, human activities, and climatic influences on the distribution of HMs in plateau soils for guidance [35].

Given the lack of widespread knowledge on the distribution of soil HMs in the YTRB and considering its hidden health effects, the main objectives of this study were: (1) to present the distribution characteristics of HMs in soils of the YTRB; (2) to identify the main sources of HMs by using Geo-detectors and environmental factor data (including remote sensing data, digital elevation model, geological information, and other geospatial data); and (3) to reveal the regional characteristics of potential susceptibility to pollution and provide a scientific basis for soil management on the TP and in similar high altitude areas.

2. Materials and Methods

2.1. Study Area

The collision of the Indian and Eurasian plates resulted in the formation of the TP. The YTRB is a densely populated and agriculturally developed area in the TP and an important highland barley-producing area in China. Since the Quaternary, the YTRB has experienced strong extrusion and formed complex fracture structures [36]. The Yarlung Tsangpo River (YTR) is one of the highest large rivers in the world, with an average elevation above 3000 m and a total length of 2840 km [33]. The climate of the YTRB is a highland temperate semi-arid monsoon climate zone. Due to its abundant water volume and large and concentrated drop, the YTRB is rich in hydraulic resources, ranking second in China after the Yangtze River.

2.2. Data Source and Method

The framework of prediction methods for soil HM sources and risk characteristics is shown in Figure 1. First, the content of heavy metals in soil was detected. The main influencing factors, interactions, and potential risk areas of heavy metals in soil were analyzed using multi-source data and Geo-detectors. In addition, a soil HM pollution assessment method was used to analyze the soil pollution level.

2.2.1. Sample Collection and Processing

In July 2020, 153 soil samples were collected from the YTRB (Figure 2). A five-point mixed sampling (200 m × 200 m grid) method was used to collect approximately 1000 g of surface soil (0–20 cm), which was thoroughly mixed and put into sample bags. The geographical coordinates and elevation of the sample points were recorded. The collected soil samples were returned to the laboratory for air-drying, and the debris (plant roots, sticks, biological residues) were removed. The soil samples were ground, passed through a 0.15 mm nylon sieve and bagged for further testing. An X-ray fluorescence (XRF) spectrometer (Niton XL 3t, Niton, Greenwich, CT, USA) was used for the determination of the HM content in the soil samples. The detection limits of As, Cd, Pb, chromium (Cr), nickel (Ni), and zinc (Zn) were 0.7, 0.05, 0.7, 5, 1, and 0.5 mg/kg, respectively. The data precision is expressed as relative standard deviation (RSD), and the reliability as the relative error (RE). The RE of As, Cd, Pb, Cr, Ni, and Zn in standard soil (GBW07403), as determined by the XRF method, was 5.79, −5.71, 9.74, 2.25, 9.94, and 15.4%, respectively. The RSD of As, Cd, Pb, Cr, Ni, and Zn, as determined by the XRF method, was 6.96, 7.12, 0.52, 8.41, 9.77, and 9.55%, respectively. Both XRF and inductively coupled plasma mass spectrometry were used to determine As, Cd, Pb, Cr, Ni, and Zn levels in the soil, and the test results indicated that the XRF method could achieve a higher accuracy in the determination of HM levels in the soil samples [37,38].

2.2.2. Remote Sensing Data

The Chinese land use remote sensing monitoring dataset is based on the Landsat remote sensing images of the United States as the main information source. China’s national-scale land use/land cover thematic database (2020) was constructed through manual visual interpretation. The data are classified into the 6 categories of arable land, forest land, grassland, water area, construction land, and unused land using the first level.

2.2.3. Geological Information

Geological factors include rock type (RT) and geological age (GA) [4]. The soil is weathered by the parent rock, and the crustal trace elements are not uniformly distributed. Additionally, the weathering of the rocks may have varied in different GA. Thus, these two factors are mainly considered in terms of geological factors. The map information of RT and GA is integrated from the Geological Atlas of China (2001).

2.2.4. Other Geospatial Data

Topographic and hydrological characteristics specifically include a DEM and precipitation, which influence runoff and erosion processes and ultimately the transport of HMs. Mining activities are considered critical anthropogenic factors influencing the distribution of HMs [4]. pH and soil texture (clay content, sand content, silt content) affect the migration of HMs [3]. Mining activities, soil type, soil content, land use (LUCC), and pH were obtained by the Resource and Environmental Science Data Center (RESDC) of the Chinese Academy of Sciences (http://www.resdc.cn). All influencing factors were plotted and classified using original classification or classification algorithms (Figure 3), such as Natural Break (Jenks), K-means, and equal intervals, ensuring large local variance between sub-regions and slight local variance within sub-regions. Secondly, the number of stratifications should be reasonable to improve the computational efficiency of the geographical detector.

2.3. Geo-Detector Method

The Geo-detector is a spatial analysis model that detects the intrinsic linkages according to the spatial heterogeneity among variables. It contains four components, the factor detector, interaction detector, risk detector, and ecological detector [14]. The detailed mechanism of the Geo-detector refers to the literature [39]. The primary hypothesis of the geographic detector approach is such that if a factor (D) is related to the spatial heterogeneity of HMs (H), then the spatial distribution of H is analogous to D. The calculation formulas are as follows:

$$q=1-\frac{1}{n{\sigma }_H^2}{\displaystyle \sum }_{i=1}^m\left(n_{D,i}·{\sigma }_{H_{D,i}}^2\right)$$

where D means the influencing factor (e.g., pH); m is the number of classifications of factor D (stratification); H represents the spatial variance of soil HM concentrations; q refers to the force of the ranks D on H; n denotes the total number of study samples; and σ² is the variance of the influencing factor (D). Generally, a high q implies a greater influence of the influencing factors on the spatial variation of HM concentrations. The interaction detector determines the effect of a combination of two different factors on HMs that tend to enhance, weaken, or become independent of each other. The risk detector compares the differences in mean values of HMs in soil for different influencing factors and looks for potential contamination zones for specific characteristics.

2.4. Moran Index

Moran indices have been widely used in analyzing spatial correlations and can be divided into global Moran indices and local Moran indices [40,41]. The global Moran index is applied to identify whether the attribute values are spatially clustered. In contrast, the local Moran index can serve as a further indication of clustered regions with high or low attribute values. The calculation formulas are as follows:

$$\begin{gathered} I=\frac{n\times {{\displaystyle \sum }}_{i=1}^n{{\displaystyle \sum }}_{j=1}^n{\omega }_{i,j}\left(x_i-\overline{x}\right)\left(x_j-\overline{x}\right)}{{{\displaystyle \sum }}_{i=1}^n{{\displaystyle \sum }}_{j=1}^m{\omega }_{i,j}\times {{\displaystyle \sum }}_{i=1}^n{\left(x_i-x_j\right)}^2} \\ {I}_i=\frac{\left(x_i-\overline{x}\right){{\displaystyle \sum }}_{j=1}^n{\omega }_{i,j}\left(x_i-\overline{x}\right)}{{{\displaystyle \sum }}_{i=1}^n{\left(x_i-\overline{x}\right)}^2} \end{gathered}$$

where I represent the global Moran index and I_i represents the local Moran index; a detailed description can be referred to Wu, et al. [37].

2.5. Methods for Assessment of HM Pollution

The contamination factor (CF), geo-accumulation index (I_geo), and potential ecological risk index (PERI) were used to evaluate HM pollution. These three indices have been adopted for risk evaluation of soil HM contamination [30].

2.5.1. Contamination Factor (CF)

CF = C_Sample / C_Standard

where C_Sample is the concentration (mg/kg) of HM in the sample and C_Standard is the Environmental Class 2 Standard (mg/kg). The detailed CF classification is shown in

Table 1. Grading criteria for contamination factor (CF), geoaccumulation (I_geo), and potential ecological risk index (PERI).

CF	Contamination Factor	Igeo	Pollution Class	PERI	Risk Class
CF < 1	low contamination	Igeo ≤ 0	practically unpolluted	PERI < 110	low
1 ≤ CF < 3	moderate contamination	0< Igeo ≤1	unpolluted to moderately polluted	110 ≤ PERI < 220	moderate
3 ≤ CF ≤ 6	considerable contamination	1< Igeo ≤2	moderately polluted	220 ≤ PERI < 440	strong
CF > 6	very high contamination	2 < Igeo ≤ 3	moderately to strongly polluted	440 ≤ PERI < 880	very strong
		3 < Igeo ≤ 4	strongly polluted	PERI > 880	highly-strong
		4 < Igeo≤ 5	strongly to extremely polluted
		Igeo > 5	extremely polluted

.

2.5.2. Geoaccumulation Index (I_geo)

$$I_{geo}=lo{g}_2\frac{C_n}{1.5\times B_n}$$

The I_geo shows the accumulation of HMs in soil [42]. C_n and B_n are the concentrations of HMs and geochemical background values, respectively. The detailed I_geo classification is shown in Table 1.

2.5.3. Potential Ecological Risk Index (PERI)

$$PERI=T_r^i\frac{C^i}{C_n^i}$$

The PERI was established by Hakanson [43] to assess the potential ecological risk of HMs in soil. $T_r^i$ is the toxicity factor for HMs. Cⁱ and $C_n^i$ are the concentration of HMs in the sample and the reference value for HMs. The toxicity response factors for Pb, As, Cd, Cr, Zn, and Ni were 1.5, 10, 30, 2, 1, and 5 [43]. The detailed PERI classification is shown in Table 1.

2.6. Statistical Analysis

The data were organized, classified, and calculated using Microsoft Excel software. The spatial distribution and spatial autocorrelation of HMs in the soil were analyzed using the ArcGIS 10.2 software. The inverse distance weighting method was used to analyze the spatial distribution, and the global Moran’s I index was used to estimate the spatial autocorrelation. Plots were generated using Origin Pro 2022 and Microsoft PowerPoint software.

3. Results and Discussion

3.1. Statistics of the HMs Present in the Agricultural Soil of the YTRB

The statistical results of HM contents in YTRB soils are shown in

Table 2. Descriptive statistics analysis results for soil heavy metals.

	Range	Mean	Median	Standard Deviation	Variable Coefficient	Soil Environment Quality a	Soil Background Value/(mg·kg−1)
	mg·kg−1						TP b	China c	World c
As	1.7–81.3	17.9	15.5	11.2	62.6%	25	18.7	11.2	6
Pb	10.7–225.6	32.4	29.0	20.8	64.3%	170	28.9	26	35
Cd	0.04–0.47	0.14	0.12	0.09	62.3%	0.6	0.08	0.097	0.35
Cr	0.8–270.2	63.1	57.4	40.7	64.5%	250	77.4	61	70
Ni	4.4–446.9	43.7	32.1	53.3	121.7%	190	32.1	26.9	50
Zn	25.0–300.9	100.9	52.4	52.4	51.9%	300	73.7	74.2	90

^a Risk-screening values for soil contamination of agricultural land come from the Soil Environmental Quality Risk Control Standard for Soil Contamination of Agricultural Land (GB 15618—2018). ^b The background values of the Tibet Plateau are from Zhang et al., 1994 (in Chinese). ^c Data are derived from Background values of soil elements in China, 1990 (in Chinese).

. The average values of the six HMs in the YTRB were lower than the national soil environmental quality standards. The concentrations of Pb, Zn, Ni, and Cd, exceeded the soil environmental background values of the TP. Compared with the background values of Chinese soils, the contents of HMs in the YTRB were all higher, especially for As and Ni. In addition, only As (17.9 mg/kg) surpassed the world soil background value (6.0 mg/kg). The coefficients of variation of topsoil HMs in the YTRB were Ni > Cr > Pb > As > Cd > Zn in order of magnitude. Among them, the variability of Ni was extreme (>100%) and Cr, Pb, As, Cd, and Zn showed high variability (>50%), probably influenced by complex geology, geomorphology, climate, and human activities [29]. Therefore, soil HMs have obvious spatial heterogeneity and point source pollution. According to the results of the I_geo (Figure 4c), the soil HMs in the YTRB ranged from unpolluted to extremely polluted (−5.97–4.76), with an average I_geo value of 2.39 for As, followed by Ni (0.32), Zn (0.27), Pb (0.16), Cd (−0.12), and Cr (−0.13). Arsenic pollution levels (−0.85–4.76) were the highest, with moderate to strong pollution accounting for 55.8% of soil samples.

3.2. Spatial Distribution of the Soil HMs in the YTRB

The spatial distribution characteristics of the HMs (As, Ni, Cd, Pb, Cr, and Zn) present in the YTRB varied greatly, as shown in Figure 4a. High levels of soil As are distributed in the middle and upper reaches of the study area, which is similar to the spatial distribution of As in the surface water of the YTR [44]. It may be mainly influenced by natural factors, such as the perennial outflow of high-As geothermal water [45]. The high values of Cr and Ni are mainly distributed in the middle reaches and the upper reaches, respectively, which may be related to the natural weathering of rocks with high background content. The spatial distribution of Cd, Pb, and Zn are analogous, with high values mainly concentrated in the midstream and downstream areas, which may be related to human activities, such as mining [27]. Pearson correlation analysis showed that Cd, Pb, and Zn showed significant correlations (Figure 4b). Additionally, Cr and Ni also showed high positive correlations. However, the correlation between As and other HMs is lower.

The results of the global Moran’s I index analysis of the soil HMs in the YTRB are shown in

Table 3. Moran’s I of soil heavy metal content.

Elements	Moran’s I	Z	p
Pb	0.285	3.742	0.001
Zn	0.465	5.92	0.001
Ni	0.283	3.67	0.001
Cr	0.325	4.21	0.001
As	0.106	1.45	0.147
Cd	−0.045	−0.446	0.655

. The Monte Carlo iterative method was selected to verify the significance of the obtained Moran’s I index values. Conversely from those for Cd and As, the Z values of Pb, Cr, Ni, and Zn were greater than 1.96 (p < 0.05), indicating that the spatial autocorrelation for the Moran’s I index of Pb, Cr, Ni, and Zn was statistically significant (Table 3). The Moran’s I index values for these four elements were greater than 0 (0.2853 for Pb, 0.3255 for Cr, 0.2827 for Ni, and 0.4652 for Zn). From a global perspective, Pb, Cr, Ni, and Zn have positive spatial correlation distribution in the YTRB, indicating that these four HMs have a clustered distribution in the study area. In contrast, As and Cd have a discrete distribution. The results for the local Moran’s I index of Pb, Cr, Ni, and Zn in the YTRB are shown in Figure 5. High–high areas were mainly distributed in the central region for Ni, in the middle and lower reaches for Pb and Zn, and the middle and upper reaches for Cr. Therefore, the spatial distribution of the HM content in the soil predicted by the interpolation method was consistent with the spatial correlation distribution of the local Moran’s I index.

3.3. Importance Evaluation of Influencing Factors

The explanatory powers of 10 environmental variables on soil HM concentration in the YTRB, generated by the factor detector model, were ranked, as shown in Figure 6. All factors were statistically significant at the 0.05 level, indicating that each driving factor had a significant contribution to the spatial distribution of the HMs. The top three cumulative explanatory powers of As, Ni, Cd, Zn, Cr, and Zn were 49.8%, 48.4%, 73.2%, 48.4%, 77.8%, and 96.2%, respectively. In terms of the important sequence of the explanatory power of the variables, precipitation, pH, GA, and RT are the most important influencing factors. It is noteworthy that the DEM has strong explanatory power for As, while the effect on other HMs is relatively weak. In addition, soil HMs were relatively unaffected by LUCC and clay.

The influence of precipitation may cause the diffusion of atmospheric HMs to the surface. In the case of the Pearl River Delta, higher As levels under rainfall conditions were due to more elevated atmospheric As concentrations [46]. In addition, precipitation affects rocks and organic matter, increasing the migration of eroded material to the ambient water system, thereby diluting HMs in the topsoil [23,47]. Precipitation conditions vary greatly from upstream (<100 mm) to downstream (>800 mm) within the YTRB (Figure 3), and differences in weathering of rocks and melting of permafrost may also contribute to heterogeneity in the distribution of HMs in soils [48,49]. In general, pH value is the critical factor affecting the spatial distribution of HM concentration in soil. The soil pH contributes to the dissolution and absorption of HMs, thus affecting their ability to migrate [50]. Additionally, our results showed that GA and RT also have a considerable effect on the HMs in the YTRB. Studies have demonstrated that the weathering of the parent rock is the primary source of HMs in the soil [4]. The YTRB has undergone tectonic evolution, hydrothermal activity, and diagenesis during its long geological history, resulting in different GAs and RTs, leading to differences in trace elements [51]. For example, it was reported that the concentration of As in Permian siltstone is as little as 5.8 mg/kg, and in Pliocene gabbro, Miocene granite, and Cretaceous limestone, it is 10.5–22.5 mg/kg; however, the concentration of As in Eocene gabbro and andesite is even as high as 142 mg/kg [45]. In addition, the DEM also had a high explanatory power (12%) for the spatial differentiation of soil As. This result is aligned with the findings by Shi, et al. [18], who pointed out that the spatial distribution of As in urban soils is primarily governed by geological and topographical factors. Other influencing factors may lead to a more complicated spatial distribution of the HM concentrations, resulting in greater spatial variability.

Table 4. Effects of the interaction between each pair of driving factors on the enrichment of the As concentration in the YTRB.

(X1 ∩ X2)	DEM	Precipitation	LUCC	Clay	Sand	Silt	pH	GA	RT	Mining
DEM	0.12
Precipitation	0.27	0.22
LUCC	0.18	0.25	0.05
Clay	0.13	0.25	0.09	0.05
Sand	0.17	0.25	0.09	0.10	0.04
Silt	0.16	0.26	0.10	0.09	0.09	0.03
pH	0.26	0.24	0.20	0.21	0.20	0.21	0.16
GA	0.24	0.30	0.14	0.15	0.16	0.16	0.27	0.08
RT	0.22	0.28	0.16	0.15	0.15	0.15	0.24	0.23	0.08
Mining	0.19	0.27	0.10	0.10	0.11	0.12	0.22	0.18	0.19	0.04

,

Table 5. Effects of the interaction between each pair of driving factors on the enrichment of the Ni concentration in the YTRB.

(X1 ∩ X2)	DEM	Precipitation	LUCC	Clay	Sand	Silt	pH	GA	RT	Mining
DEM	0.01
Precipitation	0.35	0.24
LUCC	0.05	0.27	0.02
Clay	0.03	0.27	0.03	0.00
Sand	0.05	0.30	0.04	0.02	0.01
Silt	0.06	0.30	0.05	0.02	0.03	0.01
pH	0.21	0.30	0.17	0.16	0.17	0.19	0.14
GA	0.14	0.36	0.13	0.11	0.13	0.14	0.28	0.10
RT	0.14	0.32	0.15	0.12	0.15	0.16	0.28	0.22	0.10
Mining	0.11	0.32	0.11	0.06	0.09	0.10	0.22	0.23	0.24	0.05

,

Table 6. Effects of the interaction between each pair of driving factors on the enrichment of the Cd concentration in the YTRB.

(X1 ∩ X2)	DEM	Precipitation	LUCC	Clay	Sand	Silt	pH	GA	RT	Mining
DEM	0.03
Precipitation	0.43	0.37
LUCC	0.06	0.40	0.02
Clay	0.06	0.38	0.08	0.02
Sand	0.15	0.40	0.11	0.09	0.07
Silt	0.12	0.39	0.09	0.08	0.09	0.06
pH	0.32	0.41	0.27	0.27	0.30	0.31	0.24
GA	0.19	0.43	0.18	0.16	0.22	0.21	0.35	0.13
RT	0.17	0.43	0.14	0.15	0.21	0.19	0.33	0.27	0.10
Mining	0.09	0.44	0.06	0.06	0.11	0.11	0.32	0.19	0.17	0.02

,

Table 7. Effects of the interaction between each pair of driving factors on the enrichment of the Pb concentration in the YTRB.

(X1 ∩ X2)	DEM	Precipitation	LUCC	Clay	Sand	Silt	pH	GA	RT	Mining
DEM	0.02
Precipitation	0.34	0.29
LUCC	0.07	0.35	0.01
Clay	0.05	0.31	0.05	0.01
Sand	0.09	0.32	0.07	0.05	0.05
Silt	0.08	0.32	0.07	0.03	0.07	0.03
pH	0.24	0.32	0.18	0.16	0.17	0.16	0.11
GA	0.12	0.45	0.13	0.11	0.15	0.14	0.23	0.08
RT	0.10	0.43	0.10	0.08	0.14	0.12	0.20	0.15	0.06
Mining	0.12	0.34	0.08	0.04	0.09	0.07	0.16	0.17	0.17	0.01

,

Table 8. Effects of the interaction between each pair of driving factors on the enrichment of the Cr concentration in the YTRB.

(X1 ∩ X2)	DEM	Precipitation	LUCC	Clay	Sand	Silt	pH	GA	RT	Mining
DEM	0.03
Precipitation	0.49	0.43
LUCC	0.09	0.45	0.03
Clay	0.05	0.44	0.04	0.00
Sand	0.15	0.45	0.11	0.08	0.07
Silt	0.12	0.45	0.10	0.08	0.10	0.06
pH	0.31	0.46	0.26	0.25	0.27	0.29	0.23
GA	0.17	0.49	0.14	0.11	0.21	0.21	0.33	0.08
RT	0.21	0.51	0.18	0.16	0.22	0.21	0.35	0.28	0.12
Mining	0.17	0.54	0.17	0.13	0.24	0.19	0.36	0.27	0.28	0.10

and

Table 9. Effects of the interaction between each pair of driving factors on the enrichment of the Zn concentration in the YTRB.

(X1 ∩ X2)	DEM	Precipitation	LUCC	Clay	Sand	Silt	pH	GA	RT	Mining
DEM	0.10
Precipitation	0.54	0.50
LUCC	0.14	0.52	0.07
Clay	0.12	0.51	0.09	0.03
Sand	0.26	0.56	0.19	0.16	0.11
Silt	0.21	0.55	0.17	0.12	0.16	0.09
pH	0.40	0.55	0.35	0.35	0.39	0.39	0.33
GA	0.27	0.57	0.24	0.21	0.29	0.27	0.42	0.13
RT	0.24	0.57	0.19	0.16	0.24	0.22	0.40	0.27	0.10
Mining	0.18	0.59	0.15	0.13	0.24	0.22	0.45	0.28	0.25	0.08

represent “effects of interaction between each pair of driving fac-tors”—combined effects are listed outside the main diagonal. On the main diagonal of Table 4, Table 5, Table 6, Table 7, Table 8 and Table 9, for the same pairs of driving factors (DEM and DEM), their individual factors from Figure 6 are listed. The interaction detector showed that the combined effect of two drivers is greater than that of the individual factors. For example, the q values of precipitation, pH, and DEM for As in the YTRB were 21.5%, 16.3%, and 12.0% respectively, but the q values for their interactions were 24.0% for precipitation ∩ pH and 27.0% for precipitation ∩ DEM (Table 4). This result indicates that the distribution characteristics of the HMs in the soil are not simply superimposed by driving factors, but the interaction between different factors influences their spatial variation.

3.4. Identification of Risk Characteristics of HMs

The results of the ecological risk evaluation model are shown in Figure 7b. The ecological risk of Cd was the highest, followed by that of As, and the lowest was for Zn. Overall, the average ecological risk of HMs was low. The results of the pollution index are shown in Figure 7c. The ranking order of the six HMs was As (20.3%) > Cr, Zn, and Ni (2%) > Pb (0.7%) > Cd (0). HMs with a higher proportion of pollution indices did not indeed lead to higher ecological risks; mainly, the differences were in toxicity coefficients, resulting in different risks to the ecosystem. Generally, the soil HMs in the study area are still in a “clean” state, but some points have exceeded the reference values.

To identify potential vulnerable areas for HMs pollution, this study selected the type or range of the top five dominant factors with the most significant influence on HM accumulation based on risk detectors, as shown in Figure 7a. For example, for As, precipitation (100–200 mm), pH (7–7.5), the DEM (>5500 m), RT (ophiolite), and GA (Ordovician) had the highest average concentration of As in each type (range) of the dominant factors. Considering the results, although the dominant factors of some HMs are the same, there are differences in the types or ranges. For example, the most dominant driver of HMs is precipitation, and the dominant range for high values of As is 100–200 mm, while the dominant range for Pb, Zn, and Cd is 400–600 mm. In other words, As is highly concentrated in areas with low rainfall, while Pb, Zn, and Cd are highly concentrated in areas with high rainfall. Previous studies have shown that the abundance of hot springs on the TP influences the distribution of As in the plateau environment [52]. Thus, the dilution effect of rainfall makes the YTR show a spatial pattern of high As in upstream water and low As in downstream water [44]. Therefore, areas with low rainfall may be more prone to the enrichment of As. The YTRB is rich in mineral resources, and mining activities are relatively concentrated in the middle and lower reaches [53]. Abundant rainfall may promote the release of Pb, Zn, and Cd. Additionally, there is a significant correlation between soil Pb, Zn, and Cd in the YTRB (Figure 4b), and the high-value areas are mainly distributed in the middle and lower reaches (Figure 4a). Therefore, to reduce the long-term enrichment of Pb, Zn, and Cd, the tailings in the middle and lower reaches should be appropriately managed to avoid the infiltration of hazardous elements into the surrounding environment through precipitation. The primary type of other dominant factors of HMs, such as RT, GA, and sand, are also not entirely consistent (Figure 7a). These results can effectively locate the more vulnerable areas. They can be extrapolated to the entire TP, providing scientific support to benefit the health of the local population and to cope with the possible effects of future climate change.

The present model considers most of the factors related to soil HMs, but there are limitations in data representatives. For example, the compositional data closure effect and spatial variability of samples may have an impact on the possibility of our research results. In addition, the number of influencing factors stratified is based on previous research and data processing experience, which increases the uncertainty of the results. However, in general, the model of this study has high accuracy and reliable interpretation, which provides a scientific basis for soil management for environmental managers.

4. Conclusions

In this work, we analyzed the contamination data of six soil HMs (As, Pb, Cd, Cr, Ni, and Zn), combined with remote sensing data, geological information, and other geospatial data to identify the primary influencing factors of HM. The results for the I_geo showed that the pollution of As was the most serious. From the perspective of ecological risk, HMs are at low ecological risk, with Cd at the highest risk level, followed by As. The precipitation, pH, GA, and RT were the most important influencing factors for HMs. Based on the risk detector, the key influential types (ranges) of the predominant factors and their effects on the enrichment of the HM concentration were identified. This work provides new data on the extent of soil HM pollution on the TP. The results can be used to identify vulnerable areas of HM pollution in advance, providing an essential scientific basis for carrying out monitoring of HMs in remote areas at high altitude, and providing scientific support to safeguard the health of the local population and to cope with the possible effects of future climate change.