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

: The Loess Plateau of China suffers from severe erosion, which results in a great variety of economic and ecological problems. For scientiﬁc 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 2 · 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 signiﬁcantly 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 speciﬁcity of the study area, the Geo-detector model performed better in a relatively ﬂat 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 scientiﬁc use of the Geo-detector model.


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

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 windblown dust and brown soil. The region is subject to arid to semi-arid climatic conditions with a mean annual maximum temperature of 22.9℃ 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).

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, midelevation 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 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, midelevation 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).

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).

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: where A is the estimated soil erosion modulus (t/km 2 ·a), R is the rainfall-runoff erodibility factor (MJ·mm·km −2 ·h −1 ·a −1 ), K is the soil erodibility factor (t·km 2 ·h·MJ −1 ·mm −1 ·km −2 ), 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 30 (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]: 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]: 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]: 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).
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 landuse types and the P values of different land-use types are shown in the following tables (Tables 2 and 3) [44].  Due to different data sources and spatial resolutions, all these factors in the RUSLE model were resampled to a 250 m spatial resolution.

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: 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 2 and σ 2 ; 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). Table 4. Types of interaction between two covariates.
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).
The average erosion modulus of Yan'an City decreased from 9246.57 t/km 2 ·a in 2010 to 7319.21 t/km 2 ·a in 2017; the overall decrease was 1927.36 t/km 2 ·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.  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).

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 2 ·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 average erosion modulus of Yan'an City decreased from 9246.57 t/km 2 ·a in 2010 to 7319.21 t/km 2 ·a in 2017; the overall decrease was 1927.36 t/km 2 ·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 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) (Figure 4c). 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).  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).

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 smallrelief 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 landuse types were the dominant factors ( Figure 5). The explanatory power of vegetation coverage gradient increased with elevation, whereas a "turnpoint" was observed in midelevation 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%). 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) (Figure 4c). 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).  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 midrelief 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%). 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

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). 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 smallrelief 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.

Identification of High-Risk Areas of Soil Erosion
Based on the principle of risk detection, stratification of influencing factors by highrisk of soil erosion (soil erosion modulus ≥ 15,000 t/km 2 ·a) can be determined (Table 7, 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 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). Table 7. High-risk stratification of factors in different geomorphological types and the average erosion modulus (t/km 2 ·a). Land 2021, 10, x FOR PEER REVIEW Figure 6. High-risk areas of soil erosion in the study area.

Geomorphological Types Land-Use Types Proport
Small-relief and low-elevation forest land 2.0

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 landuse types were considered the independent variables (Table 9). 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) 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.

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 (Figures 4 and 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.

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 2 ·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 2 , respectively [49]. The average erosion modulus in the study area in 2010 and 2017 was 9246.57 and 7319.21 t/km 2 , 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.

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 2 , 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 2 ·a, respectively, and the overall erosion degree was relatively lower (Table 11). 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.

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.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to privacy concerns.