Soil Texture and Its Relationship with Environmental Factors on the Qinghai–Tibet Plateau

Abstract

Soil texture data are the basic input parameters for many Earth System Models. As the largest middle–low altitude permafrost regions on the planet, the land surface processes on the Qinghai–Tibet Plateau can affect regional and even global water and energy cycles. However, the spatial distribution of soil texture data on the plateau is largely unavailable due to the difficulty of obtaining field data. Based on collection data from field surveys and environmental factors, we predicted the spatial distribution of clay, silt, and sand contents at a 1 km resolution, from 0–5, 5–15, 15–30, 30–60, 60–100, and 100–200 cm soil depth layers. The random forest models were constructed to predict the soil texture according to the relationships between environmental factors and soil texture data. The results showed that the soil particles of the QTP are dominated by sand, which accounts for more than 70% of the total particles. As for the spatial distribution, silt and clay contents are high in the southeast plateau, and low values of silt and clay mainly appeared in the northwest plateau. Climate and NDVI values are the most important factors that affect the spatial distribution of soil texture on the QTP. The results of this study provide the soil texture data at different depths for the whole plateau at a spatial resolution of 1 km, and the dataset can be used as an input parameter for many Earth System Models.

1. Introduction

Soil texture is one of the most basic physical properties of soil [1,2]. The composition and content of coarse and fine particles can affect the hydraulic characteristics of water retention and conductivity in soil, as well as thermal parameters such as the heat transfer coefficient and heat capacity [3,4,5]. Soil texture can also affect evaporation and infiltration rates [6,7]. Soil texture has an important impact on the soil’s moisture. Studies have shown that the spatial heterogeneity of a soil’s texture is responsible for the spatial variation of moisture, which further increases the differences in the spatial distribution of moisture over time [8,9]. As a relatively stable natural attribute, soil texture also has a close relationship with soil carbon stocks [10,11,12,13]. Therefore, soil texture is an important input parameter for land surface process models, hydrological models, and atmospheric models [14,15]. In the warming climate, there is an increasing need for soil texture data in studies that deal with global and regional issues such as climate change, soil degradation, water resource shortage, environmental pollution, agriculture, and sustainable development.

The Qinghai-Tibet Plateau (QTP) has an average elevation of more than 4000 m and has been called the roof of the world [16]. The total area of the plateau is about 2.5 million km² [17]. As a typically cold region, the QTP hosts the largest permafrost region in the low and middle latitudes, with a permafrost area of 1.06 × 10⁶ km² [18,19]. The QTP is sensitive to global climate change, and increases in air temperature since the 1980s is at a rate of twice the global average [20,21]. As an important influencing factor and process parameter of the global water and energy cycle, soil moisture affects the distribution of surface energy, directly affecting the formation and development of soil [22]. In the past few decades, the permafrost on the QTP has been experiencing a rapid degradation, such as an increasing active layer, increasing ground temperature, and decreasing permafrost area [23,24,25]. Permafrost degradation not only presents threats to engineering construction [23,26], but also stimulates the decomposition of the soil organic matter and further affects ecological succession and the carbon cycle [21]. To quantify the environmental changes in the QTP and their feedback to climate, many land surface models have been used to simulate the changing processes, and have attempted to predict the future changes in the water, carbon, and energy cycles on the QTP. Most of these models require multiple input parameters, and soil texture data are one of the most important of these parameters.

Due to the importance of knowledge concerning soil texture data, many efforts have been conducted to create the spatial distribution of soil texture on the QTP. Previous studies have used spatial interpolation methods to upscale data from a site scale to a regional scale. The Food and Agriculture Organization of the United Nations (FAO) developed a global soil texture dataset with a data ratio of 1:100,000 [27], which has been widely used in many land surface models. This dataset contains very rare soil sample data in China, but there are almost no data in the QTP. According to the field survey data from the Second National Soil Survey of China, the conventional linkage method was used to obtain the spatial distribution of sand, silt, and clay in China, with a spatial resolution of 1 km [28]. The dataset compiled many observation data, and the accuracy was improved. However, the sampling sites on the QTP during the Second National Soil Survey were mainly distributed in seasonally-frozen ground areas, and the data in rural permafrost regions were largely lacking. Since the permafrost regions occupy nearly half of the QTP, there is much room for improvement in obtaining the spatial distribution of soil textures on the QTP.

