Quantitative Analysis of Factors Influencing Spatial Distribution of Soil Erosion Based on Geo-Detector Model under Diverse Geomorphological Types

Abstract

The Loess Plateau of China suffers from severe erosion, which results in a great variety of economic and ecological problems. For scientific control of soil erosion, the key questions urgently to be addressed in this paper are: (1) Which are the driving factors under diverse geomorphological types? (2) Do these driving factors operate independently or by interacting? (3) Which zones under diverse geomorphological types suffer from severe erosion and need more attention? In this paper, the RUSLE model was applied here to demonstrate the spatio-temporal variations in soil erosion from 2010 to 2017 in Yan’an City, and the Geo-detector model proved to be a useful tool to solve the questions mentioned above. The results showed that the average erosion modulus in Yan’an City decreased by 1927.36 t/km²·a from 2010 to 2017. The intensity of soil erosion in the northern Baota District, central Ganquan County, Luochuan County, Ansai County, and Zhidan County decreased to varying degrees. The effect size of driving factors affecting soil erosion varied significantly in diverse geomorphological types. The effect size of interaction between land-use types and vegetation coverage, land-use types and slope, slope and precipitation was higher than that of a single factor. High-risk zones with severe erosion were closer to cultivated land and forest land with steep slopes (>25°) in the mid-elevation hills of Yan’an City. Additionally, based on the specificity of the study area, the Geo-detector model performed better in a relatively flat region, and factors with macroscopic spatial distributions weaken its explanatory power on soil erosion on a regional scale. Based on data selection, data of different accuracy sparked the issue of “data coupling”, which led to the enormous deviation of model results in mid-elevation plains. Results from our analysis provide insights for a more ecologically sound development of Yan’an City and provide references for the scientific use of the Geo-detector model.

1. Introduction

Healthy soil is fundamental to life on Earth and the food system, as a medium that can provide human beings and the whole ecological system with basic materials, such as organic matters [1,2]. According to the Food and Agriculture Organization of the United Nations (FAO), the majority of the world’s soil resources are facing many challenges and threats [3]. The analysis results of soil sustainability have demonstrated that soil erosion is reported to be increasing, causing the main threat to healthy soil conditions across the globe [1,4,5].

Supported by numerous studies, soil erosion is a complex process under the influence of multiple factors and it is often the result of multiple driving factors and their interactions [6]. Extensive studies have identified that vegetation coverage [7,8], land-use types [9,10,11], precipitation [12,13], slope [14], and soil types [3,15] are important natural factors affecting soil erosion. Human activities in pursuit of rapid economic development lead to the over-exploitation of the soil environment and also to an imbalance in the human–landscape system, thus accelerating soil erosion [16,17]. Some authors concluded that the increase in vegetation coverage can effectively reduce from 98.0% to 99.9% of soil erosion, and erosive rainstorms can greatly intensify soil loss from 2.1 to 15.6 times [18]. Many studies have shown that the changes in land use or vegetation have a large impact on soil erosion [1,19,20]. In particular, human-induced land-use changes and their effects on soil erosion are being repeatedly mentioned as a threat to the environment and agricultural productivity [21]. To effectively reduce soil erosion, the following key and difficult points need to be answered: (1) Among these driving factors, which are the dominant factors, and what is the influential degree of each factor? (2) Do these influence factors operate independently or by interacting? (3) How do we identify the areas at higher risk of soil erosion?

With the popularization of computer technology, massive statistical analysis methods, including correlation and regression analysis, principal component analysis, and redundancy analysis [22,23,24,25,26], and spatial analysis methods, including hot-spot analysis, superposition statistical analysis, etc. [27], have been applied in the study of factors driving soil erosion. By estimating a multiple regression to investigate the effect of modified rainfall erosivity, land use and soil erodibility on soil erosion, the authors concluded that rainfall erosivity and soil erodibility were the major factors [15]. As an alternative to classical regression, network and redundancy analysis has also been applied as a practical method to quantify the interactions co-dependencies between soil erosion and influential factors as well as the interactions among driving factors [28]. Nonetheless, these statistical methods can be used to quantify the relationship between soil erosion and its driving factors, often failing to demonstrate the spatial distribution characteristics of those influencing factors [15,29]. Additionally, spatial analysis methods offer the advantages of identifying zones with high or low erosive intensity, but they cannot directly quantify the effect of factors driving soil erosion [30]. Further efforts are still needed to solve the three questions mentioned above.

The Geo-detector model is a relatively new statistical approach [31,32,33]. It is based primarily on the clear assessment of spatial heterogeneity in geographical phenomena. It characterizes the similarities in spatial variations between geographic data. In principle, if two independent variables have a significant effect on a dependent variable, then the spatial distribution of the two should be similar. The model includes factor detection, interactive factor detection, risk detection, and ecological detection. It can not only quantitatively determine the dominant factors but can quantify the influence of two interacting explanatory variables on a specific target variable. Specifically, it can also detect high-risk areas of soil erosion in the study area by superimposing the maximum soil erosion modulus of each influencing factor [34,35]. In theory, the model can answer all the questions discussed above.

