**Sai Hu 1,2, Long Li 1,2,3, Longqian Chen 1,2,\*, Liang Cheng 1,2,4, Lina Yuan 1,2, Xiaodong Huang <sup>5</sup> and Ting Zhang 1,2**


Received: 8 July 2019; Accepted: 26 August 2019; Published: 29 August 2019

**Abstract:** It is generally acknowledged that soil erosion has become one of the greatest global threats to the human–environment system. Although the Revised Universal Soil Loss Equation (RUSLE) has been widely used for soil erosion estimation, the algorithm for calculating soil erodibility factor (*K*) in this equation remains limited, particularly in the context of China, which features highly diverse soil types. In order to address the problem, a modified algorithm describing the piecewise function of gravel content and relative soil erosion was used for the first time to modify the soil erodibility factor, because it has been proven that gravel content has an important effect on soil erosion. The Chaohu Lake Basin (CLB) in East China was used as an example to assess whether our proposal can improve the accuracy of soil erodibility calculation and soil erosion estimation compared with measured data. Results show that (1) taking gravel content into account helps to improve the calculation of soil erodibility and soil erosion estimation due to its protection to topsoil; (2) the overall soil erosion in the CLB was low (1.78 Mg·ha−1·year<sup>−</sup>1) the majority of which was slight erosion (accounting for 85.6%) and no extremely severe erosion; and (3) inappropriate land use such as steep slope reclamation and excessive vegetation destruction are the main reasons for soil erosion of the CLB. Our study will contribute to decision-makers to develop soil and water conservation policies.

**Keywords:** soil erosion; RUSLE; soil erodibility; gravel content; Chaohu Lake Basin

## **1. Introduction**

Soil erosion has become a major global environmental hazard [1], posing a serious threat to the ecological environment, natural resources, and socio–economic development [2,3]. It often has a negative impact on downstream areas, including lowered water levels in reservoirs [4], threats from floods and mudslides [5], damage to habitats of species [6], and reduced agricultural productivity [7,8]. For example, the Mediterranean vineyards are among the most degraded agricultural ecosystems affected by extreme soil erosion rates [9,10], and this situation has also been observed in citrus plantations [11,12]. As one of the most severe areas of soil erosion in the world, the Loess Plateau in China has caused great damage to the natural environment, and the economic and social development

of the region [13]. It is acknowledged that soil erosion has become a major threat to sustainability due to the immediate damage it causes to the soil system and the acceleration of the land degradation process [12]. Soil erosion is a natural process that is controlled by a variety of environmental factors, such as topography, soil, climate, and vegetation [14]. Human activities, for example, deforestation, agricultural production, and construction, accelerate the rate of soil erosion [2,15]. The United Nations (UN) has highlighted soil and water protection as a key land-use policy issue, which is an effective manner to address the challenges of the UN Sustainable Development Goals [16,17]. In general, estimation of soil erosion provides a basis for soil and water conservation [18].

Quantitative modeling, either physically based or empirically based, has become a widely accepted approach in soil erosion estimation research [19,20]. At present, many empirical models have been developed for estimating soil erosion [1]. Among them, the Revised Universal Soil Loss Equation (RUSLE) has become the most commonly used in different environmental conditions and on varying scales [21–25]. The parameter factors in the RUSLE model include rainfall erosivity, soil erodibility, slope length and steepness, cover fraction, and support practice [26], as soil erosion is the result of a combination of natural and human factors [27].

Topography is one of the factors determining the amount of soil erosion, and slope plays an important role in increasing soil loss [28]. Smith and Wischmeier first found that soil loss was a polynomial function of slope (θ) [29] and later modified this function by creating the slope factor, which is a polynomial function of the sine of slope (sin θ; see Equation (9)) as it could improve the accuracy of soil erosion prediction on steep areas [30]. McCool et al. proposed an algorithm for calculating the slope factor using a piecewise function and included it in the RUSLE model [31] while Chinese researchers later modified the algorithm of slope factor for slopes above 10◦ using measured data, which has been widely used in the context of China [32]. However, most of the algorithms were based on measurements obtained from runoff plots below 15◦, which may result in less accurate soil erosion prediction on steep areas. Therefore, it is important to improve the slope factor calculation, in order to improve the accuracy of soil erosion prediction, particularly in the areas with high slopes [33]. In this study, we obtained the fitting formula in the form of piecewise function that describes the functional relationship between the sine of slope θ and the slope factor, using the measured soil erosion data with the slopes ranging from 10 to 25◦ and above 25◦.