Recently, several machine learning methods have been widely used for soil texture mapping. For example, the spatial distribution of soil texture fractions in Denmark, predicted at six depth layers at a resolution of 30 m, was performed using the Cubist decision tree algorithm [1]. A 100 m resolution model of the soil particle size distribution in Chile has been created using the classification and regression tree algorithm [29]. A soil texture map of different depths was also created for China, with a spatial resolution of 90 m, and used the random forest model [6]. However, there are still few field sites in the permafrost regions on the QTP.

In order to obtain the spatial distribution of soil texture data on the QTP, we collected a dataset across vast areas of both permafrost and nonpermafrost regions on the QTP. By combining high-resolution, remotely-sensed data and interpolated meteorological datasets, we simulated the spatial distribution of soil texture on the QTP (At depths of 0–5, 5–15, 15–30, 30–60, 60–100, and 100–200 cm) using the random forest method. We also examined the relationship between the soil texture and environmental factors. The results can serve as a basic input parameter for Earth System Models to study the processes of land surfaces on the QTP.

2. Materials and Methods

2.1. Study Area

The QTP (Figure 1) spans more than 3000 km from east to west and more than 1500 km from north to south [17]. The plateau is located in the semi-arid and arid alpine climate zone, with a mean annual precipitation of about 250 mm [30]. The mean annual average temperature on the QTP is below 2.8 °C [31]. The air temperature on the QTP shows great heterogeneities, with the coldest values in the north of the plateau and the warmest values in the southeast of the plateau [30].

2.2. The Observation Data

There were two sources for the data used in this study: (1) The field observation data of typical permafrost regions on the QTP from 2009 to 2012. The sampling sites were selected based on environmental gradients such as altitude, slope, aspect, and geophysical prospecting. During the field work, soil pits (a) were excavated and soil samples were collected according to the soil genetic horizons. (2) The data collected from the National Soil Series Survey and Compilation of Soil Series of China, conducted from 2009 to 2019 (Soil pits (b)). Most of these samples were located in the seasonally-frozen soil area of the QTP (Figure 1). The two sources of soil sample data mentioned above included valuable samples collected both in the permafrost area of the QTP and in seasonally frozen areas, which provide supporting data for the study of the spatial distribution of soil texture on the QTP. All of the soil profiles were recorded by the genetic horizons, and soil samples were collected at each horizon. The soil variables for different depths were resampled using the mathematical mean method. The soil properties were measured, including the soil texture, soil organic carbon, bulk density, pH, cation exchange capacity, total nitrogen, total phosphorus, total potassium, available phosphorus, and available potassium. Among them, we determined the soil texture fractions, namely the percentages of sand (2–0.05 mm), silt (0.05–0.002 mm), and clay (<0.002 mm), and these percentages were determined using a combination of the wet sieving technique and a particle size analyzer (Sedimat 4–12, UGT, Müncheberg, Germany) [32].

2.3. Data Processing and Analysis

As the measured particle composition is based on the occurrence layer, the number of layers divided by different profiles and the depth of each layer are both different, thus the depth of the acquired profile layers needs to be standardized (0–5, 5–15, 15–30, 30–60, 60–100, and 100–200 cm). The method is as follows:

$${\mathrm{PC}}_{\mathrm{t}–\mathrm{b}}={\displaystyle {\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{PC}}_{\mathrm{i}}\times {\mathrm{H}}_{\mathrm{i}}}$$

(1)

PC_t–b is the content of a certain fraction of the normalized level t–b, where t and b represent the top and bottom depths of the level, respectively; n represents the profile and the number of occurrence layers. PC_i represents the content of a specific fraction of the occurrence layer i, and H_i represents the proportion of the depth of layer i in the normalized layer t–b.

The soil particle size data were resampled at different depth intervals, i.e., at 0–5 cm, 5–15 cm, 15–30 cm, 30–60 cm, 60–100 cm, and 100–200 cm. Among the data, there were 316 samples for 0–5 cm and 5–15 cm, 314 samples for 15–30 cm, 302 samples for 30–60 cm, 266 samples for 60–100 cm, and 184 samples for 100–200 cm.

