Assessment of Soil Erosion Dynamics Using the GIS-Based RUSLE Model: A Case Study of Wangjiagou Watershed from the Three Gorges Reservoir Region, Southwestern China

Abstract

The adjacent agricultural watershed is a vital component of the Three Gorges Reservoir Region (TGRR); however, it is affected by serious soil erosion. Assessing soil erosion dynamics in such watersheds is useful for identifying its causes and tendencies to develop, in turn providing scientific information for soil and water conservation at the regional scale. In the present study, the spatial and temporal patterns of soil erosion of a small agricultural watershed in central TGRR were investigated from 2002 to 2014 using the Revised Universal Soil Loss Equation (RUSLE) model, combined with Geographic Information Systems (GIS). The trends and processes of the overall soil erosion intensity were analyzed using spatial overlay analysis and the Markov transition matrix model, respectively. The spatial distribution of soil erosion rates within this watershed was relatively consistent during the study period. Erosion intensity was moderate, with a mean soil loss of 35.1 t·ha⁻¹·year⁻¹. Precipitation was a dominant factor influencing the intensity of soil erosion. Moreover, most erosion intensities shifted closely to middle grades from 2002 to 2008, and declined from 2008 to 2014, indicating that soil erosion in the Wangjiagou watershed has recently decreased. These results suggest that recently implemented integrated soil management practices were responsible for the recently observed erosion patterns.

1. Introduction

Soil erosion is a global concern [1]. Under the background of the climate change and continuing decline in the ratio of available natural resources to population, practices such as unreasonable tillage methods, overgrazing, and construction expose soils to greater risk of erosion [2]. Accelerated erosion has negative social and ecological impacts as it not only decreases water quality, deposits silt in rivers, destroys habitat and causes flooding, but also causes nutrient loss and land degradation [3,4,5,6]. Appropriate policies to reduce the effects of soil erosion within watersheds rely on the analysis of spatiotemporal erosion patterns [7,8,9]. Therefore, understanding the spatial and temporal dynamics of soil erosion is of great significance for local land management and erosion control [10].

Soil erosion is influenced by diverse factors, such as the climate, soil texture, topography, land cover, and human activity. While each variable can be measured independently, it is difficult to predict the overall impact of these interconnected factors [8]. Thus, qualitative, semi-quantitative, and quantitative approaches have been used to simulate erosion processes on watershed and regional scales (including empirical, conceptual and physical-process-based erosion models) [11,12]. Among these models, the Revised Universal Soil Loss Equation (RUSLE) is widely considered the best empirical model for quantitatively evaluating and monitoring soil erosion [13]. On the one hand, RUSLE can be applied to simply predicting soil erosion rates using the watershed features and local hydro-climatic conditions [14]. On the other hand, when combined with Geographic Information Systems (GIS) RUSLE can predict the spatial heterogeneity of soil erosion [15], which is of prime importance for soil and water conservation planning. As a simple, efficient, and robust soil erosion model, RUSLE has also been deemed fit for assessing soil erosion for ungauged watersheds in the Three Gorges Reservoir Region (TGRR) [16]; however, it also suffers from several shortcomings such as spatial scale effects, limited extrapolation capabilities, and challenges dealing with the complexity of the overall erosion process [17].

The TGRR in China refers to the 21 counties and cities in Hubei Province and Chongqing Municipality with an aggregate area of 58,800 km² [18]. The agricultural watershed around the TGRR has a significant disturbance from human activities, which have increased the soil erosion rates [19]. The large proportion (95.7%) of mountainous and hilly regions of TGRR has necessitated widely-distributed slope croplands (78% of total cropland), which are an important source of soil erosion at the regional scale [20,21]. The primary soil type in this region is the purple soil, which is typically eroded during rainfall-runoff events [22]. Thus, the TGRR has been characterized by serious soil erosion in the upper parts of the Yangtze River [23]. Annual soil loss in the TGRR amounts to 157 million t, and approximately 25% of that is delivered to the Yangtze River [18]. Moreover, nutrients and pollutants are transported to surface waters via erosion and runoff, potentially causing eutrophication and trace metal pollution [19,24,25,26,27]. Given this serious erosion of the TGRR, erosion control would greatly improve the sustainability of agriculture and improve water quality.