Soil erodibility is also a key factor related to soil erosion estimation. Among many soil erodibility factor algorithms is the widely used Erosion Productivity Impact Calculator (EPIC) proposed by Williams et al. [34]. However, application of this algorithm to multiple erosion-prone areas of China resulted in soil erodibilities that were greater than measured values for all soil types [35]. This suggests that the EPIC algorithm might not be well suited for soil erodibility estimation in the context of China. Gravel on the soil surface or the top layer of soil has a direct or indirect effect on soil erosion [36]. A series of laboratory-based and runoff plot-based experiments show that rock fragment and gravel content are negatively related to soil erodibility [37–40]. Therefore, it is necessary to modify the prediction of soil erodibility. For improved soil erodibility estimation accuracy, gravel content, an important parameter used to modify the calculation of soil erodibility, should be considered in the RUSLE model. For the purpose of improving the accuracy of soil erodibility estimation, Shi [41] proposed a new algorithm to modify soil erodibility and constructed a piecewise function between relative soil erosion (ratio of soil erosions with different gravel content under the same condition) and gravel content. However, this algorithm has not been tested and applied in soil erodibility and soil erosion studies. In this study, this algorithm was used to modify the soil erodibility for the first time, and then we evaluated its accuracy and estimated soil erosion. China has been tackling soil erosion for decades in the Loess Plateau [42]. However, such an eco–environmental issue also occurs in many other areas, such as the Chaohu Lake Basin (CLB). In recent years, the contradiction between economic development, population, resources, and the eco–environment has become increasingly prominent in this region [43]. The decline of the eco–environmental quality in the CLB, such as eutrophication and increased soil erosion, has seriously affected the sustainability of regional development. It is interesting that the mountainous

areas in the eastern part of the CLB have high gravel content, reaching more than 20% in certain places. In this study, we therefore selected the CLB as the study area and used the RUSLE model for soil erosion estimation to contribute to the general goal of water and soil conservation. Specific objectives are as follows:


#### **2. Materials and Methods**

#### *2.1. Study Area*

At the center of the east Chinese province of Anhui (116◦20 –118◦0 E, 29◦01 –33◦16 N), the CLB consists of 11 administrative districts—2 urban districts and 9 counties—covering a geographical area of 2.04×10<sup>4</sup> km2 (Figure 1). In 2017, the CLB had a population of around 11.52 million, with an urbanization rate of 68.5%, and produced a GDP of approximately 8345 billion CNY (Chinese yuan), according to the Statistical Yearbook of Anhui Province. Among the many water systems in the CLB is the Chaohu Lake, which is profiled as one of the five largest freshwater lakes in China, is one of the main drinking water sources in the CLB, and is replenished by surface runoff and rainfall. With an average elevation of 50.35 m, the CLB features highlands in the southwest and lowlands in the northeast. Influenced by geomorphic types and parent materials, the soil types in the CLB are particularly complex. The rainfall is mainly concentrated in summer and autumn, which accounts for more than 60% of the total annual rainfall, according to the rainfall data provided by the National Meteorological Information Center.

**Figure 1.** The study area: (**a**) the location of Anhui province in China; (**b**) the location of the Chaohu Lake Basin (CLB) in Anhui province; (**c**) the digital elevation model (DEM) of the CLB.

#### *2.2. Data*

Remote sensing images, digital elevation model (DEM) dataset, soil data, rainfall data, and vector data of the study area were used to generate the input variables for the RUSLE model (Table 1).



#### *2.3. Methods*

In this study, soil erosion was estimated using the RUSLE model, where the soil erodibility factor was modified by incorporating gravel content. The technical flowchart of this study is presented in Figure 2.