2.4. Environmental Covariates

The topography of the QTP is complex, meaning that the effects of environmental factors on soil texture are complicated. The soil is a product of climate, biology, topography, parent material, and time [33,34]. This theory provides a theoretical basis to select environmental factors to predict the soil texture, i.e., the climate, vegetation, topography, and parent materials. Although it largely remains unknown which climate variables are most important in affecting soil textures, the annual mean temperature and annual mean precipitation data are closely associated with soil textures. Therefore, we selected these variables as climate variable input parameters (

Table 1. Description of the selected variables.

Variables	Description	Resolution	Soil Forming Factors
Elevation	Elevation above sea level (km)	1 km	r
Slope	Slope gradient	1 km	r
Aspect	Aspect gradient	1 km	r
Plan	Plan curvature	1 km	r
Profile	Profile curvature	1 km	r
TWI	Topographic wetness index	1 km	r
TCA	Total catchment area	1 km	r
RSP	Relative slope position	1 km	r
CI	Convergence index	1 km	r
VD	Valley depth	1 km	r
CND	Channel network distance	1 km	r
CNB	Channel network base level	1 km	r
LS	Slope length and steepness factor	1 km	r
MAT	Annual mean temperature (°C)	1 km	c
MAP	Annual precipitation (mm)	1 km	c
NDVI	Mean NDVI during the growing season	1 km	o, c
QGT	Quaternary geological type	1 km	p
ST	Soil type	1 km	o

).