The Loess Plateau in China exhibits some of the most severe soil erosion problems in the world. According to the “Comprehensive Scientific Investigation of Soil Erosion and Ecological Security in China”, the Loess Plateau has a total area of 640,000 km², of which 390,000 km² was characterized by eroded soils [36]. Yan’an City, located in the hinterland of the Loess Plateau, is one of the most severely water-eroded areas in the middle reaches of the Yellow River Basin; its annual average soil erosion modulus was 9000 t/km²·a [37]. In this paper, the RUSLE model was used to estimate the temporal and spatial variations in soil erosion from 2010 to 2017 and the Geo-detector model was applied to determine the dominant factor(s) and their interactions driving soil erosion, as well as to identify the high-risk areas of soil erosion based on different geomorphological types. The study will provide references for ecological protection and sustainable development in the study area.

2. Materials and Methods

2.1. Study Area

The landscape of Yan’an City consists of a combination of ridges, plateaus, and mounds. The northern and western parts of the city are located at relatively high elevations and are characterized by hilly gullied areas. The southern and eastern parts of the city are lower in elevation and consist of plains. Soil types mainly include loessal soil, the clastic and predominantly silt-sized sediment that is formed by the accumulation of wind-blown dust and brown soil. The region is subject to arid to semi-arid climatic conditions with a mean annual maximum temperature of 22.9 °C and annual precipitation of 400–500 mm. Summer storms are common and can deliver over 40 mm of precipitation in a single event. Yan’an City has serious hydraulic erosion problems over a wide area, and the erosion is the main source of sediment to the lower Yellow River (Figure 1).

2.2. Statistics of Geographical Factors among Various Geomorphological Types

A reference study divided the geomorphology of China into 28 basic landscape types based on altitude and topographic relief [38]. Yan’an City was classified into six specific types according to elevation, including mid-elevation plains, mid-elevation terraces, mid-elevation hills, low-elevation mountains, and mid-elevation mountains. The mountains were further classified into areas of small and medium relief according to undulations in the topography. The geographic environmental factors in various geomorphological types (Figure 1c) were spatially heterogeneous. The average slope, elevation, and terrain niche index, an index that combines elevation and slope, were all high in mountainous and hilly areas. The annual mean precipitation varied little between the different geomorphological types. Cultivated land was mainly concentrated in small-relief and low-elevation mountains, mid-elevation plains, mid-elevation terraces, and mid-elevation hilly areas (

Table 1. Statistics of geographical factors in different geomorphological types.

Geomorphological Types	Average Elevation (m)	Average Slope (°)	Average Terrain Niche Index	Annual Mean Precipitation (mm)	Proportion of Cultivated Land in 2017 (%)
Small-relief and low-elevation mountains	766.74	13.27	0.32	617.60	0.40
Mid-elevation plains	1066.64	12.41	0.40	623.90	0.38
Mid-elevation terraces	1071.03	13.88	0.41	634.13	0.38
Small-relief and mid-elevation mountains	1205.90	15.44	0.46	640.30	0.23
Mid-elevation hills	1265.97	16.17	0.48	621.15	0.31
Mid-relief and mid-elevation mountains	1366.39	16.88	0.51	626.18	0.12

).

2.3. Data Sources