**Figure 2.** Flowchart depicting methodology of the study.

#### 2.3.1. RUSLE

As a widely used soil erosion prediction model, the Revised Universal Soil Loss Equation (RUSLE) can quantify the soil erosion modulus in different scenarios and reflect the relationship between soil erosion and various impact factors [15]. Use of the RUSLE was made in this study based on the following equation [30,44]:

$$A = \mathbb{R} \cdot \mathbb{K} \cdot \text{LS} \cdot \mathbb{C} \cdot P\_{\text{\textquotedblleft}} \tag{1}$$

where *A* is the average annual soil loss (Mg·ha−1·year−1), *R* is the rainfall–runoff erosivity factor (MJ·mm·ha−1·hr−1·year−1), *K* is the soil erodibility factor (Mg·ha·h·ha−1·MJ−1·mm−1), *LS* is the slope length and steepness factor (unitless) accounting the effect of topography on soil erosion, *C* is the cover fraction factor (unitless), and *P* is the support practice factor (unitless). These factors are detailed below.

• Rainfall–Runoff Erosivity Factor (*R*)

Rainfall erosivity is the potential possibility of soil erosion induced by rainfall, which is the most important external driving force and the dynamic indicator of soil separation and transportation [45,46]. Many of the existing classic and straightforward algorithms for estimating this factor are mainly based on annual rainfall, monthly, daily, or hourly rainfall [47–49]. Since using annual and monthly rainfall data results in less accurate estimation and it is often challenging to access hourly rainfall data, daily rainfall data was used in this study to calculate the rainfall erosivity in the CLB through the algorithm presented in the First National Census for Water [50]:

$$\overline{R}\_k = \frac{1}{N} \sum\_{i=1}^{N} \left( \alpha \sum\_{j=1}^{M} P\_{d\_{kj}}^{\emptyset} \right) \tag{2}$$

$$a = 21.293 \beta^{-7.3967},\tag{3}$$

$$
\beta = 0.6243 + \frac{27.346}{\overline{P}\_{d\_0}} \text{,} \tag{4}
$$

$$\overline{P}\_{d\_0} = \frac{1}{N} \sum\_{i=1}^{N} \sum\_{k=1}^{12} \sum\_{j=-1}^{M} P\_{d\_{ikj}} \tag{5}$$

$$
\overline{\mathcal{R}} = \sum\_{k=-1}^{12} \overline{\mathcal{R}}\_{k\prime} \tag{6}
$$

where *Rk* is the rainfall erosivity of the *k*th month (MJ·mm·ha−1·hr−1), *N* is the sequence length of the calculated data; *M* is the frequency of erosive rainfall in the *k*th month of the *i*th year; *Pdikj* is the rainfall of the *j*th erosive rainfall in the *k*th month of the *i*th year (mm), and daily rainfall ≥ 12 mm is defined as erosive rainfall; α and β are model parameters calculated using Equations (3) and (4); *Pd*<sup>0</sup> is the multi-year average of erosive rainfall (mm); and *R* is the average annual rainfall erosivity (MJ·mm·ha−1·hr−1·year<sup>−</sup>1).

• Soil Erodibility Factor (*K*)

The soil erodibility factor *K* reflects the sensitivity of soil to erosion. The *K* value can be estimated by making observations at a large number of test plots; however, such an approach is difficult to apply to a large watershed [51,52] such as the CLB. To solve this issue, we calculated soil erodibility in this study based on the EPIC model [34] using the following equation:

$$\begin{array}{l} K = \left\{ 0.2 + 0.3 \exp\left[ -0.0256 \mathcal{W}\_d \left( 1 - \frac{W\_i}{100} \right) \right] \right\} \times \left( \frac{W\_i}{W\_i + W\_l} \right)^{0.3} \\\ \times \left[ 1 - \frac{0.25 \mathcal{W}\_c}{\mathcal{W}\_c + \exp\left( 3.72 - 2.95 \mathcal{W}\_c \right)} \right] \times \left[ 1 - \frac{0.7 \mathcal{W}\_n}{\mathcal{W}\_n + \exp\left( -5.51 + 22.9 \mathcal{W}\_n \right)} \right] \end{array} \tag{7}$$