A large body of literature has been devoted to soil and water conservation and has proposed methods to understand the soil erosion processes and assess the soil quality to minimize the environmental impacts of land degradation [9,18,28,29,30,31,32,33]. Specifically, RUSLE combined with the GIS application model has been used to both quantify soil loss and analyze its spatial distributions for different slopes and land use patterns in the TGRR through a watershed-based approach [16]. The spatiotemporal dynamics of soil erosion in the TGRR have also been evaluated using GIS techniques in 1995 and 2004 [34] and concluded that soil erosion was under control during those years. Nevertheless, these studies either only provided a snapshot of the current situation, or just assessed the dynamics of soil erosion in the past decades. Therefore, a current assessment of soil erosion dynamics within the TGRR is needed to provide updated information to improve sustainable agriculture and water quality.

The current study was undertaken in a small agricultural watershed (Wangjiagou watershed) in central TGRR. The RUSLE model with GIS and Remote Sensing (RS) was utilized to investigate the spatial and temporal patterns of soil erosion risk from 2002 to 2014. The Markov transition matrix model was also applied to analyze the trends and processes underlying changes in soil erosion intensity [35,36,37,38]. The results of this study will provide scientific data to assist decision makers, thereby facilitating the application of beneficial land management practices in this region.

2. Methodology

2.1. Case-Study Area

The Wangjiagou agricultural watershed (29°54′ N, 107°30′ E in the World Geodetic System 1984 (WGS-84)), lies in the middle of the TGRR and covers an area of 67.1 ha (Figure 1) [39]. This area is characterized by the typical hilly topography with elevations ranging from 186 m up to 307 m above sea level. The terrain of the watershed declines from Northeast to Southwest and extends to the Yangtze River. The inclination varies from 0° to 59° with a mean of 15° (Figure 2a). The whole area is surrounded by mountain ridges and thus can be considered as an “enclosed” catchment with only one water flow outlet existing at the shore adjacent to the Yangtze River [39]. This watershed is located in a subtropical zone with a monsoonal climate [39]. The annual average temperature is 22.1 °C, with an average summer high of 39 °C in July and an average winter low of 3 °C in January. Similarly, annual precipitation is 920 mm with 70% rainfall occurring from May to July [39]. The major soil type is purple soil, classified as Calcaric Purpli-Udic Cambisols according to the World Reference Base for Soil Resources (WRB), which developed from sandstones deposits of the Jurassic period [19] with a relatively thin soil layer and short hydraulic retention time. The soil pH ranges from 5.6 to 8.5, and the organic matter content in the soil ranges from 4.4 to 28 g kg⁻¹, with a mean of 5.9 ± 0.36 g kg⁻¹. There are no industrial and mining facilities within the watershed [39]. The primary crops are tuber mustard (var. tumida Tsenetlee), rice (Oryza sativa), and corn (Zea mays) with their seeding and harvest period of early October and late February for mustard and April and middle August for corn and rice [39]. Along with the implementation of the Integrated Small Watershed Management (ISWM) project around the TGRR [18], for soil erosion reduction and local development, several land management practices and control techniques (i.e., level terraces and retention ponds) have also been developed from 2008 to 2009 in the Wangjiagou watershed [40]. There are several land use patterns including paddies, forests, drylands, meadow, and other types (Figure 2b). Owing to the natural resources and full land use patterns, the Wangjiagou watershed can be regarded as characteristic of the surrounding region.

2.2. Materials

In this study, 13-year precipitation data were collected from Fuling District Bureau of Meteorology, Chongqing, China. Sources of data were categorized and listed in

Table 1. The dataset lists.

Dataset	Content	Format	Data Source
Precipitation	Monthly precipitation data in the target area from 2002 to 2014	Doc	Fuling district bureau of meteorology, Chongqing, China
Soil	Subsoil sand fraction, silt fraction and clay fraction of brown-purple soil	Dataset	http://vdb3.soil.csdb.cn/ (China soil database)
Normalized Difference Vegetation Index (NDVI)	Monthly MODIS/Aqua dataset from 2002 to 2014	HDF	http://glovis.usgs.gov/
Digital Elevation Model (DEM)	Dataset with 2.5 m spatial resolution	Grid	College of resources and environment of Southwest University, China
Land use	Dataset with 2.5 m spatial resolution including 11 kinds of land use patterns in the target area	Grid	College of resources and environment of Southwest University, China
Basic geographical information	Watershed boundary at 1:1000 scale	Vector	College of resources and environment of Southwest University, China