The daily precipitation data were obtained from the China Meteorological Data Network (http://data.cma.cn/, accessed on 1 June 2021). The UN Food and Agriculture Organization and the International Institute for Applied Systems Analysis provided soil mechanical composition data with a spatial resolution of 1 km (http://westdc.westgis.ac.cn/, accessed on 1 June 2021). DEM data at a spatial resolution of 30 m were derived from the geospatial data cloud on the computer network information center of the Chinese Academy of Science (CAS) (http://www.gscloud.cn/, accessed on 1 June 2021). MOD13Q1 data products at a 250 m spatial resolution were obtained for 2010 and 2017 from the National Aeronautical and Space Administration (NASA) (https://www.nasa.gov/, accessed on 1 June 2021). The land-use and coverage data at a 30 m resolution in 2010 were gathered from the National Geographic Information Public Service Platform (NGIPSP) (http://www.globelandcover.com, accessed on 1 June 2021), whereas 10 m resolution data for 2017 were obtained from a database of global land-use coverage created by Tsinghua University (http://data.ess.tsinghua.edu.cn, accessed on 1 June 2021). The soil texture data, soil type data, geomorphic distribution data, population data, and vector boundary data were obtained from the Resources and Environmental Science Data Center of the Chinese Academy of Sciences (http://www.resdc.cn/, accessed on 1 June 2021).

2.4. Research Methods

2.4.1. RUSLE Model

The Revised Universal Soil Loss Equation (RUSLE) [39] was applied to estimate soil erosion modulus in Yan’an City. The mathematical expression is as follows:

$$\mathrm{A}=\mathrm{R}·\mathrm{K}·\mathrm{L}·\mathrm{S}·\mathrm{C}·\mathrm{P}$$

(1)

where A is the estimated soil erosion modulus (t/km²∙a), R is the rainfall-runoff erodibility factor (MJ·mm·km⁻²·h⁻¹·a⁻¹), K is the soil erodibility factor (t·km²·h·MJ⁻¹·mm⁻¹·km⁻²), LS is the slope length factor, C is the coverage management factor, and P is the soil and water conservation measures factor. The range of C and P values is [0, 1].

(1)

Rainfall-runoff erodibility factor (R):

The R factor is the external influencing factor driving soil erosion [13]. The classical method of EI₃₀ (energy and 30 minutes of intensity of rainfall) and simple algorithms for conventional meteorological data were often used to calculate R [40], and both methods have limitations in terms of data acquisition. We adopted a simple algorithmic model based on daily rainfall amount and the equation is as follows [41]:

$$\mathrm{R}=\propto {\displaystyle \sum }_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{p}}_{\mathrm{j}}^{\mathsf{\beta }}$$

(2)

$$\mathsf{\alpha }=21.586{\mathsf{\beta }}^{-7.1891}$$

(3)

$$\mathsf{\beta }=0.8363+18.144{\mathrm{P}}_{\mathrm{d}12}^{-1}+24.455{\mathrm{P}}_{\mathrm{y}12}^{-1}$$

(4)

where P_j is the precipitation in half a month (mm). P_d12 is the average daily precipitation with a daily rainfall over 12 mm. P_y12 is the average annual precipitation with a daily rainfall over 12 mm. α and β are the model parameters.

(2)

Soil erodibility factor (K):

The K factor is a parameter to characterize the sensitivity of soil properties to erosion. We adopted Erosion-Productivity Impact Calculator (EPIC) model and the equation is as follows [42]:

$$\mathrm{K}=\left\{0.2+0.3{\mathrm{e}}^{\left[-0.0256\mathrm{SAN}\left(1-\frac{\mathrm{SIL}}{100}\right)\right]}\right\}\times {\left(\frac{\mathrm{SIL}}{\mathrm{CLA}+\mathrm{SIL}}\right)}^3\times \left(1-\frac{0.25\mathrm{C}}{\mathrm{C}+{\mathrm{e}}^{=2.95\mathrm{C}+3.72}}\right)\times \left(1-\frac{0.7\mathrm{SN}}{\mathrm{SN}+{\mathrm{e}}^{-5.51+22.9\mathrm{SN}}}\right)$$

(5)

where SAN, SIL, CLA, and C are the sand fraction (%), silt fraction (%), clay fraction (%), and content of soil organic carbon (%), respectively (SN=1˗SAN/100).

(3)

Slope length factor (LS):

The LS factor is the immediate trigger driving soil erosion [13]. L and S were calculated based on the methods and the equations are as follows [43]:

$$\mathrm{L}={\left(\mathsf{\lambda }/22.13\right)}^{\mathsf{\alpha }}$$

(6)

$$\mathsf{\alpha }=\mathsf{\beta }\left(\mathsf{\beta }+1\right)$$

(7)

$$\mathsf{\beta }=\left(\sin \mathsf{\theta }/0.0896\right)/\left[3\times {\left(\sin \mathsf{\theta }\right)}^{0.8}\right]+0.56$$

(8)