where *Wd* is the fraction of sand (ø 2–0.05 mm) in %, *Wi* is the fraction of silt (ø 0.05–0.002 mm) in %, *Wl* is the fraction of clay (ø < 0.002 mm) in %, *Wc* is the soil total organic carbon content in %, and *Wn* <sup>=</sup> <sup>1</sup> <sup>−</sup> *Wd* <sup>100</sup> .

• Slope Length and Steepness Factor (*LS*)

The slope length and steepness factor (*LS*) can accelerate soil erosion, representing the influence of topographic features on soil erosion [53]. A steeper slope and a longer slope length lead to more serious soil erosion [54]. The *L* value can be calculated by the following equation [55]:

$$L = \text{ (cell size/22.13)}\text{''},\tag{8}$$

where *cell size* = grid cell size (30 m in this study); the value of m varies between 0.2 and 0.5 (0.2 for slopes less than 1%, 0.3 for 1–3%, 0.4 for 3–4.5%, and 0.5 for slopes exceeding 4.5%).

As mentioned in the Introduction (Section 1), Wischmeier and Smith modified the relationship between soil erosion and slope by creating the slope factor *S*, which is the quadratic function of the sine of slope θ (Equation (9)) and applied it in the USLE (Universal Soil Loss Equation) model [30]:

$$S = \, 65.4 \sin \theta \, ^2 + 4.56 \sin \theta \, + 0.0654.\tag{9}$$

Then, the slope factor formula was modified using observation data in the RUSLE model [32]:

$$S = \begin{cases} 10.8\sin\theta + 0.03 & \theta \le 5^{\circ} \\ 16.8\sin\theta - 0.50 & 5^{\circ} < \theta \le 10^{\circ} \\ 21.9\sin\theta - 0.96 & \theta > 10^{\circ} \end{cases} \tag{10}$$

This algorithm improves the prediction accuracy of soil loss on slopes above 10◦. However, it would not be appropriate for soil erosion estimation in the context of the CLB because soil erosion mainly occurs in mountainous and hilly areas with high slopes in the CLB [56]. As such, it would be necessary to modify this algorithm in the case of slopes greater than 10◦. By extraction from scientific research articles, reports, and books, we complied a dataset of measured soil erosion for a number of plots (see Table A1, Appendix A). This dataset covers slopes ranging from 10◦ to 45◦, which are highly representative. We used these sample data to establish the functional relationship between the sine of slope θ and the slope factor.

• Cover Fraction Factor (*C*)

The *C* factor characterizes the restriction of surface vegetation cover on soil erosion as vegetation helps to retain soil and water [1]. Although both spectral vegetation indices and spectral mixture analysis modeled vegetation fractions [57–59] can be used to calculate the *C* factor values [60,61], it is easier to extract spectral vegetation indices than vegetation fractions. In this study, the mostly commonly used vegetation index, Normalized Difference Vegetation Index (*NDVI*), which is the ratio of the difference between spectral reflectance in near infrared and red regions [62,63], was used to calculate *C* factor values according to the following equation [64]:

$$\mathcal{C} = \exp\left[-a \cdot \frac{NDVI}{(\beta - NDVI)}\right],\tag{11}$$

where α and β are parameters that determine the shape of the *NDVI*-*C* curve, and the α-value of 2 and β-value of 1 provide reasonable results of *C* values compared with those estimated as suming a linear relationship [65].

• Support Practice Factor (*P*)

The *P* factor is defined as the ratio of soil loss with a specific support practice to the corresponding loss with upslope and downslope tillage [66]. It is a dimensionless factor with a value between 0 and 1 [44]: 0 means that no soil erosion will occur while 1 means that no soil and water conservation measures have been taken or the measures have completely failed. The *P* factor is closely related to land use type and land use change [67]. Using the maximum likelihood classifier, we classified the land use/cover types in the CLB into six categories, namely farmland, forestland, grassland, water body, construction land, and unused land. The high-resolution satellite images in 2017 provided by Google Earth Pro were used to assess classification accuracy [68]. In total, 500 sample points were randomly generated in the classified image in ArcGIS 10.1 and then imported into Google Earth Pro to retrieve the ground-truthing data. By constructing a confusion matrix, we obtained the overall accuracy (90.1%) and Kappa coefficient (0.86) for this classification. These high values indicate that the classification was well performed and that the classification map could be used for further analysis in this study.