. All of the raster data were used in a spatial resolution of 2.5 m in the ArcGIS 10.2 software. Note that MODIS Normalized Difference Vegetation Index (NDVI) dataset (United States Geological Survey, Reston, Virginia, VA, USA) needs to be converted to the same resolution of 2.5 m for further spatial analyses (Detailed information in Section 3.1.4). Moreover, the land use data in this study were derived through digitizing the 1:1000 land use patterns mapped based on a field survey [41].

3. Methods

3.1. RUSLE Model and Its Processing

RUSLE, one of the most widely-used empirical models for predicting soil erosion, is described by the following Equation (1), and it gives an explicit perspective form to understand the relationship between rainfall and soil erosion.

$$A=R\times K\times LS\times C\times P$$

(1)

where $A$ is the computed spatial average soil loss per year (t·ha⁻¹·year⁻¹); $R$ is the rainfall-runoff erosivity factor (MJ·mm·ha⁻¹·h⁻¹·year⁻¹); $K$ is the soil erodibility factor (t·ha·h·MJ⁻¹·mm⁻¹·ha⁻¹); $LS$ is the slope length-steepness factor (dimensionless); $C$ is the vegetation and management factor (dimensionless); and $P$ is the erosion control practice factor (dimensionless).

In this study, the input parameters of RUSLE were processed under ArcMap GIS 10.2 software environments. It is noteworthy that RS technologies were utilized to derive the NDVI data from the Global MYD13Q1 images under an ENVI 5.2 software environment, which are samples of the MODIS/Aqua Vegetation Indices 16-Day L3 Global 250 m SIN Grid product. Moreover, according to the Chinese professional standard “Classification criteria for soil-erosion intensities” (SL190-2007) [42], annual soil erosion intensity will be commonly classified as 6 grades including very slight erosion (<5 t·ha⁻¹), slight erosion (5–25 t·ha⁻¹), moderate erosion (25–50 t·ha⁻¹), severe erosion (50–80 t·ha⁻¹), very severe erosion (80–150 t·ha⁻¹), and destructive erosion (>150 t·ha⁻¹).

3.1.1. Rainfall-Runoff Erosivity Factor (R)

The $R$ factor reflects the potential erosivity associated with the rain for soil erosion, and it is determined by the volume, intensity and duration of rainfall [43]. The $R$ factor is traditionally obtained based on the EI₃₀ index and its corresponding individual rainfall datasets, and the EI₃₀ index refers to the product of rain kinetic energy (E) times the biggest rain density in 30 min (I₃₀). In order to overcome the data limitations, simpler methods based on the monthly or annual precipitation data have been proposed to estimate the rainfall erosivity [43,44]. In our study, the annual $R$ was calculated as Equations (2) and (3) based on the monthly precipitation data [45].

$$R_y=\alpha \times F^{\beta }$$

(2)

$$F={\displaystyle \sum }_{\mathrm{i}=1}^{12}\frac{P_i^2}{P_y}$$

(3)

where $R_y$ is the annual rainfall erosivity (MJ·mm·ha⁻¹·h⁻¹·year⁻¹); $P_i$ is the monthly rainfall amount (mm); $P_y$ is the annual rainfall amount (mm); $F$ is the modified Fournier index; $\alpha$ and $\beta$ are indeterminate coefficients in the model. In our study, $\alpha$ and $\beta$ were respectively given values of 5.249 and 1.205 based on the previous research in the Fulling District [46]. In the estimation of rainfall-runoff erosivity factor ($R$), the spatial distribution of rainfall erosivity was assumed to be uniform due to the relatively small size of the study area (67.1 ha).

3.1.2. Soil Erodibility Factor (K)

Soil erodibility is a major index to evaluate the ability of soil suffering from erosion by water. $K$ is closely associated with the soil texture, structural integrity, cohesiveness, drainage potential, and organic matter. In this study, the $K$ factor is computed by applying Equations (4) and (5) [47].

$$K=7.594\left\{0.0017+0.0494\exp \left[-\frac{1}{2}{\left[\frac{\log (D_g)+1.675}{0.6986}\right]}^2\right]\right\}$$

(4)