$$\mathrm{S}=\{\begin{array}{cc}10.8\times \sin \mathsf{\theta }+0.036& (\mathsf{\theta }<5.1428)\\ 16.8\times \sin \mathsf{\theta }-0.5& (5.1428\le \mathsf{\theta }<14.0362)\\ 21.9\times \sin \mathsf{\theta }-0.96& \mathsf{\theta }\ge 14.0362\end{array}$$

(9)

where λ is slope length (m). α is a variable length-slope exponent. β is a factor relative to slope gradient. θ is the slope.

(4)

Coverage management factor (C):

The C factor represents the effects of vegetation coverage and agricultural management. Based on MOD13Q1 data, the vegetation coverage (f) was extracted and C was calculated by Equation (10).

$$\mathrm{C}=\{\begin{array}{cc}1& 0\le \mathrm{f}<0.096\\ 0.6508-0.3436{\log }_{10}\left(\mathrm{f}\right)& 0.096\le \mathrm{f}<0.783\\ 0& \mathrm{f}>0.783\end{array}$$

(10)

$$\mathrm{f}=\left(\mathrm{NDVI}-{\mathrm{NDVI}}_{\min }\right)/\left({\mathrm{NDVI}}_{\max }-{\mathrm{NDVI}}_{\min }\right)$$

(11)

where f is vegetation coverage. NDVI is the Normalized Difference Vegetation Index. NDVI_min is the value of an area completely covered with bare land. NDVI_max is the value of an area completely covered with vegetation.

(5)

Soil and water conservation measures factor (P):

The P factor is defined as the ratio of soil loss after taking soil and water conservation measures to that of planting along the slope. It is closely related to the effect of land-use types and the P values of different land-use types are shown in the following tables (

Table 2. P values of different land-use types.

Land-Use Types	P
Forest land	1
Grassland	1
Water	0
Urban land	0
Unused land	1

and

Table 3. P values of cultivated land in different slopes.

Slope (°)	P
0 ~ 5	0.100
5 ~ 10	0.221
10 ~ 15	0.305
15 ~ 20	0.575
20 ~ 25	0.705

) [44].

Due to different data sources and spatial resolutions, all these factors in the RUSLE model were resampled to a 250 m spatial resolution.

2.4.2. Geo-Detector Model

The Geo-detector model has two main advantages over other spatial analysis methods. First, it can utilize both quantitative and qualitative data. Second, it can determine the influence of two interacting explanatory variables on a specific target variable. The model outcome includes factor detection, detection of interactions, risk detection, and ecological detection [31,32,33]. To answer the three questions mentioned above, the first three detection methods were used.

Factor detection: This is used to determine to what extent do independent variables (X) affect the dependent variable (Y). The degree of explanatory power is generally expressed in statistical terms as q-value, and it is calculated as follows:

$$\mathrm{q}=1-\frac{\mathrm{SSW}}{\mathrm{SST}}, \mathrm{SSW}={\sum }_{\mathrm{h}=1}^{\mathrm{L}}{\mathrm{N}}_{\mathrm{h}}{\mathsf{\sigma }}_{\mathrm{h}}^2, \mathrm{SST}={\mathrm{N}\mathsf{\sigma }}^2$$

(12)

where h is the index to denote each of the strata related to the dependent variable Y as well as the independent variables X; N_h and N are the units of stratum h and whole areas, respectively, and similarly in the following for the variances of σ_h² and σ²; SSW and SST are the sum of intra-layer variances and the total variances of whole areas, respectively; q-value is [0, 1], and the greater the q-value, the stronger the explanatory power of the independent variables X to the dependent variable Y.

Detection of interactions: This estimates whether the interaction of any two independent variables (X1 and X2) enhances or weakens their respective explanatory power of the dependent variable (Y). By comparing the q-value between a single factor (q (X1) and q (X2)) and paired factors (q (X1∩X2)), interaction types of this detection can be put into 5 categories (

Table 4. Types of interaction between two covariates.

Criterion	Interactive Forms
q(X1∩X2) < Min(q(X1),q(X2))	Weakened, nonlinear
Min(q(X1),q(X2)) < q(X1∩X2) < Max(q(X1),q(X2))	Weakened, single factor nonlinear
q(X1∩X2) > Max(q(X1),q(X2))	Enhanced, double factors
q(X1∩X2) = q(X1) + q(X2)	Independent
q(X1∩X2) > q(X1) + q(X2)	Enhanced, nonlinear

).

Risk detection: This can identify the high-risk regions of soil erosion based on stratification differences of factors.

The input variables to the Geo-detector model are required to be stratified; thus, continuous variables need to be discretized (Figure 2). Additionally, the model must connect dependent variables and independent variables at specific locations or points. A total of 69,088 points collected at a sampling interval of 1 km served as the operational data for the Geo-detector model in ArcGIS. The soil erosion map in 2017 was used to run the Geo-detector model.

Note: vegetation coverage (X1), terrain niche index (X2), precipitation (X3), slope (X4), the population density (X5), land-use types (X6), soil types (X7), soil texture (X8).

3. Results

3.1. Spatio-Temporal Changes in Soil Erosion Intensity from 2010 to 2017

According to the classification standards of soil erosion issued by the Ministry of Water Resources, the mean erosion modulus (t/km²·a) within the study area can be uniformly divided into six categories: slight (A ≤ 1000), minor (1000 < A ≤ 2500), moderate (2500 < A ≤ 5000), intense (5000 < A ≤ 8000), very intense (8000 < A ≤ 15,000), and extreme (A ≥ 15,000) (

Table 5. The proportion of study area characterized by soil erosion intensity categories.

Erosion Intensity Level	2010 (Year)	2017 (Year)
Slight (%)	38.44	54.53
Minor (%)	5.63	4.13
Moderate (%)	7.92	5.83
Intense (%)	7.12	4.91
Very intense (%)	13.61	9.83
Extreme (%)	27.34	20.84
Mean erosion modulust (/km2·a)	9246.57	7319.21

).

The average erosion modulus of Yan’an City decreased from 9246.57 t/km²·a in 2010 to 7319.21 t/km²·a in 2017; the overall decrease was 1927.36 t/km²·a. The proportion of the study area classified by erosion intensity level as “Slight” increased by 16.09%, whereas the proportion related to other erosion intensity levels (from “Minor” to “Extreme”) decreased by 1.50%, 2.09%, 2.21%, 3.78%, and 6.50%, respectively (Table 5). The mean erosion modulus among various geomorphological types decreased to varying degrees from 2010 to 2017. The obvious result was that the proportion of “Slight” increased and that of “Extreme” decreased in each geomorphological type (Figure 3). The observed trends indicated that the degree of soil erosion in Yan’an City has appreciably reduced.

In 2010, areas of “Slight” erosion were distributed mainly in the southern part of the study area and occurred along distinct slices or belts of terrain. Importantly, the belt-like distribution was concentrated within cultivated land with gentle slopes, including southern Luochuan County, central Fuxian County, and central Ganquan County. Areas of “Extreme” erosion were mainly distributed in the northern and eastern parts of the study area. These areas also formed “belts” where erosion intensity increased with slope (Figure 4a,c). By 2017, the spatial distribution of soil erosion in Yan’an has changed dramatically. The areas of “Slight” covered a large part of the study area and had spread northward. Areas of “Intense, Very Intense, and Extreme” erosion were greatly reduced. In short, the intensity of soil erosion in Luochuan County, Baota District, Ganquan County, Ansai County, and Zhidan County greatly declined between 2010 and 2017 (Figure 4b).

3.2. Geodetector-Based Quantitative Analysis of Soil Erosion Heterogeneity among Various Geomorphological Types in 2017

3.2.1. Analysis of Simple Effect by Influential Factors of Soil Erosion

Based on factor detection, the comprehensive results demonstrate that the major factors controlling soil erosion, in order of q-values were: vegetation coverage (79.42%) > land-use types (52.43%) > population density (22.78%) > precipitation (21.70%) > slope (14.73%) > soil types (10.90%) > soil texture (9.71%) > terrain niche index (6.55%), and that the effect degree of influencing factors driving soil erosion was notably different in areas of different geomorphological types. In small-relief and low-elevation mountains, the type of land-use was the dominant factor. In mid-elevation plains, land-use types, slope and population density were the dominant factors. In mid-elevation terraces, vegetation coverage, population density and land-use types were the dominant factors. In small-relief and mid-elevation mountains, vegetation coverage was the dominant factor. In mid-elevation hills, vegetation coverage and land-use types were the dominant factors. In mid-relief and mid-elevation mountains, vegetation coverage, precipitation and land-use types were the dominant factors (Figure 5). The explanatory power of vegetation coverage gradient increased with elevation, whereas a “turnpoint” was observed in mid-elevation hills, down to 14.20%. The following patterns were observed for the explanatory power of those factors, including the population density, slope, and terrain niche index. The explanatory power of those factors was higher in relatively flat areas than that in mountainous and hilly areas, but it was generally low. Precipitation had a general low explanatory power affecting soil erosion among various geomorphological types, except in mid-relief and mid-elevation mountains (13%).

3.2.2. Analysis of Interaction Effect of Factors of Soil Erosion

Results of detection of interactions show that the combined influence of any two factors was significantly higher than that of a single factor and that the dominant interaction differed between areas of different geomorphological types. This paper summarized the interactions with the five highest effect degrees (

Table 6. Interactive determination of dominant factors under different geomorphological types.

	Small-Relief and Low-Elevation Mountains	Mid-Elevation Plains	Mid-Elevation Terraces	Small-Relief and Mid-Elevation Mountains	Mid-Elevation Hills	Mid-Relief and Mid-Elevation Mountains
Dominant interaction1 (q-value)	X1∩X6 (0.1257)	X5∩X6 (0.2157)	X1∩X6 (0.2300)	X1∩X6 (0.1893)	X1∩X6 (0.1945)	X1∩X6 (0.3207)
Dominant interaction2 (q-value)	X6∩X7 (0.1195)	X4∩X5 (0.1939)	X1∩X4 (0.1839)	X1∩X4 (0.1737)	X1∩X5 (0.1631)	X1∩X3 (0.3079)
Dominant interaction3 (q-value)	X4∩X6 (0.1138)	X4∩X6 (0.1912)	X5∩X6 (0.1721)	X1∩X8 (0.1733)	X1∩X3 (0.1581)	X1∩X7 (0.2957)
Dominant interaction4 (q-value)	X5∩X6 (0.1107)	X3∩X6 (0.1825)	X1∩X5 (0.16)	X1∩X5 (0.1723)	X1∩X4 (0.1542)	X1∩X4 (0.2936)
Dominant interaction5 (q-value)	X6∩X8 (0.1059)	X1∩X4 (0.1714)	X1∩X7 (0.1661)	X1∩X7 (0.1721)	X1∩X2 (0.1474)	X1∩X8 (0.2934)

Note: Vegetation Coverage (X1), Terrain niche index (X2), Precipitation (X3), Slope (X4), The population density (X5), Land-use types (X6), Soil types (X7), Soil texture (X8).

). The higher the effect degree of any two factors, the more obvious the inter-layer difference between the two factors. The highest explanatory power in areas of all the geomorphological types, except for the mid-elevation plains, was the interaction between vegetation coverage (X1) and land-use types (X6). This meant that the inter-layer difference of soil erosion varied significantly between X1 and X6. For example, the amount of soil erosion between unused land with low vegetation coverage and forest land with high vegetation coverage differed greatly. The top five interactions were mostly combinations of land-use types (X6) and other factors in small-relief and low-elevation mountain areas. Among the terraces, mountains and hills, the top five interactions were mostly combinations of vegetation coverage (X1) and other factors, but there were obvious differences. Except for small-relief and mid-elevation mountains, the combination of slope (X4) and land-use types (X6) had a relatively high explanatory power (over 10% among various geomorphological types and as high as 19.12% in the mid-elevation plains). The significantly higher explanatory power of the interaction of slope (X4) and precipitation (X3) was observed in mid-elevation plains (the least undulating platforms) and mid-relief and mid-elevation mountains (the most undulating platforms) than other geomorphological types.

3.2.3. Identification of High-Risk Areas of Soil Erosion

Based on the principle of risk detection, stratification of influencing factors by high-risk of soil erosion (soil erosion modulus ≥ 15,000 t/km²·a) can be determined (

Table 7. High-risk stratification of factors in different geomorphological types and the average erosion modulus (t/km²·a).

	Small-Relief and Low-Elevation Mountains	Mid-Elevation Plains	Mid-Elevation Terraces	Small-Relief and Mid-Elevation Mountains	Mid-Elevation Hills	Mid-Relief and Mid-Elevation Mountains
Vegetation Coverage	0.5–0.75	0.5–0.75	0.5–0.75	0.5–0.75	0.1–0.75	0.5–0.75
Average	11,549.73	11,179.80	11,236.70	10,103.10	15,909.79	12,207.72
Terrain niche index	0.2–0.6	0.2–0.8	-	-	0.6–0.8	-
Average	10,699.89	15,769.80	-	-	12,445.71	-
Precipitation	590–620	590–660	500–530, 590–620	560–590	560–590	560–590
Average	10,715.86	10,853.52	9522.95	8756.87	12,164.26	12,141.47
Slope	>6°	>15°	>15°	-	>25°	-
Average	12,506.80	14,428.36	9890.23	-	11,386.47	-
population density	0–200	250–350	50–100	200–250	50–100	-
Average	10,900.69	11,432.38	11,971.24	18,463.1	13,155.69	-
Land-use types	Unused land, Forest land, Grassland	Unused land, Forest land, Grassland	Unused land, Grassland	-	Unused land, Cultivated land, Grassland	Unused land
Average	15,057.63	12,437.11	17,999.81		19,611.94	12,796.70
Soil types	Skeletal soil, Red clay, Loessal soil, Dark loessial soil	Skeletal soil, Red clay, Alluvial soil	Skeletal soil, Red clay, Loessal soil	-	Skeletal soil, Red clay, Loessal soil	-
Average	11,047.01	11,361.38	10,060.6	-	8730.23	-
Soil texture	Loam clay; Clay loam; Sandy clay loam; Loam; Sandy loam	Sandy clay loam, Loam	Sandy clay loam, Sandy loam	-	Loam clay; Sandy clay loam; Loam; Sandy loam	-
Average	10,040.31	9855.87	9367.34	-	9986.68	-

Note: - means that soil erosion modulus is less than 15,000 t/km²·a.

, and high-risk areas of soil erosion can be identified (Figure 6). In all geomorphological types, vegetation coverage was at a critical value in the range of 0.5–0.75, where the soil erosion modulus peaked and then decreased. The soil erosion modulus increased with the slope in all geomorphological types, except for small-relief and mid-elevation mountains and mid-relief and mid-elevation mountains. The land-use types prone to erosion were unused land, forest land and grassland. The soil types prone to erosion were skeletal soil and loessal soil, both of which are soft and loose. The soil texture types prone to erosion were sandy clay loam and loam with sand percentages of 55–85% and 40–55%, respectively (Table 7). The high-risk areas of soil erosion were mainly distributed in the northern and eastern edges of the study area, mainly located in the cultivated land (62.57%) and forest land (26.06%) in mid-elevation hills (Figure 6) (

Table 8. The proportion of land-use types in high-risk areas of different geomorphological types.

Geomorphological Types	Land-Use Types	Proportion (%)
Small-relief and low-elevation mountains	forest land	2.08
	grassland	1.31
	unused land	0.48
Mid-elevation plains	forest land	0.31
	grassland	0.32
	unused land	0.27
Mid-elevation terraces	cultivated land	0.29
	grassland	0.37
Mid-elevation hills	forest land	26.06
	cultivated land	62.57
Mid-relief and mid-elevation mountains	grassland	5.11
	unused land	0.84

). The rehabilitation regions of small-relief and low-elevation mountains were concentrated in forest land (2.08%), grassland (1.31%) and unused land (0.48%) with steep slopes (>6°) and precipitation (590–620 mm). Regions of forest land (0.31%), grassland (0.32%) and unused land (0.27%) with steep slopes (>15°) and precipitation (590–660 mm) were focal areas of soil erosion prevention in mid-elevation plains. The rehabilitation regions of mid-elevation terraces were cultivated land and grassland with steep slopes (>15°). The rehabilitation regions of mid-relief and mid-elevation mountains were grassland and unused land with precipitation (560–590 mm).

3.3. Temporal Analysis of Dominant Factors of Soil Erosion

Soil erosion and its dominant factors varied temporally. To explore the changes in the explanatory power of controlling factors, the average erosion modulus was defined as the dependent variable, and vegetation coverage, precipitation, population density, and land-use types were considered the independent variables (

Table 9. Changes in the explanatory power (q-value) of time-dependent factors between 2010 and 2017.

	2010 (Year)	2017 (Year)
Vegetation Coverage	0.0290	0.1460
Precipitation	0.0050	0.1290
Population density	0.0380	0.0180
Land-use types	0.0300	0.0310

). Factor detection was used to explore their effect degree variation on soil erosion. Except for the population density factor, the explanatory power of the other influencing factors increased. Additionally, the effect size of vegetation coverage was always higher than other factors. The land coverage change patterns in Yan’an City from 2010 to 2017 (

Table 10. Changes in the proportion (%) of land-use types between 2010 and 2017.

	2010 (Year)	2017 (Year)
Proportion of cultivated land	24.39	26.27
Proportion of forest land	38.82	46.20
Proportion of grassland	35.71	17.92
Proportion of water area	0.22	0.22
Proportion of urban land	0.62	1.80
Proportion of unused land	0.08	6.57

) showed that the area of grassland decreased by 17.19%, while forest land, unused land and construction land increased by 7.38%, 6.49% and 1.18, respectively. These changes indicated that forest land and grassland were negatively related to severe soil erosion, while unused land was positively related to soil erosion. Precipitation also varied through time.

4. Discussion

4.1. Analysis of Controlling Factors of Soil Erosion

The effect degree of influencing factors, including vegetation coverage, land-use types, slope, terrain niche index, precipitation, population density, soil types and soil texture, and their interaction were analyzed quantitatively. We found that the interactions between vegetation coverage, land-use types and slope had a higher effect degree on soil erosion in each geomorphological type. Many scholars have noted that the effects of various land-use types on soil erosion differed significantly, and the soil erosion rate increased significantly as the slope of the terrain increased [1,9,45]. These conclusions emphasized the necessity of increasing forest coverage and constructing engineering measures on steep slopes. Compared to 2010, the intensity of erosion in most parts of central Yan’an City in 2017 has been greatly reduced, but soil erosion in the northern and eastern parts of the city remained relatively severe (Figure 4 and Figure 6). In recent years, two measures have been implemented to reduce soil erosion in the Loess Plateau. First, biological measures (mainly afforestation) have been applied to control soil erosion on slopes composed of loess. Second, engineering measures, such as the construction of silt dams, have been used to alleviate soil erosion within channels. Afforestation and restoration efforts using natural vegetation mainly reduce soil erosion by increasing vegetation coverage, whereas the development of terraces mainly controls soil erosion on sloping cultivated land by reducing slope length and by trapping silt [46,47]. Interventions to alleviate soil erosion have focused on vegetation and slope gradients in the study area. However, the tolerance of loessal soil to the continuous planting of forests and the cultivation of land, as well as other factors, especially soil properties of the study area, that are difficult to change over a short period, were ignored. Therefore, numerous studies and projects are needed to prevent the soil erosion of the loess gullies in the study area.

4.2. Model Validation

4.2.1. RUSLE Model

To carry out the validation of the RUSLE model, the average erosion modulus in this paper was compared with previously published results. Authors of a similar study estimated the average erosion modulus of 1360.11 t/km²·a in the Yangou watershed of Loess plateau [48]. Recently, another study suggested that the average erosion moduli in the Chinese Loess Plateau in 1999, 2000 and 2011 were 12,270.11, 13,050.11 and 11,460.11 t/km², respectively [49]. The average erosion modulus in the study area in 2010 and 2017 was 9246.57 and 7319.21 t/km², respectively (Table 5), which was within reasonable bounds. Therefore, the results of the RUSLE model in the study area show good reliability.

However, the RUSLE model has general validity and fails to take gravity erosion into consideration [18]. Gravity erosion mainly occurs on steep slopes with soft soil in hilly and gully regions and there was an underestimation of erosion modulus in the north of Yan’an City. The modification of the RUSLE model is needed to suit the hilly and gully regions in future research.

4.2.2. Geo-Detector Model

Based on calculation results in Section 3.2.1, the selection of the Geo-detector model does not take the specificity or spatial scale of the study area into account, resulting in a low effect degree of the terrain factors or the factors having microscopic spatial distribution, which is consistent with the findings of previous studies [50]. This finding is related to the limitations of the Geo-detector model, which relies on clear spatial differences between geographical phenomena [32]. However, the topography of the land in the study area is heterogeneous with areas of hilly and gully terrain, and its degree of fragmentation is higher. Some authors found that gully density in Ganquan County and Yanchuan County was as high as 5.61 and 6.78 km/km², respectively [51]. The difficulty of extracting clear spatial variations in slope in the unusually fractured topographic characteristics using DEM data at 30 m resolution weakens the explanatory power of the terrain factors on soil erosion. Moreover, precipitation data with the macroscopic spatial distribution and less microscopic differences brought by fewer sampling points, as well as the interpolation techniques adopted, weaken the explanatory power of this factor on soil erosion. Based on the discussion above, we concluded that this model works better in relatively flat regions and factors with too macro spatial distribution are more suitable for large spatial scales, such as nationwide or global scale, but not for regional scale.

The “data coupling problem” has become an important issue in academic circles due to the superposition of data with different accuracy. The statistical results that the erosive regions in mid-elevation plains focused on unused land, forest land and grassland with a steep slope (>15) (Table 7) differed from the field survey showing that urban land and cultivated land with gentle slopes and alluvial soil suitable for cultivation were the main land-use types in mid-elevation plains. The average erosion moduli of urban land and cultivated land were 0 and 5768.23 t/km²·a, respectively, and the overall erosion degree was relatively lower (

Table 11. Average erosion modulus of different land use types in mid-elevation plains (t/km²·a).

Land-Use Types	Average Erosion Modulus
Cultivated land	5768.23
Forest land	9709.05
Grassland	11,895.08
Unused land	14,123.74
Water area	0
Urban land	0

). Data accuracy inconsistency, such as the geomorphological types at a resolution of 1000 m and slope at a resolution of 30 m, led to errors in the extraction of factors, which resulted in inaccurate results of the Geo-detector model. Thus, the model should consider data accuracy consistency. In summary, the model is still of great significance in the application in soil erosion studies, but it also has a certain scope of application. All related issues mentioned above need further improvements and in-depth discussions in future studies.

5. Conclusions

This study analyzed the temporal and spatial changes of erosion intensity from 2010 to 2017, investigated the effect degree of factors and their interactions affecting soil erosion, and identified risk regions of each geomorphological type. The results are:

(1)

From 2010 to 2017, soil erosion in parts of the central counties within Yan’an City greatly improved, while in the northern and eastern regions, soil erosion remained a serious problem.

(2)

The analysis of time-dependent factors indicated that an increase in forest land can effectively improve soil erosion. Thus, transforming unused land in hilly gully regions to the forest and grassland is a feasible strategy to achieve sustainable development. It is worth mentioning that, in recent years, most people have a deeper understanding of the concept of environmental protection and urban planners are becoming more concerned about it. This means that as urbanization increased, the aggregation of the population makes economic and productive activities more efficient, which is beneficial to the region’s ecological conservation.

(3)

The selected controlling factors used in this study have different explanatory powers on soil erosion among different geomorphological types. The vegetation coverage and land-use types are the strongest variables, while terrain factors (slope and terrain niche index) and factors with macroscopic spatial distribution (the population density and population) have lower explanatory power on soil erosion and exhibit better performances for relatively flat platforms. The application of the Geo-detector model needs to take the specificity of the study area and spatial scale into consideration.

(4)

Interactions of factors can enhance the effect degree of single factors. The most influential synergistic groups were vegetation coverage with other factors, especially land-use types in different geomorphological settings. The combination of slope and rainfall also had a significant effect on soil erosion, demonstrating the importance of increasing forest coverage and constructing engineering features in regions of steep slopes.

(5)

Risk detection results indicated that the management of soil erosion needs to strictly follow the idea of “adjusting to local conditions”. Regions of cultivated land and forest land with steep slopes (>25°) in mid-elevation hills were focal areas of soil erosion prevention in Yan’an City.