In this study, we used the land use classification of the CLB in 2017 (Figure 3) and assigned the *P* factor values from Zha et al. [56] and Xu et al. [4] for each land use type (Table 2).

**Figure 3.** Land use classification map of the CLB in 2017.

**Table 2.** The *P* factor value of different land use types.


#### 2.3.2. Modifying Soil Erodibility (*K*)

In order to modify soil erodibility, gravel content was considered as a key parameter. The piecewise function proposed by Shi [41] was used for the first time to modify soil erodibility, which describes the functional relationship between relative soil erosion and different gravel content ranges. The modification coefficient *M* (relative soil erosion) was determined by the following equation:

$$M = \begin{cases} 0.0781e^{-0.0249R\_m} & R\_m > 20\% \\ 0.294 - 0.0123R\_m & 10\% < R\_m \le 20\% \\ 1 - 0.0829R\_m & R\_m \le 10\% \end{cases} \tag{12}$$

where *M* is the coefficient for modifying the soil erodibility, and *Rm* is the gravel content. The modified soil erodibility can be obtained by the following equation:

$$K\_{\mathbb{T}} = K \times M\_{\prime} \tag{13}$$

where *Kr* is the modified soil erodibility (Mg·ha·h·ha−1·MJ−1·mm−1), and *K* is the soil erodibility calculated by the EPIC model (Mg·ha·h·ha−1·MJ−1·mm<sup>−</sup>1). Using Equations (7), (12), and (13), the soil erodibility *K* and the modified soil erodibility *Kr* can be calculated. In order to assess the accuracy of this modified algorithm, five sets of measured soil erodibility data extracted from Zhang et al. [35] were used to compare with our modified soil erodibility.

#### **3. Results**

#### *3.1. Fitting Equation of Slope Factor*

Using the measured soil erosion dataset, we obtained the fitted equation with high confidence level between slope factor and sine value in the ranges of 10–25◦ and above 25◦, respectively. The regression analysis shown in Figure 4 shows a strong linear relationship between them.

**Figure 4.** Relationship between the sine of slope and the slope factor: (**a**) slopes between 10◦ and 25◦; and (**b**) slopes above 25◦.

As such, based on the algorithm proposed by Liu et al. [33], the estimation of slope factor in the case of slopes higher than 10◦ can be calculated by the following equation:

$$S\_{-}=\begin{cases} 10.8\sin\theta+0.03 & \theta \le 5^{\circ} \\ 16.8\sin\theta-0.50 & 5^{\circ} < \theta \le 10^{\circ} \\ 16.10\sin\theta-0.89 & 10^{\circ} < \theta \le 25^{\circ} \\ 23.82\sin\theta-2.64 & \theta > 25^{\circ} \end{cases} \tag{14}$$

#### *3.2. Accuracy Assessment of Modified Soil Erodibility*

The calculated *K* values and *Kr* values were shown in Figure 5. We found that the range of the original *K* values was smaller than that of the modified soil erodibility *Kr* values. The means of modified soil erodibility *Kr* were remarkably reduced. This suggests that gravel content has an important effect on the calculation of soil erodibility.

It is clear that the measured *K* values for the five soil types were all within the range of *Kr* (Table 3). By comparing the ratios of the calculated *K* and *Kr* values to the measured values, we noticed that the calculated *Kr* values were by far closer to the measured values than the calculated *K* values. It is therefore believed that this modified algorithm can result in better soil erodibility estimation and has the potential to improve the accuracy of soil erosion prediction.

**Figure 5.** The *K* values and *Kr* values for each soil types.



Note: <sup>1</sup> *R*<sup>1</sup> is the ratio of the mean of*K* to the measured value and*R*<sup>2</sup> is the ratio of the mean of*Kr* to the measured value.