The DEM data were downloaded from the International Scientific Data Service Platform (http://datamirror.csdb.cn, accessed on 5 August 2021). Using the SAGA GIS software (http://www.saga-gis.org, accessed on 29 May 2022), we calculated the derivative variables of the DEM data (Table 1). Normalized difference vegetation index (NDVI) data were obtained from the United States Geological Survey (http://modis.gsfc.nasa.gov/, accessed on 10 August 2021). To avoid the possible errors that could be caused by rapid climate change since 2000, we used the multi-year (2000–2015) average data as the input parameters. The mean annual air temperature (MAT) and mean annual precipitation (MAP) data were downloaded from WorldClim version 2.1 (https://www.worldclim.org, accessed on 15 August 2021). These datasets were observed from global meteorological stations in the years 1970–2000. The parent material data were retrieved from the Quaternary geological map, which were download from the National Glacier, Permafrost and Desert Science Data Center (http://www.ncdc.ac.cn, accessed on 20 August 2021). Soil type data were retrieved from the Resource and Environment Science and Data Center (https://www.resdc.cn/, accessed on 22 August 2021).

2.5. Random Forest Model

During the past decades, various machine learning methods have been applied to digital soil mapping research. The classic machine learning method used is the neural network, despite being a computationally expensive method. The classification tree algorithm was established in the 1980s, and this method greatly reduced the amount of calculation needed. The random forest (RF) algorithm is a widely used machine learning algorithm based on the decision tree model [35,36,37]. The RF algorithm mainly consists of classification and regression models [38,39]. The RF algorithm has a high prediction accuracy, good tolerance to outliers and noise, and is not prone to overfitting [40,41]. Therefore, the algorithm is a natural nonlinear modeling tool, and it is currently one of the most useful tools in data mining and geographic information systems [42,43]. The random forest (RF) algorithm has been implemented into the R software (https://www.r-project.org, accessed on 22 February 2019). On the R platform, the random forest regression model uses the so-called ‘randomForest’ and ‘caret’ packages to study spatial distribution. The random forest regression model uses the %IncMSE index to describe the importance of the variables [44]. The greater the value of the %IncMSE, the greater the importance of the environmental variables in the model.

2.6. Accuracy Validation

In order to verify the random forest model, a ten-fold cross-validation was performed to evaluate its accuracy. The coefficient of determination (R²), mean absolute error (MAE), and root mean square error (RMSE) were used to assess the prediction’s accuracy. Higher R² and correspondingly lower MAE and RMSE values indicate a higher prediction accuracy. The calculation methods are as follows [6,45]:

$${\mathrm{R}}^2=1-\frac{{\displaystyle {\sum }_{\mathrm{i}=1}^{\mathrm{n}}{({\mathrm{P}}_{\mathrm{i}}-{\mathrm{O}}_{\mathrm{i}})}^2}}{{\displaystyle {\sum }_{\mathrm{i}=1}^{\mathrm{n}}{({\mathrm{P}}_{\mathrm{i}}-\overline{\mathrm{O}})}^2}}$$

(2)

$$\mathrm{MAE}=\frac{1}{\mathrm{n}}{\displaystyle {\sum }_{\mathrm{i}=1}^{\mathrm{n}}|({\mathrm{O}}_{\mathrm{i}}-{\mathrm{P}}_{\mathrm{i}})|}$$

(3)

$$\mathrm{RMSE}=\sqrt{\frac{1}{\mathrm{n}}{\displaystyle {\sum }_{\mathrm{i}=1}^{\mathrm{n}}{({\mathrm{P}}_{\mathrm{i}}-{\mathrm{O}}_{\mathrm{i}})}^2}}$$

(4)

The ${\mathrm{P}}_{\mathrm{i}}$ and ${\mathrm{O}}_{\mathrm{i}}$ are the predicted and observed values of a soil texture fraction at the sample site i; $\overline{\mathrm{O}}$ is the average observation value.

3. Results

3.1. Statistical Summary of Soil Texture Samples

For all of the soil particle data, the mean values of sand were significantly higher than those of silt and clay, and the values of silt and clay decreased while sand increased with depth (

Table 2. Statistical analysis of the clay, silt and sand percentages at different depths.

Depth (cm)	Mean (%)	SD (%)	CV	Skewness	Kurtosis
Clay
0–5	11.56	7.25	0.63	0.66	0.01
5–15	11.55	7.36	0.64	0.76	0.26
15–30	11.33	7.69	0.68	0.93	0.84
30–60	10.81	7.99	0.74	1.03	0.81
60–100	9.90	8.13	0.82	1.25	1.19
100–200	8.00	7.11	0.89	1.70	2.59
Silt
0–5	16.94	12.72	0.75	0.87	−0.08
5–15	16.71	12.16	0.73	0.71	−0.73
15–30	16.51	12.58	0.76	0.85	−0.33
30–60	15.56	12.19	0.78	0.89	−0.28
60–100	14.37	12.05	0.84	1.19	0.55
100–200	12.41	11.96	0.96	1.88	3.18
Sand
0–5	71.48	17.14	0.24	−0.60	−0.41
5–15	71.73	16.80	0.23	−0.49	−0.83
15–30	72.14	18.02	0.25	−0.63	−0.63
30–60	73.61	18.23	0.25	−0.77	−0.43
60–100	75.64	18.55	0.25	−0.99	0.29
100–200	79.58	17.15	0.22	−1.66	0.36

). The standard error (SD) values of clay were lower than those of sand and silt, ranging from 7.11 to 8.13. The SD values of silt and sand were 11.96–12.72 and 16.80–18.55, respectively. The coefficient of variation (CV) for sand ranged 0.22–0.25, which was a lower range than those of silt and clay. The skewness and kurtosis of silt, clay, and sand were all near zero, indicating that the soil particle distribution largely satisfies the normal distribution.

3.2. Model Accuracy Evaluation

The ten-fold cross-validation showed that the R² values of sand and silt were significantly higher than that of clay (

Table 3. Accuracy assessment of the soil particle size distribution predictions.

Depth (cm)	R2	RMSE	MAE
Clay
0–5	0.36	5.85	4.30
5–15	0.33	6.10	4.45
15–30	0.30	6.49	4.83
30–60	0.28	6.91	5.06
60–100	0.32	6.84	4.86
100–200	0.50	5.02	3.49
Silt
0–5	0.49	9.11	6.46
5–15	0.50	8.64	6.49
15–30	0.47	9.29	6.87
30–60	0.44	9.30	7.02
60–100	0.49	8.77	6.19
100–200	0.57	8.05	5.43
Sand
0–5	0.52	12.01	8.68
5–15	0.51	11.95	8.91
15–30	0.45	13.54	10.15
30–60	0.48	13.29	10.12
60–100	0.49	13.28	9.47
100–200	0.62	10.67	7.21

). The R² range of values for silt and sand were 0.44–0.57 and 0.45–0.62, respectively. This indicates that the random forest model can explain 44–57% of the total variances in silt and 45–62% of the total variance in sand. The simulation of sand and silt had higher accuracies than those of clay, which had R² values in the range of 0.28–0.50.

3.3. Spatial Patterns of Predictions

The contents of silt and clay in soil are largely high in the southeast and low in the northwest, while sand showed the opposite spatial patterns, with high values in the northwest and low values in the southeast (Figure 2, Figure 3 and Figure 4). Similar with the statistical results from the field observation data, the sand content increases with depth (Table 2).

3.4. Controlling Factors of QTP Soil Texture Patterns

The importance of environmental factors for soil texture in each layer were also calculated by the RF algorithm. We listed the relative importance of the first 12 important variables at multiple depths (Figure 5, Figure 6 and Figure 7). For the clay content in the soils, the most important environmental factors affecting the spatial distribution of clay for the depth intervals within 0–60 cm were the annual mean temperature, channel network base level, and NDVI. For the 60–100 cm, the most important environmental factors were the channel network base level, annual mean temperature, NDVI, and elevation. For the 100–200 cm depth, the channel network base level, elevation, and annual mean temperature were the most important factors. For the silt content in the soil, the annual precipitation and annual mean temperature were the two most important factors in each layer. In the depth layer intervals within 0–60 cm, soil type, channel network base level and NDVI were also relatively important in addition to annual mean temperature and annual precipitation. In the depth layer intervals within 60–200 cm, elevation also played an important role. For the sand content in the soil, the annual mean temperature, annual precipitation, and NDVI were the main environmental factors that affect the spatial distribution of sand. For the 100–200 cm layer, soil type and elevation showed the greatest impact in sand contents.

4. Discussion

Our study showed that the soil particles on the QTP were mainly composed of coarse particles of sand, which accounted for more than 70% of the total content, along with low contents of silt and clay. As for the spatial distribution, silt and clay were higher in the east of the plateau and lower in the west of the plateau. Sand showed an opposite spatial pattern. This finding is in agreement with previous reports that the soil on the QTP largely consist of coarse particles [32,46,47]. The QTP is the youngest plateau of the world, and the harsh natural conditions have hampered vegetation growth [48]. As a result, the weak pedogenesis leads to the coarse soil particle sizes distribution. Comparing with the permafrost map of the QTP [18], it can be found that higher sand contents mainly appear in the permafrost regions, while silt and clay are mainly distributed in seasonally-frozen ground regions. This can be explained as the fact that the permafrost regions on the QTP have a low air temperature and precipitation, but a high evapotranspiration [49,50]. Under these conditions, pedogenesis is usually weak [51], resulting in the coarse soil particle distribution that was observed.

In some studies that are based on regression models, collecting sample data in the study area is relatively easy, and can be obtained in large quantities. Because of this, a part of the sample data can be reserved for the accuracy evaluation of the mapping results and can be used in the comparative analysis of previous research results [6,52,53,54]. However, the QTP has a large altitude gradient and a harsh environment, especially in the permafrost regions, making data collection difficult. Therefore, our sample points were few and unevenly distributed relative to the huge QTP. We used all of the sample points of each standard layer for the establishment of the random forest model, and used a ten-fold cross-validation to evaluate the model’s accuracy. The ten-fold cross-validation showed that the random forest model has a good accuracy in the prediction of sand and silt. In many studies, the simulation for sand and silt usually had higher accuracies than that of clay [53,55,56]. This is reasonable, as the clay contents usually were low and had large heterogeneities in horizontal and vertical distribution. For the predictions of silt and sand, our results had R² values in the range of 0.44–0.57 and 0.45–0.62 for silt and sand, respectively, and these values were high. For the soil texture mapping of China, the ten-fold cross-validated R² values of sand and silt are 0.45–0.50 and 0.44–0.49, respectively [6]. Using the Cubist decision tree algorithm, the accuracy of soil texture predictions in Denmark had the R² values of 0.27–0.44, 0.29–0.55, and 0.26–0.54 for sand, silt and clay, respectively [1]. In France, the R² values of silt, clay, and sand predictions were are 0.21–0.42, 0.19–0.33, and 0.25–0.44 for sand, silt, and clay, respectively [53]. Through the horizontal comparison with the regression model R² of the above studies in different regions worldwide, we found that we achieved good results in the prediction and mapping of the soil texture stratification in the QTP using the random forest model.

For clay, within the depth layers of 0–30 cm, CNB, MAP, and NDVI are the main control environmental factors, but in the layers of 30–60 cm, 60–100 cm, and 100–200 cm, the main control environmental factors have no obvious regularity. For silt, within the depth layers of 0–30 cm, MAT, MAP, and ST are the main control environmental factors, and at 30–60 and 60–100 cm, MAT, NDVI, and MAP play the most important role. For sand, within the 0–30 cm depths, MAT, MAP, and NDVI play a decisive role, while at 30–60 and 60–100 cm, NDVI and MAT are the main control factors. Therefore, MAT, MAP, and NDVI can generally be regarded as the most important environmental factors affecting the spatial distribution of soil texture on the QTP. Studies have shown that the precipitation decreases from east to west, and that the temperature decreases when moving from southeast to northwest on the QTP [57,58], which has great synergy with the spatial distribution of soil texture that we predicted. Climate can control the formation and development of soil by affecting vegetation type [59]. Under arid and semi-arid climates, the QTP has a sparse soil vegetation coverage, and the surface soil is vulnerable to wind and water erosion, meaning that the organic matter accumulation and profile development are weak. Therefore, vegetation develops better in areas on the QTP with higher temperatures and more precipitation, which promotes the accumulation of organic matter and the development of soil profiles, and results in finer soil particles; however, in areas with low temperature and low precipitation, the vegetation development is very poor, the surface is mostly bare, the soil development degree is very low, and the soil particles are mostly composed of large particles. Studies have shown that the grassland, which occupies nearly 75% of the plateau, is the dominant land cover type of the QTP, and the vegetation coverage on the plateau decreases when moving from the southeast to the northwest [30,60]. Consequently, the NDVI was found to be an important factor affecting the distribution of soil particles. We stress that our study shows the importance of mean annual precipitation and mean annual air temperature in the determination of the soil texture; it is possible that other climate variables, such as monthly mean air temperature and precipitation, may be better predictors of soil texture. Future studies should be conducted to find more useful climate indices to increase the accuracy of soil texture predictions.

The QTP has long been recognized as an important geomorphic unit on Earth [61]. The energy, water, and carbon cycles in this region has global implications [21]. Affected by factors such as altitude, temperature, and precipitation, the composition of soil particles contains obvious regional characteristics. However, the existing soil texture data on this plateau prior to this study were only part of data on a national or global scale [6]. Previous reports only included a few sample sites from the QTP, and especially in permafrost regions where collection is limited by terrain, traffic, etc., sample sites were few. In addition, the models to predict soil particle size from a national or global scale likely result in a lower prediction accuracy on the QTP, which can be attributed to the fact that the QTP has very different climatic and geological conditions. Thus, the results obtained by machine learning have a great level of uncertainty in the QTP region. In this study, we collected a large amount of field soil profile data on the QTP to overcome the shortcomings of previous studies.

5. Conclusions

We collected field observation data on the QTP and examined the spatial distribution of the soil texture on the QTP using the random forest regression model. The ten-fold cross-validation results show that the simulation of silt and sand have high accuracy, with the R² values being 0.44–0.57 and 0.45–0.62 for silt and sand, respectively. On the QTP, the soil particles were dominated by sand, which accounted for more than 70% of the total composition. The sand contents were mostly high in the northwest and low in the southeast portion of the plateau. Conversely, clay and silt contents were high in the southeast and low in the northwest. The high sand contents were mainly found in permafrost areas with high altitudes and harsh climatic conditions. In contrast, the high values of silt and clay mainly appeared in the seasonally-frozen ground areas. The climate and NDVI values were determined to be the most important environmental factors affecting the spatial distribution of soil texture on the QTP. It is worth noting that it is difficult to obtain soil samples in the ‘no man’s land’ hinterland of the QTP due to the limited traffic access. Therefore, more field observation data studies are needed to improve the simulation accuracy at the plateau scale.