$$D_g=\exp \left[0.01{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{f}_i\ln m_i\right]$$

(5)

where, $D_g$ is the geometric mean particle diameter [48]; $f_i$ is the primary particle (sand, silt, and clay) size fraction in percent according to the USDA soil particle-size classification system, while the obtained particle proportions were classified by China soil particle-size system. Thus, the cubic spline interpolation method was applied to uniform these two different classification systems [49]; $m_i$ is the arithmetic mean of the particle size limits of that size.

This $K$ factor algorithm was modified when it was used in China according to the Chinese local soil texture, and the revised $K$ is as following Equation (6) [50] which was also employed in TGRR [51].

$$K_m=-0.00911+0.55066\times 0.1317K$$

(6)

where, $K_m$ is the revised soil erodibility factor (t·ha·h·MJ⁻¹·mm⁻¹·ha⁻¹). In this study, only one estimated $K$ value was applied for the entire watershed due to the single soil type (purple soil) in the study area.

3.1.3. Slope Length and Steepness Factor (LS)

The terrain is a fundamental natural geographical element influencing the soil erosion by specifically affecting the formation and development of soil and vegetation, and deciding the motion state and direction of the surface runoff. The effect of topography on soil erosion in RUSLE model is performed by the $LS$ factor, which is the product of slope length ($L$) and steepness ($S$). In this study, $L$ is estimated as Equations (7) and (8) [52]. $S$ is calculated by taking step coupling methods as Equation (9) [17,53,54].

$$L={\left(\raisebox{1ex}{$\lambda $}\!\left/ \!\raisebox{-1ex}{$22.13$}\right.\right)}^m$$

(7)

$$m=\{\begin{array}{cc}0.2& \theta \le 1°\\ 0.3& 1°<\theta \le 3°\\ 0.4& 3°<\theta \le 5°\\ 0.5& \theta >5°\end{array}$$

(8)

$$S=\{\begin{array}{cc}10.8sin\theta +0.03& \theta <5°\\ 16.8sin\theta -0.50& 5°\le \theta <10°\\ 21.9sin\theta -0.96& \theta \ge 10°\end{array}$$

(9)

where $\theta$ is the slope of the Digital Elevation Model (DEM) (°); $\lambda$ is the horizontal slope length; $m$ is the index related to slope.

3.1.4. Vegetation and Management Factor (C)

$C$ is defined as the ratio of soil loss from land with specific cover and management to the corresponding loss from continuous bare fallow soil. It has been one of the most difficult RUSLE coefficients for estimating over broad geographic areas [55]. The $C$ factor may range from 0 to 1 with specific values in the study region mainly depending on the ground cover [56].

In previous publications, the $C$ factor was assigned directly from literature or field data into a classified crop profile [57], or estimated from NDVI data. In our study, $C$ was computed according to the Equation (10) based on the NDVI [58].

$$C=1-\frac{\mathrm{NDVI}-{\mathrm{NDVI}}_{min}}{{\mathrm{NDVI}}_{max}-{\mathrm{NDVI}}_{min}}$$

(10)

where ${\mathrm{NDVI}}_{max}$ and ${\mathrm{NDVI}}_{min}$ represent the maximum and minimum NDVI values within the study area, respectively, which a certain remote sensor can detect from a fully vegetated surface and bare soil in its resolution. Altogether, the C factor was derived from ENVI 5.2 and ArcMap GIS 10.2 software environments based on MODIS NDVI 16-day images (HDF files). It should be noted that the MODIS Conversion Toolkit (MCTK) is an HDF file conversion and projection utility extension for all known MODIS data products which helps to extract and project the data we need under ENVI 5.2 environment.

3.1.5. Erosion Control Practices Factor (P)

The erosion control practices factor ($P$) reflects the impact of individual support practices that are utilized in combination to reduce soil loss. Such practices commonly include terracing, contour cultivation, and strips among others. In general, uniform $P$ values derived from certain land use patterns in one study site were always applied in different researches regardless of irreplaceable soil condition for each studied area. Thus, $P$ values derived based on land use patterns from the previous study [17,25] were synthesized and a weighted area-mean-value of $P$ was used in this study. The derived $P$ values are represented in Table S1.

3.2. Markov Transition Matrix Model

In the current study, soil erosion transformation maps were firstly obtained to reveal the spatial dynamics of soil erosion intensities based on the theorem of the map algebra method following Equation (11) performed on the ArcMap 10.2 [37].

$${\mathrm{M}}_{\mathrm{ij}}={\mathrm{A}}_{\mathrm{ij}}^{\mathrm{k}}\times 10+{\mathrm{A}}_{\mathrm{ij}}^{\mathrm{k}+1}(\mathrm{i}\le 10\mathrm{and}\mathrm{j}\le 10)$$

(11)

where ${\mathrm{M}}_{\mathrm{ij}}$ expresses soil erosion transformation map between different years, ${\mathrm{A}}_{\mathrm{ij}}^{\mathrm{k}}$ and ${\mathrm{A}}_{\mathrm{ij}}^{\mathrm{k}+1}$ express individual soil erosion map.

The Markov process model is a special type of stochastic movement processes with discrete states and times [59]. It demonstrates that the state of one motor system at time T + 1 is linked to the state at time T rather than its previous states [35,36,60]. Given as a two-dimensional matrix, this model is feasible to assess dynamic changes in soil erosion [61]. In the current study, soil erosion intensities were classified into six grades including very slight erosion (1), slight erosion (2), moderate erosion (3), severe erosion (4), very severe erosion (5), and destructive erosion (6). Thus, firstly, a 6 × 6 Markov transformation matrix was obtained following Equation (12) based on the statistical data derived from the soil erosion transformation maps, then the transition probability matrix model for soil erosion dynamics was established following Equation (13), to analyze the transition trends and processes of soil erosion between different time intervals.

$$P_{ij}=\left[\begin{array}{ccccc}{P}_{11}& P_{12}& P_{13}& \cdots & P_{1j}\\ P_{21}& P_{22}& P_{23}& \cdots & P_{2j}\\ P_{31}& P_{32}& P_{33}& \cdots & P_{3j}\\ \cdots & \cdots & \cdots & \cdots & \cdots \\ P_{i1}& P_{i2}& P_{i3}& \cdots & P_{ij}\end{array}\right]$$

(12)

$$D_{ij}=\frac{P_{ij}}{P_i}$$

(13)

where, the parameter $D_{ij}$ here is the transition probability that represents the $i$ soil erodibility in year T converts to the $j$ soil erodibility in the next year T + 1 [10]. $P_{ij}$ is the transition area from $i$ soil erodibility in year T to the $j$ soil erodibility in the next year T + 1. $P_i$ is the gross transition area converted from $i$ soil erodibility in year T.

4. Results and Discussion

4.1. RUSLE Factors

Annual rainfall erosivity $R$ values (Table S2) showed a significant temporal variation ranging from 1320.41 to 3829.42 MJ·mm·ha⁻¹·h⁻¹·year⁻¹ from 2002 to 2014 (average of 2154.81 ± 741.52 MJ·mm·ha⁻¹·h⁻¹·year⁻¹). The relationship between annual rainfall amount and rainfall erosivity was illustrated by Figure 3. Higher $R$ values of 3497.22 and 3829.42 MJ·mm·ha⁻¹·h⁻¹·year⁻¹ were observed in 2003 and 2013. Even though the total annual precipitation amounts in these two years were comparable with those in other years, the occurrence of most rainfall (about 71% and 63% in 2003 and 2013) was concentrated in three months (from May to July), suggesting that rainfall intensity is probably a primary factor influencing R.

In this study, only one estimated $K$ value (0.02 t·ha·h·MJ⁻¹·mm⁻¹·ha⁻¹) was applied (Table S3), which is comparable with the previous study where $K$ for purple soil in the TGRR was calculated to be 0.0184 t·ha·h·MJ⁻¹·mm⁻¹·ha⁻¹ [51]. As regarding the purple soil, it is weathered from purplish rocks and defined as a soil without any profile development other than A and C horizons and has a relatively high erodibility due to its loose crumb structure. Combined with erosive precipitations and steep topography, as well as the poor permeability of the underlying soil layer, high soil erosion rates can thus be induced [22].

Topography factors (slope length factor $L$ and steepness factor $S$) were calculated through the algorithms mentioned above by using DEM as well as examining the interactions between topography and flow accumulation (Figure 4). The results showed that $L$ and $S$ factors ranged from 0 to 3.49 and from 0.03 to 17.79, respectively. Moreover, these two factors showed a similar distribution trend with an increased value appearing in the mountains and high elevation area, and decreased value presenting near the main drainage network. The average topographic $LS$ factor of the Wangjiagou watershed was calculated to be 8.90 ± 7.88.

Figure S1 described the NDVI variations from July to October during the study period. Four months (July, August, September and October) of NDVI data per year in this study were used to estimate the C factor due to local cropping system [39], as NDVI data from these four months could represent the annual coverage characteristics. Most NDVI values in the study period showed a decrease from July to October, which was associated with the corresponding coverage change at that time period. As mentioned above, corn and rice bloomed in July and harvested in August, and mustard seedlings generally transplanted in early October and fully grown harvested in late February, resulting in an entire coverage decline in the watershed from July to October. Moreover, no distinct changes of the annual $C$ values were observed during the study period (Figure S2), suggesting that the vegetation coverage remains consistent during the study period.

As the ubiquitous soil and water conservation practice in the study area, level terrace played an important role in erosion control. In this study, a weighted area-mean-value of 0.324 was utilized in $P$ estimation. It should be noted here that several soil management practices (i.e., level terrace and retention ponds) were conducted in the studied area from September of 2008 to November of 2009. Specifically, these practices have significant effects on erosion and surface runoff in this agricultural watershed [40]. Thus, the improvement of the resolution of these practices in RUSLE would be beneficial to the accuracy of soil erosion simulation.

4.2. Characteristics of Soil Erosion Dynamics

4.2.1. Temporal Dynamics of Soil Erosion in the Wangjiagou Watershed

Figure S3 showed the spatial soil loss maps at the pixel level from 2002 to 2014 in the Wangjiagou agricultural watershed. The annual soil erosion rates had a similar distribution pattern. The severest soil erosion appeared in high elevation area with steep slopes of 15–25° or over 25° and light soil erosion was found in most of the lowland near the drainage network. Soil loss was calculated by the Zonal Statistics function to be from 22.96 to 62.55 t·ha⁻¹·year⁻¹ with an average of 35.08 t·ha⁻¹·year⁻¹ from 2002 to 2014 (

Table 2. The soil loss from 2002 to 2014 in the Wangjiagou agricultural watershed (t·ha⁻¹·year⁻¹).

Year	2002	2003	2004	2005	2006	2007	2008	2009	2010	2011	2012	2013	2014	average
Soil Loss	26.77	62.55	41.09	39.13	24.55	29.64	30.84	36.95	23.79	22.96	33.58	60.16	24.09	35.08

), suggesting that this studied watershed generally suffered from moderate erosion intensity. From

Table 3. The variations of the area ratio of each soil erosion intensity in the Wangjiagou agricultural watershed from 2002 to 2014.

Erodibility (t·ha−1·year−1)	2002	2003	2004	2005	2006	2007	2008	2009	2010	2011	2012	2013	2014
Very slight (<5)	14.77%	6.78%	10.25%	10.23%	14.83%	12.61%	12.34%	10.22%	14.93%	15.30%	11.04%	6.52%	15.55%
Slight (5–25)	46.10%	19.34%	32.10%	31.84%	46.75%	39.56%	40.31%	33.45%	45.96%	49.88%	37.12%	21.55%	48.36%
Moderate (25–50)	23.40%	24.97%	25.90%	29.00%	25.63%	28.65%	27.88%	28.81%	28.61%	24.45%	29.95%	25.91%	24.78%
Severe (50–80)	11.45%	19.10%	18.97%	16.48%	10.87%	15.03%	13.26%	18.32%	9.52%	8.41%	15.04%	21.01%	9.31%
Very severe (80–150)	3.95%	23.63%	10.84%	11.75%	1.90%	4.12%	6.00%	8.83%	0.98%	1.96%	6.13%	18.34%	1.59%
Destructive (>150)	0.32%	6.19%	1.94%	0.70%	0.01%	0.03%	0.21%	0.36%	0.0019%	0.00%	0.72%	6.67%	0.40%

, the inter-annual variability of different soil erosion intensities generally had a similar trend. Specifically, slight erosion areas accounted for the largest proportion, followed by moderate erosion areas, very slight erosion areas, severe erosion areas, while very severe erosion areas and destructive erosion areas occupied the lowest proportion except for 2003 and 2013. Note that the area ratio of intensities above very severe in 2003 and 2013 accounted for 29.82% and 25.01%, respectively. As for other years, the corresponding ratios were significantly lower than those in 2003 and 2013.

Precipitation and vegetation coverage have been reported to be two primary factors influencing the dynamics of soil erosion [62]. Figure S4 showed the annual variations of rainfall-runoff erosivity, C value (for evaluating vegetation coverage) and soil loss in the study area from 2002 to 2014, in which rainfall-runoff erosivity and soil loss followed a similar trend. Soil erosion can thus be deduced to obtain a weaker response to vegetation coverage than the precipitation which may be a predominant contributor to soil erosion in this watershed. This observation was supported by a previous study that soil erosion in some hilly regions in Loess Plateau was predominantly related to precipitation [13] Thus, we reasoned that very severe and destructive erosion grades occurred in 2003 and 2013 were possibly due to the concentrated precipitation in these two years.

4.2.2. Spatial Dynamics of Soil Erosion in the Wangjiagou Watershed

In order to reveal the spatial dynamics of soil erosion intensity, the spatial overlay analysis method was applied to analyze the space transformation of soil erosion intensity in the study area at two time intervals (from 2002 to 2008 and from 2008 to 2014) (Figure 5). Specifically, Figure 5a indicates that the area with unchanged soil erosion intensity level (0) accounting for 58.06% of the total study area from 2002 to 2008, mainly distributed in a low altitude area and the drainage network area. The percentages of the area receiving a grade reduction (−1) and a grade increase (1) are 12.42% (mainly distributed in the high altitude area in the northeast region) and 25.72% (in the central and southeast part of the study area with sharp slopes), respectively. However, only 2.45%, 0.15%, and 1.20% of the study area received a two and three grades increase (2 and 3) and two grades reduction (−2). Note that these areas with a two or three grades increase from 2002 to 2008 were mainly distributed in orchards. It can be explained that orchards were largely developed in those years accompanying with adjustment of agricultural industrial structure in TGRR. Soil erosion in orchards became serious due to the poor soil conditions of newly constructed slope orchards, a thin soil layer, and the large area of the idle land. As a result, it triggered double pressure for the ecosystem and agricultural sustainable development within this region [63]. From 2008 to 2014, the unchanged area (0) and the area with one grade reduction (−1) occupied a fairly large proportion (57.53% and 30.99%) of the whole area (as shown in Figure 5b), and most of the unchanged area is still mainly located in the drainage network area. While the area with −1 was much more than that between 2002 and 2008 and was primarily distributed in the southern part and the central part with sharp slopes. It was further considered that the implementation of integrative soil management practices from 2008 to 2009 was responsible for the large proportion of area with one grade reduction from 2008 to 2014. As for the area with one grade increase mainly distributed in the southern shrubwood, it may be attributed to the intensive human activities.

To better understand the transition process of soil erosion intensity, the transition probability matrices in

Table 4. The transition probability matrix of soil erosion intensity in the Wangjiagou agricultural watershed from 2002 to 2008.

2002	2008
	Very Slight	Slight	Moderate	Severe	Very Severe	Destructive
Very slight	75.73%	24.27%	0.00%	0.00%	0.00%	0.00%
Slight	2.50%	67.63%	26.24%	3.31%	0.32%	0.00%
Moderate	0.00%	21.13%	47.81%	27.12%	3.93%	0.00%
Severe	0.00%	5.26%	35.08%	29.09%	30.56%	0.00%
Very severe	0.00%	0.00%	14.58%	51.36%	29.18%	4.88%
Destructive	0.00%	0.00%	0.00%	7.47%	86.21%	6.32%

and

Table 5. The transition probability matrix of soil erosion intensity in the Wangjiagou agricultural watershed from 2008 to 2014.

2008	2014
	Very Slight	Slight	Moderate	Severe	Very severe	Destructive
Very slight	93.87%	6.13%	0.00%	0.00%	0.00%	0.00%
Slight	9.85%	79.85%	8.95%	1.36%	0.00%	0.00%
Moderate	0.00%	52.50%	37.04%	10.37%	0.08%	0.00%
Severe	0.00%	5.92%	68.28%	19.63%	6.17%	0.00%
Very severe	0.00%	0.00%	29.93%	54.44%	11.41%	4.22%
Destructive	0.00%	0.00%	0.00%	0.00%	30.54%	69.46%

produced by Equations (12) and (13) were employed in the Wangjiagou agricultural watershed. As can be seen from Table 4, the stabilities of different erosion intensities from very slight to destructive in the period were 75.73%, 67.63%, 47.81%, 29.09%, 29.18%, and 6.32%, respectively. Interestingly, 86.21% of destructive intensity shifted to a very severe intensity, and 51.36% and 35.08% of very severe and severe intensities both decreased to the next levels. Meantime, 24.27% and 26.24% of very slight and slight intensities shifted to slight and moderate intensities, respectively. Most erosion intensities (low or high levels) transited closely to middle grades from 2002 to 2008. As regards to the transition process from 2008 to 2014, erosion intensities from slight to destructive shifted 9.85%, 52.50%, 68.28%, 54.44%, and 30.54% to the previous level, respectively. That is, all erosion intensity grades from 2008 to 2014 generally declined, suggesting that soil erosion decreased during that period (Table 5).

4.3. Consistency with Other Researches in TGRR

The finding in this study agrees well with the published studies conducted in three small agricultural watersheds in the TGRR, with soil loss values of 39.67 t·ha⁻¹·year⁻¹, 34.64 t·ha⁻¹·year⁻¹, and 33.93 t·ha⁻¹·year⁻¹ in slope croplands with slopes of 15°, 14°, and 8° to 15°, respectively [16]. Published study also indicated that soil loss of the purple sloping croplands ranged from 34.64 to 94.52 t·ha⁻¹·year⁻¹ [21]. Furthermore, a study using the ¹³⁷Cs technique showed that the erosion rates on the slope can reach up to 49.37 t·ha⁻¹·year⁻¹ [64]. Although research reported that soil erosion rate in a sloping land of TGRR was only 25.09 t·ha⁻¹·year⁻¹ on the slope of 29° by using ¹³⁷Cs technique [65], that may be attributed to the clayey texture and high erosion resistance of yellow soil as compared to the purple soil [65]. These aforesaid works are comparable and are indirect evidence for the validation of the result of the present study.

4.4. Addressing Future Research and Efforts in the Wangjiagou Watershed

Generally, RUSLE only considers rill and inter-rill erosion from hillslopes induced by rainfall-runoff events, and not processes such as gully and stream-channel erosion [66]. Although we reported a comparable soil loss with other researches in TGRR, limitations of this study are difficult to avoid. Moreover, RUSLE cannot be utilized to evaluate the sediment deposition during the soil erosion process. Thus, we highly recommend the WATEM/SEDEM distributed erosion and sediment transport model for further erosion simulation within agricultural watersheds around TGRR [18]. WATEM/SEDEM not only allows the estimation of the rill, inter-rill, and gully erosion, but also can predict sediment delivery to river channels considering the sediment transport capacity equation [67,68], which is of prime importance to erosion control. Whereas it seems that there are still some challenges to utilize WATEM/SEDEM in TGRR. For example, WATEM/SEDEM focuses much on the spatial and to a lesser degree on the temporal features [18]. Moreover, it also requires a modest number of parameters compared with RUSLE, which may bring a great challenge to ungauged agricultural watersheds like the Wangjiagou watershed in this study.

Soil erosion has been previously demonstrated to be one of the most important processes that can deliver trace metals (such as mercury) and nutrients (nitrogen and phosphorus) from terrestrial to aquatic environments [19,25,69], which considerably increases Agricultural Nonpoint Source Pollution (ANPS) and poses serious threat to the surface water by eroded sediments [20]. Notably, land management practices and the control technique for ANPS in the Wangjiagou watershed have been developed from 2008 to 2009, respectively [19,40]. We thus believe that the information about soil erosion dynamics addressed in this study will be useful to make more effective predictions about ANPS within the study area. Additionally, precipitation has been identified as a key driver for soil erosion in this area, therefore, more effective storm management practices should be implemented (i.e., contour hedgerows and perennial alfalfa mulching) [20], although the crop-mulberry management system for erosion control has been previously reported in this watershed [19]. Finally, a large body of literature has reported the relationship between the climate with landuse changes and soil erosion or sediment yields [70,71,72,73,74], whereas this kind of research in agricultural watersheds around TGRR is quite limited. Therefore, more attention should further be paid in these areas.

5. Conclusions

The present study was conducted in a small agricultural watershed in central TGRR to investigate soil erosion dynamics based on the RUSLE model combined with the GIS technique. The results showed that the study area generally suffered from a moderate erosion intensity with a similar distribution pattern of soil erosion rates in the study period. The transition probability model also highlighted the declining tendency of soil erosion intensity grades from 2008 to 2014, which can be explained by the implementation of integrative soil management practices from 2008 to 2009 within this watershed. Additionally, the concentrated precipitation was mainly distributed in three months (May, June, and July), making rainfall a dominant factor responding to soil erosion, which is valuable to formulate and implement further conservation practices for reducing soil erosion within this area.

In principle, RUSLE only considers rill and inter-rill erosion from hillslope induced by rainfall-runoff events, and not processes such as gully and stream-channel erosion. Although indirect validation of the model results has been shown in this paper, the limitations are difficult to avoid. Another concern in this study is the resolution of the estimated RUSLE factors. For instance, the relatively low resolution of NDVI images may exert high uncertainty for the spatial simulation of soil erosion. For the accurate assessment of erosion modelling, additional data such as the collection of sediments from the bounded plots or other field surveys are necessary, although it is complicated and cost and labour intensive.