With the modified algorithm of soil erodibility factor, the RUSLE model (Equation (1)) can be transformed as follow:

$$A \; = R \cdot K\_r \cdot LS \cdot \mathbb{C} \cdot P. \tag{15}$$

#### *3.3. RUSLE Factors*

The *<sup>R</sup>* values in the CLB varied from 2856.17 to 8985.13 MJ·mm·ha−1·hr−1·year<sup>−</sup>1, with a mean of 4244.71 MJ·mm·ha−1·hr−1·year−<sup>1</sup> (Figure 6a). Spatially, the values decreased from southwest to northeast with the highest and smallest values of *R* observed in the counties of Shucheng and Feixi, respectively.

The *Kr* value obtained by Equation (13) varied from 0 to 0.38 Mg·ha·h·ha−1·MJ−1·mm−<sup>1</sup> (Figure 6b). Accounting for 67.9% of the basin's area, gleyic anthrosols was characterized by soil erodibility values ranging between 0.25 and 0.28 Mg·ha·h·ha−1·MJ−1·mm<sup>−</sup>1. The *Kr* values for ferralic cambisols, endosalic cambisols, lithic leptosols, skeletic regosols, and stagnic cambisols were nearly equal to 0. The soil types with the highest *Kr* value were the humic cambisols (0.36 Mg·ha·h·ha−1·MJ−1·mm<sup>−</sup>1) and calcaric cambisols (0.38 Mg·ha·h·ha−1·MJ−1·mm<sup>−</sup>1), mainly distributed in the Dabie Mountains in the southwest, the Sangong Mountains in the south, Chaohu east of Chaohu Lake, and mountainous areas in Hanshan.

The *LS* value in the CLB changed from 0.03 to 40.57, with a mean of 1.04 (Figure 6c). The areas with *LS* values > 5 were mainly concentrated in the mountainous areas with high elevations and slopes, while the areas with *LS* values < 0.1 were mostly distributed in flat terrain that was dominated by construction land and lakes.

The *C* value in the CLB varied from 0 to 1, with a mean of 0.73 (Figure 6d). When the *C* value = 1, it represents invalid vegetation cover and management measures and a low *NDVI* value, with a high probability of soil erosion; conversely, when the *C* value = 0, it represents vegetation cover and management measures to inhibit soil erosion with a positive effect. The calculation of *C* value in the CLB shows that the vegetation cover was concentrated in the areas with low *C* value.

**Figure 6.** Maps of the Revised Universal Soil Loss Equation (RUSLE) factors: (**a**) *R*, the rainfall–runoff erosivity factor; (**b**) *Kr*, the soil erodibility factor; (**c**) *LS*, the slope length and steepness factor; (**d**) *C*, the cover fraction factor; (**e**) *P*, the support practice factor.

The results of land use/cover classification (Section 2.3.1) show that the farmland area was 1,196,187 km2 (accounting for 58.7% of the basin area), which is the largest land use type in the CLB. The forestland area reached 4490.87 km<sup>2</sup> (22%) and the construction land area was 2157.64 km<sup>2</sup> (10.6%). The areas of water body, grassland, and unused land were relatively small, comprising only 6.9%, 0.1%, and 1.8%, respectively. The *P* value was as signed according to the land use/cover classification results, varying from 0 to 1 (Figure 6e). It is shown that the areas with large *P* factor values were concentrated in the mountainous area and its surrounding area. Areas with small *P* factor values were mainly concentrated in construction land, water body, and unused land, where we as sume that the probability of their soil erosion is low.

#### *3.4. Soil Erosion Estimation*

With the factors calculated above, we used the RUSLE in the form of Equations (1) and (15) to estimate the soil erosion of the CLB in 2017. In order to show the differences between soil erosion estimation results clearly, estimated soil erosion was divided into six grades according to the Standards for Classification and Gradation of Soil Erosion issued by the Ministry of Water Resources of the People's Republic of China (SL190-2007) [69]: slight (<5 Mg·ha−1·year<sup>−</sup>1), light (5–25 Mg·ha−1·year<sup>−</sup>1), moderate (25–50 Mg·ha−1·year<sup>−</sup>1), intense (50–80 Mg·ha−1·year<sup>−</sup>1), extremely intense (80–150 Mg·ha−1·year<sup>−</sup>1), and severe (>150 Mg·ha−1·year<sup>−</sup>1). The grading maps are shown in Figure 7.

**Figure 7.** Soil erosion grading maps of the CLB: (**a**) soil erosion modulus estimated from Equation (1); and (**b**) soil erosion modulus estimated from Equation (15).

As seen in Table 4, the two equations resulted in different soil erosion moduli and our estimation (1.78 Mg·ha−1·year−<sup>1</sup> by Equation (15)) was slightly lower than the original RUSLE (1.92 Mg·ha−1·year−<sup>1</sup> by Equation (1)). Both calculation results showed that Shucheng was the largest contributor (~5 Mg·ha−1·year−1) to soil loss in the CLB. In contrast, smallest soil erosion moduli were observed in Changfeng, Feidong, and Feixi, all lower than 1 Mg·ha−1·year−1. The comparison also revealed that there was almost no change in the average annual soil erosion modulus for Changfeng, Jin'an, Urban Hefei, Feidong, and Feixi while the estimation from Equation (15) was lower than that from Equation (1) for the rest of the districts and counties.

In addition, we compared the measured data of 11 districts in the CLB provided by Anhui Provincial Water Conservancy Station with the results of Equations (1) and (15) (Figure 8). It can be observed that all the measured data were smaller than the result of Equation (1), where the largest difference was 0.59 Mg·ha−1·year−<sup>1</sup> in Shucheng and the smallest difference was 0.11 Mg·ha−1·year−<sup>1</sup> in Changfeng. In general, the result of Equation (15) was closer to the measured data than that of Equation (1). The largest difference between the result of Equation (15) and the measured data was 0.44 Mg·ha−1·year−<sup>1</sup> in Shucheng and the smallest difference was 0.08 Mg·ha−1·year−<sup>1</sup> in Hexian, both smaller than those for Equation (1). Figure 8 also shows that the differences between the results of Equations (1) and (15) with the measured data were the same in Changfeng, Jin'an, Urban Hefei, Feidong, and Feixi. For the entire study area, the differences between the two results and the measured data were 0.26 Mg·ha−1·year−<sup>1</sup> for the result of Equation (1) and 0.12 Mg·ha−1·year−<sup>1</sup> for the result of Equation (15).


**Table 4.** Average soil erosion modulus of the CLB in 2017 (Mg·ha−1·year<sup>−</sup>1).

**Figure 8.** Comparison of the results of Equations (1) and (15) with the measured data.

Table 5 shows the areas and proportions of each soil erosion grade. The area for each grade was quite similar from the two estimations. It is clear that the erosion-affected area (from slight to severe level) in the CLB in 2017 was 19,087.35 km2, accounting for 93.6% of the total area. The slight level was the largest (~85%) among all the levels, followed the light level (~8%). The areas of intense and extremely intense levels were quite small, both composing less than 1%. Despite a large erosion-affected area, there was no severe-level soil erosion in the CLB.


**Table 5.** Area for each soil erosion grade.

Figure 9 shows the soil erosion grade for each district and county of the CLB; the slight level was the largest at over 60% for each part of the CLB. Shucheng had more light-erosion (>30%) and moderate-erosion (>2%) areas than the other 10 districts and counties. While Wuwei, Shucheng, and Chaohu had the largest intense erosion areas, no extremely intense erosion was estimated in Changfeng, Jin'an, Urban Hefei, Feidong, and Feixi. Figure 9 also reveals that the two estimations resulted in quite similar distributions of soil erosion grades in these districts, except for the intense level for Hanshan (Figure 9d) and the extremely intense level for Chaohu (Figure 9e).

**Figure 9.** The distribution of soil erosion grade in the CLB and its 11 districts and counties: (**a**) Slight; (**b**) Light; (**c**) Moderate; (**d**) Intense; and (**e**) Extremely intense.

#### **4. Discussion**

In this study, we estimated the soil erosion in the CLB in 2017 using the original RUSLE model and the RUSLE model with modified soil erodibility. In addition, the distribution of soil erosion grades in the CLB and its 11 administrative districts were investigated. The interpretation of the results and their implications are given in this section.

#### *4.1. Factors of the RUSLE Model*

Soil erosion is a complex process influenced by a variety of natural and human-induced factors [70]. Five factors are considered to estimate soil erosion in the RUSLE model [71,72], namely rainfall erosivity (*R*), soil erodibility (*K*), slope length and steepness (*LS*), cover fraction (*C*), and support practice (*P*). The *R* factor, *K* factor, and *LS* factor contribute to greater soil erosion [4] while the *C* factor and *P* factor play an important role in preventing soil erosion [27]. By comparing the maps of the *R* factor, *K* factor, and *LS* factor (Figure 6) and the estimated soil erosion (Figure 7), we notice that the three factors are highly consistent with soil erosion in spatial distribution. This helps to explain why soil erosion was higher in Wuwei, Shucheng, Lujiang, Chaohu, Hanshan, and Hexian than the other parts of the CLB. Similarly, Figure 6 shows that the *C* factor and *P* factor play an irreplaceable role in the control of soil erosion, especially in the areas with topographic fluctuation. Since it is difficult to change the natural factors such as rainfall, topography, and soil properties, optimizing land use structure and improving vegetation coverage are considered the most effective measures to prevent soil erosion [18].

#### *4.2. Influence of Gravel Content on Soil Erosion Estimation*

Based on the comparison of the results of Equations (1) and (15) with the measured data, we found that all the calculated soil erosion modulus combined with the modified soil erodibility were closer to the measured data in Wuwei, Shucheng, Lujiang, Chaohu, Hanshan, and Hexian, which are characterized by mountainous area with high gravel content. However, there were no differences between the two results with the measured data in the districts with small mountainous areas, such as Changfeng, Jin'an, Urban Hefei, Feidong, and Feixi. Therefore, we consider that the accuracy of soil erodibility has an important impact on soil erosion estimation [73,74]. The *K* values calculated by the method in the RUSLE model were much larger than the measured values. In this study, the gravel content parameter was added into the modified algorithm for soil erodibility, and the estimated *Kr* values were closer to the measured ones. The accuracy of soil erosion estimation in the CLB was accordingly effectively improved using the modified soil erodibility. Therefore, we believe that this might be because the effect of gravel content on soil erosion was not fully considered.

Many previous studies about the effect of rock fragment and gravel content on soil erosion have also reached similar conclusions. Rodrigo-Comino et al. [37] carried out an investigation with 96 rainfall simulation experiments at the pedon scale and found that the soil losses are inversely proportional to rock fragment cover on the soil surface. Cerdà [38] carried out 20 experiments on bare areas of natural soils and the results showed that water and soil losses were reduced by the rock fragments. Poesen et al. [39] have reported the various effects of rock fragments on soil erosion and the key finding shows that rock fragment cover will offer protection to topsoil and have different efficiencies in different nested spatial scales. The results of two laboratory flume experiments carried out by Jomaa et al. [40] revealed that the rock fragments decreased the sediment transport capacity. These studies provide a reliable support for our views.

#### *4.3. Characteristics of Soil Erosion in the CLB*

Overall, the erosion-affected areas of the CLB were mainly distributed along the SW–NE direction. While the slight-level soil erosion was mostly found in the alluvial plains along the middle and lower reaches of the Nanfei River, Hangbu River, and Tianhe River, and the low mountain and hilly areas with high vegetation coverage (Figure 7), the areas with high soil erosion modulus in the CLB concentrated

in the northeast of Dabie Mountain, the north of Sangong Mountain, and the mountainous areas of Chaohu and Hanshan. Particularly, the population density of Longhekou Reservoir area in Shucheng was high and inappropriate land use existed in this area, such as steep slope reclamation and excessive vegetation destruction. The intensive interaction between human and nature has caused reservoir siltation, thus serious soil erosion problems [75].
