**Variations in Soil Erosion Resistance of Gully Head Along a 25-Year Revegetation Age on the Loess Plateau**

#### **Zhuoxin Chen 1 , Mingming Guo 1,2 and Wenlong Wang 1, \***


Received: 10 October 2020; Accepted: 19 November 2020; Published: 24 November 2020

**Abstract:** The effects of vegetation restoration on soil erosion resistance of gully head, along a revegetation age gradient, remain poorly understood. Hence, we collected undisturbed soil samples from a slope farmland and four grasslands with different revegetation ages (3, 10, 18, 25 years) along gully heads. Then, these samples were used to obtain soil detachment rate of gully heads by the hydraulic flume experiment under five unit width flow discharges (2–6 m<sup>3</sup> h). The results revealed that soil properties were significantly ameliorated and root density obviously increased in response to restoration age. Compared with farmland, soil detachment rate of revegetated gully heads decreased 35.5% to 66.5%, and the sensitivity of soil erosion of the gully heads to concentrated flow decreased with revegetation age. The soil detachment rate of gully heads was significantly related to the soil bulk density, soil disintegration rate, capillary porosity, saturated soil hydraulic conductivity, organic matter content and water stable aggregate. The roots of 0–0.5 and 0.5–1.0 mm had the highest benefit in reducing soil loss of gully head. After revegetation, soil erodibility of gully heads decreased 31.0% to 78.6%, and critical shear stress was improved by 1.2 to 4.0 times. The soil erodibility and critical shear stress would reach a stable state after an 18-years revegetation age. These results allow us to better evaluate soil vulnerability of gully heads to concentrated flow erosion and the efficiency of revegetation.

**Keywords:** soil erosion; gully erosion; vegetation restoration; soil erodibility; land use

#### **1. Introduction**

Soil erosion is recognized as a global environmental problem, which severely damages infrastructure, causes land degradation and water pollution, and threatens the safety of human production and life [1–3]. In the past few decades, many scholars have made many efforts to study the process and mechanism of soil erosion, establish many soil erosion prediction models and try to control soil erosion [4–7]. At present, a set of soil erosion control measures system integrating engineering measures, agricultural measures, and biological measures has been formed [8–10], especially vegetation measures play an extremely important role in soil erosion control [11,12].

Previous studies have shown that revegetation can effectively reduce soil erosion. For example, Wang et al. [13] found that soil detachment capacity of abandoned farmland was 1.02 to 2.29 times greater than four restored lands. Li et al. [14] reported that the ratios of the soil detachment capacity of cropland to those of orchard, shrubland, woodland, grassland, and wasteland were 7.14, 12.29, 25.78, 28.45, and 46.43, respectively. The improvement of soil erosion resistance by revegetation is mainly controlled

by the combination of soil properties and root traits [15–18]. In terms of soil properties, many studies have verified that the revegetation significantly affects soil erosion by changing the soil bulk density, organic matter content, and water-stable aggregate [19–21]. Furthermore, the vegetation root zone is the dynamic interface of soil–plant–atmosphere continuum in partitioning rain and irrigation water into evaporation, transpiration, runoff, and deep drainage [22,23], but is also the home of "green water" which is the source of plant nutrition [24]. Especially, the vegetation root systems also play a great role in protecting soil against flow scouring by affecting soil water movement [25,26]. Root-permeated soils exhibited lower erosion rates primarily through increasing the required shear stress before detachment [27]. Moreover, root growth can bind and bond soil particles and aggregate, thus, enhances soil resistance to erosion [28]. Some root parameters, for example root biomass, length density, and surface area density, were used to estimate the effect of root on soil detachment [21,28–31]. De Baets and Poesen [25] found that soil detachment rate reduced exponentially with increasing root biomass. Some studies also showed that soil detachment was related to root architecture and fibrous root was more effective than tap root in reducing soil loss [24,31]. However, the most of previous studies only focus on the impact of revegetation on soil erosion resistance of hillslopes. In the watershed dominated by gully erosion, the gully head is the main source of soil erosion, but the effect of revegetation on soil erosion resistance of gully heads remains unclear. Therefore, there is a strong need to understand the effect of revegetation on soil erosion of gully head by concentrated flow to develop a more reasonable vegetation model.

Notably, in the gully region of the Loess Plateau, about 63% of total runoff is generated from the loess tableland with a gentle slope of 1–5◦ , which can initiate gully headcut erosion and contribute 86.3% of total sediment [32]. The gully headcut erosion by concentrated flow became the main sediment resource. At present, the gully headcut erosion was controlled effectively due to the implementation of a series of control measures (e.g., the "Three Protection Belts" and the "Green for Grain" project), which, to some extent, was attributed to the fact that the revegetation improves the soil resistance of gully heads to concentrated flow [21,33]. Since some ecological restoration projects were conducted, land use has changed dramatically in the Loess Plateau [34]. Hence, the land use has changed, and the natural succession of vegetation was promoted [35]. With progression in natural restoration of grassland, soil physical and chemical properties and vegetation characteristics (e.g., coverage, community structure, species composition and diversity, and root diameter, density, and diameter distribution) varied greatly [36,37]. These changes would result in dynamic variations in soil erosion resistance. However, the response of soil erosion resistance to vegetation succession process mainly focused on the hillslope in the hilly-gully region of Loess Plateau [15,16], and few studies were conducted to explore the response of soil resistance of gully heads by concentrated flow to vegetation succession process.

Therefore, to evaluate the effect of revegetation process on soil erosion resistance of gully heads and optimize revegetation measures for controlling gully headcut erosion in the gully region of the Loess Plateau, we selected four grasslands with different revegetation ages (3, 10, 18, 25 years) along gully heads with the slope farmland as the control. This study aimed to (1) quantify the effect of revegetation age on soil detachment by concentrated flow, (2) clear the relationships between soil detachment rate and soil and root properties, and (3) confirm the dynamic variation in soil erosion resistance of gully head with revegetation age.

#### **2. Material and Method**

#### *2.1. Study Area*

The study was conducted in the Nanxiaohegou watershed in the Xifeng Research Station of Soil and Water Conservation (35◦41′–35◦44′ N, 107◦30′–107◦37′ E). The watershed has an area of 36.3 km<sup>2</sup> and altitudes ranging from 1050 to 1423 m above mean sea level (Figure 1) in the typical gully region of the Loess Plateau. The climate is temperate continental semiarid. The mean temperature is 10 ◦C,

and the frost-free period is 160–180 days. Annual precipitation is approximately 523 mm, which has the characteristic of annual variation and uneven distribution during the year. In the form of short heavy storms, 58.8% of the rainfall occurs from July to September. The soil type is yellow loamy soil. The original vegetation has disappeared due to human activities. Gully headcut erosion is the main resource of sediment yield in the watershed. Since the 1970s, some soil and water conservation projects, for example the "three protection belts" project and the "Green for Grain" project and so on, were implemented to control soil and water loss, and the vegetation cover of the Loess Plateau increased to 59.6% in 2013. Additionally, the land use has undergone tremendous changes [38]. These efforts also effectively stabilized the gully heads and thus contained the gully headcut erosion [33]. At present, the annual soil erosion module is effectively controlled at the level of 2440 t km−<sup>2</sup> a −1 in the study area, and the vegetation communities comprise mainly planted forests and shrubs and native secondary herbaceous plants [21]. projects, for example the "three protection belts" project and the "Green for Grain" project and so on,

−2 −1

**Figure 1.** Location of the study area on the Loess Plateau and the location of sampling sites in the Nanxiaohegou watershed.

#### *2.2. Sampling Sites Selection*

During our investigation of gully heads, some cracks developed near gully heads. Kompani-Zare et al. [39] and Guo et al. [21] stated that soil samples from 0 to 30 cm depth near the gully heads (the distance was less than 5 m) can represent soil properties of the gully heads. Therefore, in consideration of collapsibility and vertical joints development of loess, the sampling plots were established about 1.0 m meters from gully heads to ensure safety. As a result, four natural restoration grasslands with different ages (3, 10, 18, 25 years) were selected (Figure 1). The natural restoration age was confirmed by consulting the village elders and scientists at the scientific experimental station. The slope aspect and gradient, elevation, soil type, and previous farming practices of the selected sites were similar to minimize the effects of these factors on the experimental results. For comparison, one corn-planted farmland site, with a topography similar to that of the grasslands, was selected as a control. The basic information of the five selected sites is listed in Table 1.

**Table 1.** Basic information of the selected five sampling sites.


— *a viridis*

#### *2.3. Sampling and Measurement of Soil and Root*

In this study, seven soil and root property parameters including soil bulk density, capillary porosity, soil disintegration rate, soil water-stable aggregate, soil saturated hydraulic conductivity, organic matter content, and root mass density were measured. Firstly, three repeated sampling plots (5 m × 5 m) were established in each of gully head sampling sites with the litter layer removed, and topsoil samples (0–30 cm) were collected. Then, three cutting rings (200 cm<sup>3</sup> ) were used to randomly collect soil samples in each plot, and a total of nine samples were oven-dried at 105 ◦C for 24 h to determine the soil bulk density of each gully head site. Similarly, the other 9 soil samples were also collected by cutting rings of 200 cm<sup>3</sup> to determine the soil saturated hydraulic conductivity by applying the constant water head test method. Three cutting rings (100 cm<sup>3</sup> ) were used to collect soil samples for the measurement of soil capillary porosity [33]. Three man-made steel cubical boxes (5 cm in length) were used to collect soil samples for measuring soil disintegration rate by using a disintegration box [14,40]. Lastly, the other three samples were randomly collected in each plot to form a mixed sample. A total of 45 mixed samples were obtained and used for laboratory analyses of organic matter content and water-stable aggregate and its stability. These mixed soil samples were air-dried at room temperature, with large roots and organic residues manually removed. Sieves with apertures (0.25, 0.5, 1.0, 2.5, and 5.0 mm) were used to test the water-stable aggregate. The potassium dichromate external heating method was used to measure the soil organic matter content.

#### *2.4. Hydraulic Flume Experiments*

A hydraulic flume experiment was conducted to determine the soil resistance to concentrated flow upstream gully heads (Figure 2). The size of the flume was 2.0 m long and 0.15 m wide similar to the one used by De Baets et al. [28,31], which was enough to make water flow along the slope soil. An opening (0.5 m length and 0.1 m wide) was set at the bottom of the flume, and a metal sample box with the same size was used to collect undisturbed soil samples so that the surface of the soil sample was at the same level of the flume surface. The space between sampling box and flume edge sealed with painter's mastic to prevent boundary effects. According to the study of Guo et al. [40], the flume experiment was carried out under five different unit width flow discharges of 2, 3, 4, 5, 6 m<sup>3</sup> s −1 , and thus, a total of 100 samples (5 sites × 5 flow discharges × 4 replications) were collected to measure soil resistance of gully heads. To simulate real flow generation conditions, the soil should be saturated by using a watering pot before experiment. During the experiment, a portable flow meter instrument (LS300-A) with 1.5% accuracy was used to measure flow velocity which was regarded as the flow velocity scoring soil area. Runoff and sediment samples were collected with sampling tanks, and the sampling time was recorded. The measured flow velocity was modified according to flow regime [41]. Sampled sediment was oven-dried at 105 ◦C for 24 h to determine the soil loss amount (*SLA*, kg). Thus, the soil detachment rate (*D<sup>r</sup>* , kg m−<sup>2</sup> s −1 ) could be calculated as follows: please check font size, please check all reference citation

$$D\_r = \frac{\text{SLA}}{\text{AT}} \tag{1}$$

where *SLA* is the oven-dry mass of every sediment sample (kg), *T* is the experimental period (s) and *A* is the soil sample area (m<sup>2</sup> ). In addition, the relative soil detachment rate (*RDr*) was calculated as the ratio between *D<sup>r</sup>* for the root-permeated soil samples and that for the farmland topsoil samples, tested at the same condition [28].

−1 

−1

– – –

−1

−3 −2

**Figure 2.** Skitch of scouring flume for determining soil erosion resistance of gully heads.

In addition, the flow depth (*h*, m) and shear stress (τ, Pa) were also calculated as follows:

$$h = \frac{q}{vw} \tag{2}$$

−3

$$
\pi = \rho g h \mathcal{S} \tag{3}
$$

( ) where *h* is flow depth (m), *q* is the flow discharge (m<sup>3</sup> s −1 ), *w* is the width of the flume (m), *v* is the mean flow velocity (m s −1 ). ρ is the water mass density (kg m−<sup>3</sup> ), *g* is the gravity constant (m s −2 ), and *S* is the slope steepness (m m−<sup>1</sup> ).

 After each scouring test, a steel cubical box (10 cm in length) was used to take soil sample in the center of soil area of scouring flume, and the sample was soaked in tap water for about one hour to increase the dispersion of soil and then were placed on a 0.25 mm sieve and washed with tap water using low-pressure head. The living roots, plant debris and some pebbles were left on the sieve. Only the living roots were picked out carefully using tweezers one by one [15]. The washed roots were classified into 4 levels (0–0.5, 0.5–1.0, 1.0–2.0, and >2.0 mm) by vernier caliper and then were oven-dried for 24 h at 65 ◦C and weighed to calculate root mass density (RBD, kg m−<sup>3</sup> ).

#### *2.5. Parameter Calculation after the Experiments*

Soil particle is detached when flow shear stress exceeds the critical shear stress [6]. Soil erodibility parameter (*Kr*) and critical shear stress (τ*c*) were estimated for every natural restoration stage as the slope coefficient and intercept on the abscissa axis of the regression line between soil detachment rate and shear stress as described in the WEPP model as follow:

$$D\_r = \mathcal{K}\_r(\boldsymbol{\pi} - \boldsymbol{\pi}\_c) \tag{4}$$

Generally, soil detachment rate can be considered as zero when root reached the infinity. To quantify the relationship between detachment rate and root mass density, the Hill curve was selected to simulate the relationship between them [21,31,42]. The Hill curve is expressed as follows [43]:

$$RD\_r = \frac{KX\_r^{\
u}}{X\_r^{\
u} + b'} \ (K > 0, \ a < 0, \ b > 0) \tag{5}$$

where*RD<sup>r</sup>* is relative soil detachment rate; *X<sup>r</sup>* is root mass density; *K*, *a* and *b* are constants. The parameter *a* determines the shape of the curve, *b* determines the steepness of the curve and *K* is the asymptote of *D<sup>r</sup>* for infinitesimal *X<sup>r</sup>* values. Additionally, the Hill curve can be used to evaluate the ability of roots to increase soil resistance against concentrated flow erosion. According to Li et al. [44], *b*ˆ(1/*a*) is plant

specific and can be used as an index to compare the effectiveness of different plant roots in reducing soil erosion rates: the lower *b*ˆ(1/*a*), the more effective the plant root. When the value of *X<sup>r</sup>* is *b*ˆ(1/*a*), the soil detachment rates is reduced by 50%.

#### *2.6. Statistical Analysis and Plotting*

The analysis of one-way ANOVA followed by multiple comparisons with LSD was applied to assess the differences of soil properties (Soil bulk density, soil capillary porosity, soil disintegration rate, saturated hydraulic conductivity, organic matter content, and water-stable aggregate) and root mass density among the five revegetation ages. All soil and root variables of each revegetation ages were tested whether the data exhibited a normal distribution and variance homogeneity by Shapiro-Wilk test and Levene test, respectively. If the data failed to meet the two conditions, the Kruskal–Wallis test was performed for the above analysis. The interaction effect of flow discharge and revegetation age was detected using a two-way ANOVA. Pearson's correlation analysis was used to determine linkages among soil properties, root mass density, and soil detachment rate. Relationships among soil detachment rate, soil properties, flow shear stress and restoration age were analyzed by the regression method. The data analyses were conducted in SPSS v. 16.0 statistical software (IBM Corp., Armonk, NY, USA). The figure plotting was conducted by Origin v. 2020 (OriginLab Corp., Northampton, MA, USA).

#### **3. Results and Discussion**

#### *3.1. E*ff*ect of Revegetation on Soil and Root Properties of Gully Heads*

Figure 3 illustrates that the six soil properties of gully heads exhibited a significant increase or decrease with revegetation age. Compared with slope farmland, the soil bulk density (SBD) and soil disintegration rate (SDR) of revegetated gully heads significantly decreased by 5.7–18.6% and 28.8–80.5%, respectively (*p* < 0.05, Figure 3a,c), while the soil capillary porosity (SCP), saturated soil hydraulic conductivity (SHC), organic matter content (OMC) and water-stable aggregate (WSA) significantly increased by 3.9–13.8%, 17.4–236.2%, 34.2–221.8%, and 27.7–64.4%, respectively (*p* < 0.05, Figure 3b,d–f). Figure 4 illustrates that the roots of 1–2 mm in slope farmland had the relatively higher root mass density (RMD, 0.20 kg m−<sup>3</sup> ) and accounted for 39% of total RMD. Notably, after revegetation, the RMD of >2.0 mm was significantly greater than those of the other three root diameters, and it can account for 40–61% of total RMD. When revegetation age was greater than 3 years, there was a significant difference in RMD among four root diameter levels (*p* < 0.05). In addition, we found that the RMD of four root diameters (except for >2.0 mm) showed a non-significant increase in the first three-years and then significantly increased.

These results were similar to previous findings regarding the effects of revegetation on soil properties [45–47]. In fact, the improvement of soil properties of gully heads with revegetation age can be attributed to the accumulation of fresh plant residues in surface soil as well as roots and decomposed root residues in subsurface soil [48]. These materials can be directly transformed into soil organic matter and thus provide energy/carbon sources and nutrients for soil microorganisms [49,50], further promoting the development of soil aggregation and enhancing the cohesion of soil particles [51]. Hence, vegetation restoration would induce the formation of macroaggregates and increase the water stability of aggregates [19,52].

– – – –

– –

−

–

–

– –

'

**Figure 3.** Variation in soil properties with revegetation age. Note: Bar means the 95% confidence interval (95% CI). Different lowercase letters indicate significant difference among different revegetation ages (*p* < 0.05). (**a**) Soil bulk density (SBD); (**b**) Soil capillary porosity (SCP); (**c**) Soil disintegration rate (SDR); (**d**) Saturated soil hydraulic conductivity (SHC); (**e**) Organic matter content (OMC); (**f**) Water-stable aggregate (WSA).

**Figure 4.** Changes in root mass density (**a**) and its proportion (**b**) of different root diameters with revegetation age. Note: Bar means the 95% confidence interval (95% CI). SF refers to the slope farmland. NR3, NR10, NR18, NR25 represents the 3, 10, 18, and 25 years of natural restoration time, respectively. Different capital letters for the same restoration age indicate a significant difference among different root diameters (*p* < 0.05), and different lowercase letters for the same root diameter level indicate a significant difference among different revegetation ages (*p* < 0.05). (**a**) Root mass density; (**b**) Mass accumulated percent.

–

#### *3.2. E*ff*ect of Revegetation Age on Soil Detachment of Gully Heads*

As illustrated in Figure 5a, the *D<sup>r</sup>* of gully heads showed a significant decrease during the 25-year revegetation. This result was not agreed with the conclusion of Wang et al. [16] who stated that soil detachment capacity of sloped lands fluctuated with abandonment time, and the soil detachment capacity of the slope farmland was significantly greater than those of the abandoned farmlands. The difference was mainly attributed to the great difference in erosion environment (e.g., plant type, geomorphological feature, climate) significantly affecting the succession process [36,47]. The mean *D<sup>r</sup>* of slope farmland was 1.6 to 3.0 times greater than those of revegetated gully heads, which indicated that the revegetation played a role in enhancing the soil resistance of gully heads to concentrated flow erosion.

**Figure 5.** Change in soil detachment rate of gully heads with restoration age (**a**), and its relationships with flow discharge (**b**). Note: Bar means the 95% confidence interval (95% CI). SF refers to the slope farmland. NR3, NR10, NR18, NR25 represents the 3, 10, 18 and 25 years of revegetation age, respectively. The different lowercase letters indicate a significant difference among different revegetation ages (*p* < 0.05). (**a**) Site code; (**b**) Unit width discharge.

−4 −5 −4 −5 −5 −6 −5 −7 Figure 5b shows the *D<sup>r</sup>* of gully heads of slope farmland and four restored grasslands varied with flow discharge. The optimal relationships between *D<sup>r</sup>* and flow discharge were fitted, which can reflect the response of *D<sup>r</sup>* of gully heads to concentrated flow induced by rainstorms of different recurrence intervals. It was found that the response of *D<sup>r</sup>* of slope farmland to flow discharge could be expressed as a power function (*y* = *m* × *x n* ), and the *n*-value was greater than 1, indicating the soil loss of gully heads increases at an increased speed with increasing flow. However, for the restored gully heads, the optimal relationships between *D<sup>r</sup>* and flow discharge could be described by a series of linear functions (*y* = *p* × *x* + *q*), and the *p*-value decreased with revegetation age, indicating that the sensitivity of *D<sup>r</sup>* of the gully heads to concentrated flow erosion gradually decreased with increasing restoration age. Besides, the interacted effect of revegetation age and unit width discharge significantly affected *D<sup>r</sup>* (*p* < 0.001) (Table 2).


**Table 2.** Summary of two-way ANOVAs tests.

#### *3.3. Response of Soil Detachment to Soil Properties*

Figure 6 showed that *D<sup>r</sup>* was positively correlated with soil bulk density and soil disintegration rate (*p* < 0.01), but negatively correlated with capillary porosity, saturated hydraulic conductivity, organic matter and water-stable aggregate of >0.25 mm (*p* < 0.01). Regression analysis showed that *D<sup>r</sup>* increased with soil bulk density as a power function (Figure 7a), which showed an opposite trend with the Wang et al. [13] and Yu et al. [15]. Lower soil bulk density was caused by greater root physical and soil organisms' activities, and thus a soil with lower bulk density was harder to be detached. Additionally, *D<sup>r</sup>* decreased with capillary porosity as a logarithmic function (Figure 7b), which was caused by physically binding and chemically bonding effect of root improving soil structure and porosity and hence increasing soil resistance to erosion [28,53]. Soil disintegration rate referred to the dispersion speed of soil contacting with water, which is an important factor determining soil resistance to erosion [14]. In this study, the soil disintegration rate decreased with the revegetation time (Figure 3c) and *D<sup>r</sup>* increased linearly with an increase in soil disintegration rate (Figure 7c). This is attributed to root wedging mechanism preventing soil from detaching that roots can bind soil and tie surface soil layer into strong and stable subsurface soil layer [14,54]. Saturated soil hydraulic conductivity is an integrating parameter for several physical characteristics such as bulk density, porosity, and mechanical composition. The conclusion that the *D<sup>r</sup>* decreased with increasing soil hydraulic conductivity by a power function is reasonable (Figure 7d) because this study and previous research findings have also indicated that changes in soil bulk density and porosity influenced soil detachment and also were affected by revegetation (e.g., Neves et al. [55]; Zhang et al. [56]). A negative power function was found between *D<sup>r</sup>* and soil organic matter content (Figure 7e). The accumulation of soil organic matter in soil could promote the formation of aggregate and enhance the cohesion of soil particles [57]. Hence, water-stable aggregate also was an indicator determining soil resistance to flow erosion [19]. The *D<sup>r</sup>* decreased as a linear function of water-stable aggregate of >0.25 mm (Figure 7f). The results were in agreement with the findings of Li et al. [14]. However, in Wang et al. [13,16] studies, no significant relationships were found between *D<sup>r</sup>* and organic matter and water-stable aggregate of >0.25 mm, probably caused by small variations of the two factors in their studies and difference in land use between their studies and this study (Podwojewski et al. [58]).

– – – – **Figure 6.** Correlation matrix among soil detachment rate, soil properties, and root mass density. Note: Dr, SBD, SCP, SDR, SHC, OMC, WAS, RMD < 0.5, RMD 0.5–1.0, RMD 1.0–2.0, and RMD > 2.0 refers to the soil detachment rate, soil bulk density, soil capillary porosity, soil disintegration rate, saturated soil hydraulic conductivity, organic matter content, water-stable aggregate, root mass density of <0.5 mm, root mass density of 0.5–1.0 mm, root mass density of 1.0–2.0 mm, and root mass density of >2.0 mm, respectively.

– –

– –

**Figure 7.** Relationships between soil detachment rate and soil properties. (**a**) Soil bulk density (SBD); (**b**) Soil capillary porosity (SCP); (**c**) Soil disintegration rate (SDR); (**d**) Saturated soil hydraulic conductivity (SHC); (**e**) Organic matter content (OMC); (**f**) Water-stable aggregate (WSA).

#### *3.4. Response of Soil Loss of Gully-Head to Root Traits*

Significant correlation was found between *D<sup>r</sup>* and RMD of 0–0.5 mm, 0.5–1.0 mm, 1.0–2.0 mm and >2.0 mm (*p* < 0.01, Figure 6), of which the RMD of 0–0.5 mm had the highest correlation with *Dr* , indicating that roots of each diameter level had the significant impact on soil erosion of gully heads, especially the fibrous root of 0–0.5 mm. Furthermore, the Hill curve could well simulate the relationships between *D<sup>r</sup>* and RMD of different root diameters with *R* <sup>2</sup> varying from 0.42 to 0.57 (Figure 8). As illustrated in Figure 8, the *D<sup>r</sup>* showed a rapid decrease when RMD of 0–0.5, 0.5–1.0, 1.0–2.0 and >2.0 mm ranged from 0 to 0.25 kg m−<sup>3</sup> , 0 to 0.3 kg m−<sup>3</sup> , 0 to 0.5 kg m−<sup>3</sup> , and 0 to 1.0 kg m−<sup>3</sup> , respectively, implying that soil erosion of gully heads could be controlled once vegetation restoration or root growth in soil. Although the roots were limited in density and flexible in early revegetation stage, whereas, roots can contribute to soil cohesion and additional strength, and be crucial in reduction of soil erosion [28,59]. Additionally, root system can bind soil and tie surface soil layer into strong and stable subsurface soil layer [54]. Well-developed root system had great physical binding and chemical bonding effect that could well bind soil particles and soil aggregates together and enhance soil resistance to erosion [16,28,42].

In addition, judged by fitted efficiency (*R* 2 ), the optimal results were found in RMD of 0–0.5 mm (Figure 8a), indicating that fibrous root of 0–0.5 mm is the optimal root system reducing soil loss of gully heads. However, Li et al. [44] reported that the ability of plant roots to decrease soil erosion mainly depended on the number of fibrous roots <1.0 mm. Shangguan et al. [53] also found a similar result but recommended root surface area density as the root variable. The reason may be that plant species with contrasting root architectures have a different erosion reduction effect [25]. Additionally, Amezketa [60] and Gyssels et al. [61] reported that monocotyledonous plants are superior to dicotyledonous plants and grasses are better than cereals in stabilizing soil aggregates.

According to Li et al. [44], *b*ˆ(1/*a*) can be used as an indicator to compare the effectiveness of different diameter roots in reducing soil erosion. The lower *b*ˆ(1/*a*), the more effective the diameter root. The relatively lower *b*ˆ(1/*a*) values (0.132 and 0.131 kg m−<sup>3</sup> ) were found in the roots of 0–0.5 mm and 0.5–1.0 mm than 1.0–2.0 mm and >2.0 mm (Figure 8), indicating that the 0–0.5 mm and 0.5–1.0 mm are the most effective roots in reducing soil erosion of gully heads. However, De Baets et al. [25] study the effect of the mixed community of four grasses [*Lolium perenne* (variety: tove), L. *perenne* (variety: starlet), *Festuca rubra* (variety: echo), and *F. arundinacea* (variety: starlet)] on SDR and found the *b*ˆ(1/*a*)

value is 0.79 kg m−<sup>3</sup> that is greater than that of our study. The result fully indicated that the different plant communities had the markedly different influences on reducing soil erosion. The result also suggested us reasonably choosing plant species with different root architectures and root diameters for revegetation at gully heads.

–

–

– – –

−3 −3 −3 −3

– – –

–

– – – **Figure 8.** Relationships between relative soil detachment rate and root mass density of different root diameters. (**a**) Root of 0–0.5 mm; (**b**) Root of 0.5–1.0 mm; (**c**) Root of 1.0–2.0 mm; (**d**) Root > 2.0 mm.

#### *3.5. E*ff*ect of Revegetation Age on Soil Erosion Resistance of Gully Head*

– Rill soil erodibility parameter (*Kr*) and critical shear stress (τ*c*) were employed to characterize the soil erosion resistance of gully heads [13,21], and were determined by the WEPP model (Equation (4)). The fitted linear function between *D<sup>r</sup>* and shear stress was illustrated in Figure 9. The slope of the fitted line is equal to the *K<sup>r</sup>* , and the *K<sup>r</sup>* of the restored grasslands were 31% to 78.6% less than that of slope farmland. In addition, we found that the soil erodibility of 3-year restored grassland rapidly declined by 31% compared with the slope farmland, indicating the short-term revegetation can rapidly reduce soil erodibility of gully heads. The *K<sup>r</sup>* of grasslands in this study ranged from 0.0009 to 0.0029 s m−<sup>1</sup> , which were less than those reported by Li et al. [14]. Wang et al. [13] found that averaged *K<sup>r</sup>* of restored lands of abandoned farmland was 0.0024 s m−<sup>1</sup> that was close to those of this study. The difference was mainly caused by differences in land use, plant species and restoration time. The soil samples were taken from different vegetation restoration models (korshinsk peashrub, black locust, Chinese pine and mixed forest of amorpha and Chinese pine) in the study of Wang et al. [13], and the restoration age (37 years) was greater than that of this study (3 to 25 years). Regression analysis found that the *K<sup>r</sup>* decreased with restoration time in an exponential function and showed a slight change when restored age was greater than 18 years (Figure 10).

−1

– – – –

−1

−3

−3 –

**Figure 9.** Relationship between soil detachment rate and shear stress under different revegetation ages. (**a**) SF; (**b**) NR3; (**c**) NR10; (**d**) NR18; (**e**) NR25.

**Figure 10.** Relationships between soil erodibility and critical shear stress and revegetation age.

 In addition, the τ*<sup>c</sup>* increased with restoration age by a power function (Figure 10). However, the result was inconsistent with the finding of Wang et al. [16] that critical shear stress varied with restoration age in a nonlinear pattern, reaching the minimum at the restored age of 18. The difference in the temporal change of critical shear stress between Wang et al. [13] and this study was caused probably by differences in soil properties and vegetation characteristics of the sampling sites. Compared with slope farmland, the τ*<sup>c</sup>* of the grasslands was improved by 1.2 to 4.0 times, while τ*<sup>c</sup>* of restored land had a little decrease when restored time was more than 18 years (Figure 10). The result further indicated that revegetation can effectively improve the soil erosion resistance of gully head to concentrated flow, and the critical shear stress would reach a stable state after a 18-year revegetation.

– –

— —

– –

### **4. Conclusions**

This study was carried out to explore the effect of revegetation age on soil erosion resistance of gully heads in the gully region of the Loess Plateau. The results showed that revegetation significantly improved soil properties and promoted root accumulation of gully heads. The mean *D<sup>r</sup>* of slope farmland was 1.6 to 3.0 times greater than those of revegetated gully heads. The revegetation can effectively weaken the sensitivity of soil erosion of the gully heads to concentrated flow. The *D<sup>r</sup>* of gully heads was positively related to bulk density and disintegration rate and negatively related to soil capillary porosity, saturated soil hydraulic conductivity, organic matter content, and water-stable aggregate. Roots of 0–0.5 mm and 0.5–1.0 mm were the most effective roots in reducing soil erosion of gully head, and the native plant species with rich root of 0.5–1.0 mm and 0–0.5 mm were recommended as the first choice for revegetation to restrain gully headcut erosion. Revegetation can reduce soil erodibility of gully heads by 31% to 78.6% and improve the critical shear stress by 1.2 to 4.0 times. This study allows us to better evaluate soil vulnerability of gully head to concentrated flow erosion during revegetation. Further studies are needed to quantify the effect of the different combinations of vegetation types with different root architecture types on soil erosion resistance of gully heads.

**Author Contributions:** Conceptualization, Z.C. and M.G.; data curation, Z.C.; formal analysis, Z.C. and M.G.; funding acquisition, W.W.; investigation Z.C. and M.G.; methodology, Z.C.; project administration, W.W.; resources, W.W.; software, Z.C.; supervision, W.W.; validation, W.W.; writing—original draft, Z.C.; writing—review and editing, M.G. and W.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the National Natural Science Foundation of China (42077079; 41571275).

**Acknowledgments:** All authors thank anonymous reviewers for insightful comments on the original manuscript. **Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*

## **Estimation of Soil Erosion and Sediment Yield in the Lancang–Mekong River Using the Modified Revised Universal Soil Loss Equation and GIS Techniques**

**Pavisorn Chuenchum <sup>1</sup> , Mengzhen Xu <sup>2</sup> and Wenzhe Tang 1, \***


Received: 9 December 2019; Accepted: 29 December 2019; Published: 31 December 2019

**Abstract:** The Lancang–Mekong River basin, as an important transboundary river in Southeast Asia, is challenged by rapid socio-economic development, especially the construction of hydropower dams. Furthermore, substantial factors, such as terrain, rainfall, soil properties and agricultural activity, affect and are highly susceptible to soil erosion and sediment yield. This study aimed to estimate average annual soil erosion in terms of spatial distribution and sediment deposition by using the revised universal soil loss equation (RUSLE) and GIS techniques. This study also applied remote sensing and available data sources for soil erosion analysis. Annual soil erosion in most parts of the study area range from 700 to 10,000 t/km<sup>2</sup> /y with a mean value of 5350 t/km<sup>2</sup> /y. Approximately 45% of the total area undergoes moderate erosion. Moreover, the assessments of sediment deposition and erosion using the modified RUSLE and the GIS techniques indicate high sediment erosion along the flow direction of the mainstream, from the upper Mekong River to the Mekong Delta. The northern part of the upper Mekong River and the central and southern parts of the lower Mekong River are the most vulnerable to the increase in soil erosion rates, indicating sediment deposition.

**Keywords:** soil erosion; sediment yield; RUSLE; Lancang–Mekong River basin

#### **1. Introduction**

Soil erosion affects and challenges the world's environment and natural resources [1–7], and economic and environmental dimensions with negative impacts can affect soil erosion, further resulting in low agricultural productivity, ecological collapse and high sedimentation [6–10]. Approximately 84% of the degraded lands around the world are associated with the most relevant issues about the environment with water and wind as the main agents of erosion [7,11–13]. The average soil erosion by water is estimated to exceed 2000 t/km<sup>2</sup> /y with this type of erosion mainly occurring on croplands in tropical areas [14,15]. Human activities and climate change can also be triggered at a much higher rate thus simulating erosion [8,16–22]. Soil erosion by human activities is reportedly 10–15 times faster than any natural process [23]. For instance, approximately 80% of agricultural areas around the world face high to extreme erosion, and the amount of generated sediments can worsen the turbidity of rivers and increase further the concentration of pollutants [24–26]. Moreover, soil erosion and sediment yield can affect humans and the environment severely if sediment quantity exceeds the standard measurement value of aquatic organisms.

Soil erosion is the main part of the initial process of sediment delivery to rivers; in this initial process, displaced soil particles are transformed into sediments due to the influence of an agent of erosion. The amount of sediments can decrease the potential storage capacity of reservoirs and the performance of hydraulic structures [10,27–30]. According to Reference [31], approximately 0.5% to 1% of sediment depositions affect the annual loss of storage capacity of reservoirs around the world, indicating that most dams will likely be left with only 50% of their corresponding volumes by the 2050s. Reference [32] affirms that sediments currently occupy 40% of the reservoir storages in Asia, indicating high loss of storage capacity. These circumstances affect the long-term sustainability of water sources for hydropower dams. The supposedly low sediment yield from the trapping of dams may also cause shoreline erosion, bank erosion and loss of riparian vegetation [33–36].

Lancang–Mekong River basin, as an important transboundary river in Southeast Asia, is one of the largest rivers causing high sediment loads in Asian rivers. According to Reference [37], the average annual sediment load and the specific sediment yield in the Lancang–Mekong River basin is approximately 160 Mt/y and 200 t/km<sup>2</sup> /y, respectively. The upper Mekong basin contributes approximately 50% of the amount of sediments in the Lancang–Mekong River basin [37–39]. Moreover, the Lancang–Mekong River basin is beset by soil erosion and sediment problems because of rapid socio-economic development, population growth, land deterioration and deforestation in the last 50 years, and the problem is most especially caused by the development of hydropower dams in the region [38,40–42]. Many areas are easily vulnerable to soil erosion due to the influence of rainfall, runoff and human activities. In the last few years, the Lancang–Mekong River basin has eroded at an average rate of 5000 t/km<sup>2</sup> /y [33] which is a moderate erosion level, and it tends to increase in intensity continuously from climate change and land-use change. Conversely, sediment yield in the river basin is decreasing from 250 t/km<sup>2</sup> /y to 209 t/km<sup>2</sup> /y, because the sediment quantity is trapped by hydropower dams. Historical sediment load (1960–2013) from China to the lower Mekong River indicates clearly that the amount of sediment loads heavily decreased from 84.7 Mt/y to 10.8 Mt/y and 147 Mt/y to 66 Mt/y at Chiang Saen and Pakse stations, respectively [43].

Previous research attempted to study emphatically the sediment issue in the Lancang–Mekong River basin and some parts of the basin as a means to accumulate knowledge and information for policymakers. The study of sediments in this river can be divided into two main groups. The first group of previous research focused on the changes in sediment load from the construction and operation of dams in the upper Mekong Basin. Reference [44] considered the changing sediment load in the lower Mekong basin because of the possible effects of the cascade dams in the Lancang. Reference [45] considered the effect of sediments from the Manwan Dam in both pre-dam and post-dam stages. Reference [45] estimated the sediment load of the lower Mekong River basin by classifying the rating curve of suspended sediment concentrations obtained from adjacent stations. Reference [46] investigated the nature and magnitude of changing sediment load and their trends in the Lancang–Mekong River basin using available sediment data from 1965 to 2003. Reference [34] analysed the suspended sediment flux and the sediment supply in the lower Mekong River basin using high-frequency measurements obtained from specific stations in Vietnam. Most research in the first group reveals that the construction and operation of dams in the upper Mekong River basin affect negatively the sediment load in this river due to the trapped sediments in the reservoirs. Sediment load also appears with constantly decreasing trends. Meanwhile, the second group of previous research focused on the sediment trapping efficiency of dams in the Lancang–Mekong River basin. Reference [28] analysed and predicted the sediment trapping efficiency of reservoirs in the mainstream of Lancang River. Reference [47] developed an estimation technique for the sediment trapping efficiency of existing and planned reservoirs in the Mekong River using Brune's method. Reference [48] estimated sediment yield based on geomorphic characteristics, tectonic history and available sediment data and, subsequently, considered the cumulative sediment trapping of dams. Most research in the second group indicates that the majority of sediment loads are trapped in existing dams in the upper Mekong River basin, and they will be further trapped if planned dams are operated officially in the near future. However, most of the above studies concentrated only on sediment load data and used the trapped sediment load data of dams obtained from observation stations. Conversely, studies on soil erosion in the Lancang–Mekong River basin requiring both field surveys and other techniques are rare.

Some studies on soil erosion apply the universal soil loss equation (USLE) in combination with GIS and remote sensing techniques to analyse the spatial distributions and patterns of soil erosion in the Lancang–Mekong River basin. The method is convenient for soil erosion analysis, because it can estimate long-term soil erosion. References [49,50] estimated soil erosion in the upper Mekong River basin in Yunnan Province using USLE and analysed spatial patterns with environmental factors. Reference [51] assessed the conserved water and soil ecosystems in Yunnan Province using remote sensing techniques. Reference [52] analysed the spatial distribution of soil erosion in north-western Yunnan (Lancang River) based on the revised universal soil loss equation (RULSE) and GIS techniques. Reference [10] estimated the impact of soil erosion on the reservoirs in Yunnan Province using USLE. Reference [53] conducted a soil loss vulnerability analysis of the Mekong River basin by applying USLE. Nonetheless, the above studies identified the limitations of the USLE model, including the development of input data for new areas to satisfy the long-term data requirements, difficulties in assessing gully erosion and large-scale areas, estimation of soil loss only and insufficient computation of sediment deposition. The RUSLE model was developed accordingly to improve the estimation of potential soil erosion. The input factors in RULSE can be used by using values from the literature or adapted for empirical and statistical data in combination with GIS software. In addition, the RUSLE results are valid in terms of estimating the risks of water erosion.

Previous studies mostly investigated the changing sediment load and the sediment trapping caused by dam construction and operation. Nonetheless, the understanding of soil erosion and soil deposition is also highly important. Soil erosion, as the main part of the sediment process, can be used to plan countermeasures for the Lancang–Mekong River basin. Previous studies also emphasized that soil erosion research should focus on the simulation of sediment erosion, but they did not consider sediment deposition. Hence, this research aimed to develop methods to calculate sediment deposition and erosion based on the RUSLE model and GIS techniques and, subsequently, evaluate the impact of soil erosion on hydropower dams in the Lancang–Mekong River basin. This study only considered suspended sediment despite the limitation of the model. In addition, the factors that can influence potential and actual soil erosion in the Lancang–Mekong River basin were also determined. The simulation period of the study covered from 2000 to 2015 depending on the available data in the analysis.

#### **2. Study Area**

The Lancang–Mekong River basin is a transboundary river in Southeast Asia (Figure 1). Originating from China's Qinghai–Tibet Plateau, the source of the river is located in Yushu of Qinghai Province. By its name, Lancang River represents the upper Mekong River basin in China, while the downstream part is located in Yunnan Province. Along with the river portions in Myanmar, Laos PDR, Thailand, Cambodia, and Vietnam, the lower Mekong River basin has a length of 4909 km and a coverage area of 795,000 km<sup>2</sup> [38,54,55]. The average annual water discharge is approximately 475 km<sup>3</sup> [38,44,54,56]. Thus, approximately 24% of the total area comprises the upper Mekong River basin, with contribution rates of 15% to 20% of the water flow to the Lancang–Mekong River basin. Most areas comprise complex mountains and hills and deep valleys [10,38,53]. In addition, approximately 76% of the total area is developed by major tributary systems from the lower Mekong Basin, especially Lao People's Democratic Republic (Laos PDR) [38,53]. The elevation of the basin varies from 0 to 6549 m above sea level. The different elevations have varying distributed agriculture depending on climatic zones and temperature. Moreover, various elevations of the river have water development projects, such as cascade hydropower dams, in both the mainstream and the sub-basins. Furthermore, soil erosion in the Lancang–Mekong River basin results in sediment deposition. The sediments affect the dams in all statuses (i.e., operation, under construction, and planned). The dam system of the Lancang–Mekong River basin comprises 133 dams [38,57,58] including those in the mainstream and the sub-basins.

**Figure 1.** Location and elevation of the study area and location of dams in the mainstream river.

#### **3. Materials and Methods**

#### *3.1. RUSLE*

The RUSLE model, based on the USLE model, was developed by the US Department of Agriculture. The RUSLE is an empirical soil erosion model and has been recognised as a standard method to calculate the risk of average soil erosion on land. The RUSLE is also the most popular model for estimating average soil erosion in water [59], and it is simple to integrate with GIS and remote sensing [10,60–62]. Furthermore, RUSLE can provide international applicability and comparability for the results and methods, as the model can be adapted and applied in many regions globally. The RUSLE model can be expressed as follows:

$$A = \mathbb{R} \times \mathbb{K} \times \mathbb{L} \mathbb{S} \times \mathbb{C} \times \mathbb{P} \tag{1}$$

where:

*A* is the mean annual soil loss (t/ha·y);

*R* is the rainfall erosivity factor (MJ·mm/ha·hr·y);

· · *K* is the soil erodibility factor (t·hr/MJ·mm);

*LS* is the topographic factor (dimensionless);

*C* is the cropping management factor (dimensionless); and

*P* is the support practice factor (dimensionless).

The assessment of soil erosion in the Lancang–Mekong River basin can be classified into five levels according to the Soil Erosion Standard Document–Technological Standard of Soil and Water Conservation (SD238-87) of Reference [63].

#### 3.1.1. Rainfall Erosivity Factor

Rainfall plays an important role in the process of soil erosion and sedimentation and leads to water erosion, such as splash erosion, sheet erosion, rill erosion and gully erosion, caused by water flow. Soil particles, which are transported away from a site by the flow, are those detached by rainfall impact [64]. Therefore, high-potential erosion can be determined by rainfall intensity and storm duration. Normally, the relationship between total storm energy (E) and maximum 30 min intensity (I30) can be regarded as the *R* factor, as reported by Reference [65]. Given the limitation of precipitation data about the river, the *R* factor is derived from the Asian Precipitation Highly Resolved Observational Data Integration Towards Evaluation of the Extreme Events (APHRODITE) for the period from 2000 to 2015 which also correspond to the daily gridded precipitation data for Monsoon Asia [66]. This project is developed from the daily rain gauge data for the Asia region and cover nearly 12,000 stations. This study has selected the highest fine-gridded resolution (spatial resolution of 5 km) of available precipitation data. For the conditions in the Lancang–Mekong River basin, this study chose the formula of the *R* factor from References [41,67] which applied the assessment of the *R* factor in Southern China. Equation (2) is appropriate, because the climate and area conditions in Southern China are almost uniform to those in the Lancang–Mekong River basin.

$$R = \sum\_{i=1}^{12} \left( -1.15527 + 1.792P\_i \right) \tag{2}$$

where *R* is the rainfall erosivity factor (MJ·mm/ha·hr·y), and *P<sup>i</sup>* is the monthly rainfall (mm).

#### 3.1.2. Soil Erodibility Factor

The effect of soil characteristics and soil properties on soil erosion can be represented by the soil erodibility factor (*K*), because this factor shows the physical and chemical properties of the soil through the equations related to soil texture, soil organic matter and percentages of sand, silt, and clay. Furthermore, the *K* factor is based on soil permeability and particle size distribution. The *K* factor is strongly related with the *R* factor through the soil erosion rate per kinematic energy of rainfall erosivity index. The observed data of the local soil properties in the Lancang–Mekong River basin are extremely difficult to access. Thus, the soil data in this study were derived from the SoilGrids map which is developed and maintained by ISRIC–World Soil Information. This study used the available soil data grid with a spatial resolution of 1 km. The data on soil properties were analysed using the methods in References [68,69], in which the percentages of silt, clay, sand and organic carbon fraction were calculated by Equations (3)–(6). Soil erodibility was computed according to the method in Reference [70] as shown in Equation (7). Then, the unit of the *K* factor was transferred to the International System of Units (SI) [70]. This method is widely used for the analysis of the *K* factor for soil properties such as soil structure and particle-size distribution [10,53,68,69,71,72].

$$f\_{csand} = \left\{0.2 + 0.3 \exp\left[-0.256 m\_s \left(1 - \frac{m\_{silt}}{100}\right)\right] \right\} \tag{3}$$

$$f\_{cl-si} = \left(\frac{m\_{silt}}{m\_c + m\_{silt}}\right)^{0.3} \tag{4}$$

$$f\_{\rm org\gets} = \left\{ 1 - \frac{0.25 \text{org\gets C}}{\text{org\gets C} + \exp[3.72 - 2.95 \text{org\gets C}]} \right\} \tag{5}$$

$$f\_{\text{bias}} = \left\{ 1 - \frac{0.7\left(1 - \frac{m\_s}{100}\right)}{\left(1 - \frac{m\_s}{100}\right) + \exp\left[-5.51 + 22.9\left(1 - \frac{m\_s}{100}\right)\right]} \right\} \tag{6}$$

$$\mathcal{K} = f\_{\text{csand}} \times f\_{\text{cl\\_si}} \times f\_{\text{org}} \times f\_{\text{hisand}} \tag{7}$$

where *K* is the soil erodibility factor, *fcsand* is the function of high-coarse sand content in soil, *fcl*−*si* is the function of clay and silt in soil, *forgC* is the function of organic carbon content in soil, *fhisand* is the function of high sand content in soil, *m<sup>s</sup>* is the percentage of sand fraction content (particles with diameters from 0.05 to 2 mm) (%), *msilt* is the percentage of silt fraction content (particles with diameters from 0.002 to 0.05 mm) (%), *m<sup>c</sup>* is the percentage of clay fraction content (particles with diameters of <0.002) (%), and *orgC* is the percentage of organic carbon content of the layer (%).

#### 3.1.3. Topographic Factor

The topographic factor (*LS*) includes slope length (*L*) and slope steepness (*S*), which are the two important influencing parameters of soil erosion. Both GIS and remote sensing techniques were applied to access the *LS* factor in the RUSLE equation using the digital elevation model (DEM) [73]. For a large area, grid resolution is important for soil erosion estimation [74]. Changes in grid size affect steepness values, both directly and indirectly. The *L* factor depends on grid size and steepness, while the *S* factor affects steepness only. Hence, if the DEM data have a high resolution, then the model output can increase the accuracy of the *LS* factor in the RUSLE model [75,76]. Digital elevation model images with a 1 km resolution were downloaded from the US Geological Survey (https://earthexplorer.usgs.gov). Past researchers applied high-resolution DEM images for soil erosion determination because of these images' good accuracy and reliability [6,10,16,20,49–51,53,60,62,73,76–79]. The calculation of the *LS* factor can be based on the RUSLE principle by using the GIS software as explained in References [20,73,78,80–82]. The *S* factor was calculated in two conditions (Equations (8) and (9)), and the *L* factor was computed with Equation (10). Then, the *LS* factor in each grid cell was coupled in Equation (11).

$$S\_{\text{factor}} = 10.8 \text{sin}\Theta + 0.03; \text{slope gradients} < 9\% \tag{8}$$

$$S\_{\text{factor}} = 16.8 \sin \theta + 0.50; \text{slope gradients} \ge 9\% \tag{9}$$

$$L\_{factor} = \left(\frac{\lambda}{22.12}\right) \times \left(\frac{\frac{\left(\frac{\sin\theta}{0.0896}\right)}{\left(3\sin\theta \times 0.8 + 0.56\right)}}{1 + \frac{\left(\frac{\sin\theta}{0.0896}\right)}{\left(3\sin\theta \times 0.8 + 0.56\right)}}\right) \tag{10}$$

$$\text{LS} = \text{L}\_{\text{factor}} \times \text{S}\_{\text{factor}} \tag{11}$$

where λ is the length of the slope, *Lfactor* denotes the slope length factors, and *Sfactor* is the slope steepness factor.

#### 3.1.4. Cropping Management Factor

Vegetation cover is one of the most important factors affecting the erosion process and the development of rivers [64]. Moreover, vegetation cover can shield the soil surface from the impact of falling rain and slow down the velocity and scouring power of runoffs. Normally, vegetation cover can be depicted by the cropping and management practices in an area through the *C* factor. The range of the *C* factor is between 1 and 0. If the *C* factor is equal to 1, then no vegetation cover (i.e., bare land) exists in that area. If the *C* factor is close to 0, then strong vegetation cover exists, indicating protection against soil erosion.

The product of remote sensing data from the Moderate Resolution Imaging Spectroradiometer (MODIS), with a cell size of 250 m in spatial resolution, was applied. The MODIS is a good choice for large-area coverages. The normalized difference vegetation index (*NDVI*) was used in this study to estimate the *C* factor following the method of [83]. The detailed equations were given by Equations (12) and (13). The MODIS' remote sensing can investigate all months, from the historical period to the present (2000–2015), to investigate the study area.

$$\mathcal{C} = \frac{(-NDVI + 1)}{2} \tag{12}$$

$$NVDI = \frac{(NIR - RED)}{(RED + NIR)}\tag{13}$$

where *C* is the cropping management factor, *NDVI* is the average of the normalized difference vegetation index, *NIR* is surface spectral reflectance in the near-infrared band, and *RED* is the surface spectral reflectance in the red band. Both *NIR* and *RED* were extracted from the MODIS images. In reflecting the vegetation cover and the agricultural activities in the Lancang–Mekong River basin, the five months of January, April, July, October and December [38] were selected from 2000 to 2015. The average *NDVI* was calculated from these data covering 16 years.

#### 3.1.5. Support Practice Factor

The support practice factor was used to express the effect of land use and land cover on soil erosion. The *P* factor describes the change in potential erosion by flowing water through the effect of supporting conservation practices such as contouring, buffer strips and terraced contour farming [6,53,65,77,84]. The maximum value of the *P* factor is usually set to 1.0 to mean no erosion control solution. A decreasing value of the *P* factor means that flowing water is reduced in terms of both volume and velocity. Moreover, a decreasing *P* also means reduced intensity of sediment deposition on the surface [85]. Given the many limitations, the *P* factor was determined on the basis of the land cover type from the *C* factor (Table 1) as suggested by [86]. Land-use type was obtained from the product of the MODIS' remote sensing with a cell size of 250 m for the spatial resolution.


**Table 1.** Land cover classification and the *C* and *P* factors [86].

#### 3.1.6. Application of GIS Tools

The input data, such as rainfall, types of land use, and land cover, terrain and soil properties, in the RUSLE model were imported and calculated using the functions in ArcGIS 10.5. The five factors were analysed according to the spatial resolution and the coordinate system of their original data. The final results of the quantitative output of soil erosion were generated as the maximum grid with 5 km of spatial resolution depending on the original data. Soil erosion in the Lancang–Mekong River basin was analysed using the results of two types of erosion (i.e., potential soil erosion and actual soil erosion), as shown in Equation (1), in the spatial distribution. The *R*, *K*, *L*, and *S* factors were considered as potential soil erosion, whereas the *R*, *K*, *LS*, *C*, and *P* factors were examined as actual soil erosion.

#### *3.2. Descriptive Statistics in the RUSLE Model*

Soil erosion can be identified in each factor of the RUSLE model, indicating the influence of soil erosion on a specific area [6]. The RUSLE model is transformed into logarithmic form in Equation (15), and multiple linear regression must be applied to examine the relationships among all factors, as shown in Equation (16), and the effects on the soil erosion rate.

$$
\ln(A) = \ln(R \times K \times LS \times \mathbb{C} \times P) \tag{14}
$$

$$
\ln(A) = \ln(R) + \ln(K) + \ln(LS) + \ln(C) + \ln(P) \tag{15}
$$

$$\ln(A) = \beta\_0 + \beta\_l(\ln R) + \beta\_l(\ln K) + \beta\_k(\ln LS) + \beta\_l(\ln C) + \beta\_h(\ln P) \tag{16}$$

where ln(*A*) is the logarithm of soil erosion rate, ln(*R*, *K*, *LS*, *C*, and *P*) denotes the logarithmic value of the input factors in the RUSLE model, β<sup>0</sup> is the intercept of soil erosion rate (constant term), and βi–h is the estimated regression coefficient of each explanatory variable. Different units of the input factors are reflected through the standard coefficient (β) in Equation (16). The factors of multiple linear regression in logarithmic form can be explained as follows: if one of the factors in the RUSLE model increases by 1% in standard deviation, then βi–k percent of the standard deviation leads to an increased value of soil erosion rate (*A*). This study sets the statistical significance level at 95% confidence in SPSS. Nonetheless, given the differences in the spatial resolutions of the input factors, some factors (*K*, *LS*, *C*, and *P*) were estimated as 5 km (*A* and *R* factors) in spatial resolution using the spatially averaged values assigned in the function of ArcGIS.

#### *3.3. Technique of Sediment Yield Estimation*

References [20,87] proposed a new technique to estimate sediment yield or sediment deposition in each sub-basin of Thailand by modifying the original RUSLE model. They regarded the suspended sediment flow from one grid cell to the other grids as dependent on the sediment yield of the original grid cell (*Sy*) and the average sediment yield capacity of sub-basin (*Sc*). If *S<sup>y</sup>* is greater than *Sc*, then the sediment moves to the next site. By contrast, if *S<sup>c</sup>* is more than *Sy*, then the sediment is deposited. *S<sup>y</sup>* is calculated using the individual parameters in each grid cell (Equation (17)). In the same way, *S<sup>c</sup>* is calculated using the original RUSLE model with the area-averaged parameters (Equation (18)). This technique was only developed for the assessment of suspended sediment. It is not appropriate for analysing the total sediment form (i.e., bed load and suspended sediment).

$$S\_{\mathcal{Y}} = f(I\_1, I\_2, \dots, I\_5) \tag{17}$$

$$S\_{\mathcal{C}} = f \begin{pmatrix} \sum\_{i=1}^{n} I\_1 & \sum\_{i=1}^{n} I\_2 & \sum\_{i=1}^{n} I\_5\\ \overline{A\_{ba\sin}} \text{ } \overline{A\_{ba\sin}} \text{ } \dots \text{ } \overline{A\_{ba\sin}} \end{pmatrix} \tag{18}$$

$$D\_i \quad \text{if } S\_y < S\_c \tag{19}$$

$$T\_i \quad \text{if } S\_y > S\_c \tag{20}$$

where *S<sup>y</sup>* is sediment yield, *S<sup>c</sup>* is sediment capacity, *I<sup>i</sup>* represents the parameters in the RUSLE model (*R*, *K*, *LS*, *C*, and *P*), *Abasin* is an area of the sub-basin, *n* is the number of data in each sub-basin, *D<sup>i</sup>* is the sediment deposition in a cell *i*, and *T<sup>i</sup>* is the sediment transportation in cell *i*. *S<sup>y</sup>* is the result of actual soil erosion by computing from the RUSLE input factors. *S<sup>c</sup>* is calculated from the summation of each parameter in the RUSLE model dividing an area of the sub-basin. The five outcomes then are multiplied as *Sc*.

The above technique can show the spatial distribution of sediment yield and sediment deposition in the Lancang–Mekong River basin, indicating an integrated consideration of the sediment issue which is the main problem for water development projects in this river. Furthermore, the technique is extremely useful in studying the influence of dam construction on sediment budget, because the loss of storage capacity of dams and the reduced transport of sediments downstream are caused by sedimentation which, in turn, is the result of soil erosion [32]. Dam design and sediment management

in operations planning can be arranged properly if the sediment budget of the river is primarily determined in dam construction.

#### *3.4. Observed Sediment Data*

The results from net sediment mapping or sediment deposition and erosion mapping are estimated and compared with the observed sediment data from relevant organizations, such as MRC, and the literature for verification [32,44,56,88,89]. The present study collected, from 15 stations, the average sediment load and specific sediment yield (SSY) data for each sub-basin (Figure 2) from the years 1952 to 2011 (60 years) to cover the whole basin (see Supplementary Materials) which is the time period of the data collection. Sediment loads were estimated from the suspended sediment concentration (SSC) and water discharge using the sediment rating curve, and the SSY data in the Lancang–Mekong River basin were estimated based on historical geological and geomorphological characteristics of each sub-basin [48] and historical sediment load. The results of this study will be verified with SSY in each sub-basin only. Each observational station is a representative of a sub-basin in the Lancang–Mekong River basin for verification between observed SSY (1952–2011) and estimated SSY from the modified RUSLE model (2000–2015) (see Table 4).

**Figure 2.** Location of sediment observational stations in the Lancang–Mekong River basin.

### **4. Results**

#### *4.1. Soil Erosion Factors*

#### 4.1.1. Rainfall Erosivity Factor

The values of the *R* factor were analysed using Equation (2). Figure 3a shows the spatial distribution of the *R* factor for the Lancang–Mekong River basin. The range of the *R* factor was 65.6–524.3 MJ·mm/(ha·hr·y) with a mean of 294.9 MJ·mm/(ha·hr·y). The standard deviation was 80.3. The lowest values for the *R* factor were distributed mostly in the upper Mekong River basin or Lancang River in China. Meanwhile, the highest values for the *R* factor were distributed primarily in the sub-basins of Laos PDR and Cambodia and the Mekong Delta, because those areas are located along the direction of monsoon storms from the South China Sea in seasonal cycles. According to the results, the *R* factor increased from the lower basin to the upper basin, a scenario explaining the influence of climate and temperature on the river.

**Figure 3.** (**a**) *R* factor and (**b**) *K* factor.

#### 4.1.2. Soil Erodibility Factor

Major soil groups in the Lancang–Mekong River basin (Figure 3b) were determined using the SoilGrids database of ISRIC–World Soil Information [90,91]. The *K* factor was calculated with Equations (3)–(7). The range of the *K* factor was 0.012–0.0397 t/(hr·MJ·mm), with an average of 0.0258 t/(hr·MJ·mm). The standard deviation was 0.0012. The spatial distribution in Figure 3b indicated that the *K* factor decreased from the upper basin to the lower basin, but some areas of the Mun and Chi River basins in Thailand had high *K* values. In the Lancang–Mekong River basin, the highest elevation areas were identified by the highest *K* values, whereas the lowest elevation areas were identified by the lowest *K* values. This result corresponded with the findings in Reference [10], in which the *K* values correlated with the variation of the terrain; moreover, highly significant *K* values were found for high elevation areas such as mountains. Orthic Acrisols (Ao), Lithosoils (I) and Ferric Acrisols (Af) are the largest areas in the Lancang–Mekong River basin, and they accounted for approximately 59% of the total basin, while the other soil groups accounted for 39%.

#### 4.1.3. Topographic Factor

Topographic factor was the most influential factor of soil erosion due to the flowing water from rainfall and runoff. The *LS* factor was considered from the elevation map of the Lancang–Mekong River basin (Figure 1) and the calculations of Equations (8)–(10). The range of elevation in the study area is from 0 to 6549 m above sea level, and the elevation mean was 3274 m. The basin with high elevation is mainly located in the upper Mekong River basin, and the elevation gradually decreases in the central part of the basin. More than 65% of the natural area has a slope gradient of >9%, and this area is mainly situated in the upper Mekong River basin. Slopes from 10◦ to 70◦ account for approximately 59%. Thus, the results of the *LS* factor were in the range of 0–336 (Figure 4a), and its mean value was 168. In addition, the areas represented by the *LS* values were below 60. The slope is steep, and the slope length is short. The areas with relatively high *LS* values were located in the upper part of the river, while the those with relatively low *LS* values were situated in the central part of the Mekong Delta.

**Figure 4.** (**a**) *LS* factor and (**b**) *C* factor.

#### 4.1.4. Cropping Management Factor

The *C* factor was applied using the *NDVI* analysis from the MODIS satellite images and the calculation in Equation (11). The *C* factor varied from 0 to 0.7 (Figure 4b). The *C* mean and the standard deviation were 0.34 and 0.076, respectively. Most lands in the study area are forests in parts of China, Laos PDR and Cambodia, and they were represented by relatively low values of the *C* factor. Conversely, relatively high values for the *C* factor were shown in the upper Mekong River basin in China, Thailand, and the Mekong Delta.

#### 4.1.5. Support Practice Factor

The values for the *P* factor were determined following the suggestion in Reference [86] (Table 1). The change in *C* values to *P* values was applied with the functions in ArcGIS. The *P* values were 0.5, 0.8, and 1 (Figure 5). Nearly 52% of the *P* values were between 0.8 and 1, and they represent the largest portion. Thus, most areas in the basin are forests and lands with vegetation cover, indicating that soil is protected from agents of erosion. The areas with relatively high and low *P* values corresponded to similar areas for the *C* values (see Section 4.1.4).

**Figure 5.** *P* factor.

#### *4.2. Potential and Actual Soil Erosion*

Soil erosion was divided into two types: potential and actual soil erosion. Potential erosion (*R*, *K*, *L*, and *S*) was defined as a natural erosion process without cropping management (*C*) and support practice (*P*) factors. If potential soil erosion is combined with the *C* and *P* factors, then it can be considered actual soil erosion (*R*, *K*, *LS*, *C*, and *P*). Potential soil erosion was calculated on the basis of four factors by using ArcGIS and GIS techniques. The range of potential soil erosion rate was 5000–25,000 t/km<sup>2</sup> /y (Figure 6a). The average potential soil erosion was 15,000 t/km<sup>2</sup> /y, and the standard deviation was 4623. The findings on spatial distribution demonstrates high-potential soil erosion in most areas in the basin. Thus, all the factors were computed for actual soil erosion (Figure 6b) which is the real-world soil erosion in the Lancang–Mekong River basin. Actual soil erosion was in the range of 700–10,000 t/km<sup>2</sup> /y. The mean actual soil erosion was 5350 t/km<sup>2</sup> /y, and the standard deviation was 1470. Most of the relatively high soil erosion rates were located in the north part of upper Mekong River basin and Mekong Delta. Some parts of Thailand had values close to the mean actual soil erosion. The results of the potential erosion and actual soil erosion manifested notable differences. The potential soil erosion rate was differentiated by the *C* and *P* factors because of the forest and agricultural areas. Hence, the *C* and the *P* factors play important roles in decreasing soil erosion, and they can reduce the effect by 2.5–7 times in the basin. The *C* factor indicates the capability to absorb the impact of raindrops, reduce the velocity and scouring power of runoff and reduce the runoff volume by increasing percolation into soil. Meanwhile, the *P* factor indicates the capability to decrease the amount and rate of water runoff and soil erosion with supporting cropland practices such as cross-slope cultivation, contouring farming and strip cropping.

#### *4.3. Soil Erosion Risk Mapping*

The results of actual soil erosion can be classified into five categories (Figure 7) according to the Soil Erosion Standard Document–Technological Standard of Soil and Water Conservation (SD238-87) [63]. Table 2 shows the soil erosion in the study area ranging from minimal erosion to extreme erosion. Most of the soil erosion in the Lancang–Mekong River basin (45% of the total area) is moderate erosion. However, the soil erosion rate is higher than 5000 t/km<sup>2</sup> /y hence their classification as high erosion and

extreme erosion; the corresponding areas comprise 37% of the total area, including the high-elevation areas in China, the plateau in Thailand, Tonle Sap, and the Mekong Delta. By contrast, a low erosion rate was found mostly in Laos PDR and some parts of Cambodia because of their forest areas.

**Figure 6.** (**a**) Potential soil erosion; (**b**) Actual soil erosion.


**Table 2.** Soil erosion in the Lancang–Mekong River basin.

The analytical results on the correlation between soil erosion rate and all input factors in the RUSLE model using SPSS are shown in Table 3. The hypotheses of all factors were determined on the basis of a 95% confidence (i.e., level of statistical significance). The results were then used to build the multiple linear regression in logarithmic form for the soil erosion rate and all the RUSLE factors of the Lancang–Mekong River basin.

$$\ln(A) = 0.168 \times \ln(\text{R}) + 0.364 \times \ln(\text{K}) + 0.898 \times \ln(\text{LS}) + 0.184 \times \ln(\text{C}) + 0.246 \times \ln(\text{P})\tag{21}$$

β Equation (21) is given by the values of standardized coefficients that are strongly related with all the RUSLE factors of the soil erosion rate. The strongest influencing factor for soil erosion in the study area was the *LS* factor (β = 0.898). Therefore, slope length and slope steepness directly affect soil erosion. In other words, soil erosion likely occurs because of gravity erosion and water erosion in an area.


**Table 3.** Standardized coefficients of factors in the RUSLE model.

**Figure 7.** Soil erosion classification and location of dams in the Lancang–Mekong River basin.

#### *4.4. Estimation of Sediment Deposition Areas*

The assessment of sediment yield or sediment deposition areas in the Lancang–Mekong River basin was computed by modifying the RUSLE model according to Equations (17) and (18). The RUSLE model was determined using the spatially average parameters for the estimation of sediment yield capacity in each sub-basin. The results of the sediment yield capacity for the study area are presented in Figure 8a. Most of the sub-basins have high sediment yield capacities. Some sub-basins have low sediment yield capacity in the central part and the north part of upper Mekong River basin. The size of the sub-basin and the elevations directly result in sediment yield capacity. The average sediment yield capacity (*Sc*) was compared with the estimation of sediment yield (*Sy*) to assess the sediment deposition and sediment erosion in each grid. If the result of *S<sup>y</sup>* was higher than *Sc*, then sediment erosion appeared in each grid cell. Conversely, if *S<sup>y</sup>* was lower than *Sc*, then sediment deposition appeared in each grid cell. The results of sediment deposition and sediment erosion capacities in each grid cell are shown in Figure 8b. Sediment erosion is presented as positive values, whereas sediment deposition is presented in negative values. Potential sediment deposition and erosion are in the range from less than <sup>−</sup>3000 to more than 6200 t/km<sup>2</sup> /y. The mean was 2105 t/km<sup>2</sup> /y, and the standard deviation was 2033. Relatively high sediment erosion occurs along the flow direction of the mainstream, including the north part of upper Mekong River basin, Laos PDR, Tonle Sap and Mekong Delta. Meanwhile, most areas in Yunnan Province, Thailand and Cambodia have high sediment depositions. The sediment deposition and erosion results can be verified using the observed sediment data from the 15 stations along the Lancang–Mekong River basin. The scatter plot of the whole basin, which was based on the observed sediment data and the assessment data on sediment yield from the RUSLE model, shows good results with a correlation higher than 0.9 (Figure 9). The RUSLE model and the technique used to assess sediment deposition and erosion can be applied in the research and prediction of soil erosion and SSY.

−

**Figure 8.** (**a**) Sediment capacity in each sub-basin, (**b**) Deposited and eroded sediments in each sub-basin.

**Figure 9.** Comparative result of the whole basin based on observed data and the data of the modified RUSLE model by using the sediment estimation technique.

#### **5. Discussion**

#### *5.1. Soil Erosion Rate in the Lancang–Mekong River Basin*

This study focused on the assessment of soil erosion rate and sedimentation in the Lancang–Mekong River basin using the RUSLE model and GIS techniques with related available data. The river as the research object had many data limitations, and accessing input data to develop research on soil erosion, sediment yield capacity and sediment transport was difficult. This study attempted to utilize previous research on soil erosion in the Lancang–Mekong River basin [10,49,50,53,79], as no other evidence and information exist on how much the soil erosion rate has changed in this basin. The average soil erosion in the previous research ranges from 1400 to 8500 t/km<sup>2</sup> /y. Our results are in the range near the mean values of the previous research. The spatial pattern of soil erosion occurrence in the north part of upper Mekong River basin is generally consistent with the findings of [10,49,50,53,79], but some spatial soil erosion results differ in other areas, especially in the lower part of Mekong River. Furthermore, we presumed that the different results can be attributed to the computation of the *R* factor (the main factor in soil erosion) which is influenced by differences in rainfall data. Each of the available rainfall data were previously developed for the purpose of individual projects. Nonetheless, if the *R* factor was developed from local rainfall stations in the six riparian countries, then the soil erosion rates can be regarded accurate and be further improved. Meanwhile, the results of descriptive statics in the RULSE model clearly showed that the *LS* factor is the most influential factor for soil erosion and sediment yield in the Lancang–Mekong River basin, especially in the upper Mekong River. Most studies on soil erosion and sedimentation claim that the geographical features of the Lancang–Mekong River basin, such as its steep slopes and the slope length of its hills and mountains, are affected directly by the occurrence of soil erosion in specific areas, and these sediments are transformed when transported along the river [10,32,35,41,47,49–53,79]. Consequently, the mitigation measures currently used to reduce soil erosion need to further consider solutions related with the *LS* factor such as the implementation of check dams and the application of vegetation cover. In order to consider the analytical results on the correlation between soil erosion rate and all input factors in the RUSLE model using SPSS, the *LS* factor is the strongest influencing factor on soil erosion in the study area. Nonetheless, the analytical results may not be quite effective, because the *LS* factor varies greatly in the river basin against other factors. Therefore, this section should be considered by regarding the different geological and geomorphological characteristics of the river basin such as mountains, piedmont and lowland. Moreover, different altitudinal conditions are also important conditions that directly affect the RUSLE input factors. This issue needs to be improved correctly for understanding the influencing factor on soil erosion in each feature of the river basin. In addition, the case study on potential and actual soil erosion verifies the ability of the *C* and the *P* factors to protect and reduce soil erosion. Natural vegetation covers, such as the forests in Laos PDR and Cambodia, can decrease soil erosion at rates greater than those of agricultural areas in Thailand. Hence, if forested areas are transformed into agricultural activities, then the soil erosion rate will increase remarkably, especially in upstream areas [20].

The Soil Erosion Standard Document–Technological Standard of Soil and Water Conservation (SD238-87) [63] was applied in this study to classify the soil erosion rate in the Lancang–Mekong River basin. One of the reasons is that the river has not been evaluated using the standard on soil erosion classification. Previous research used the soil erosion classification in References [63,92]. However, the number of classifications in Reference [92] is lower than that in Reference [63] and, thus, does not correspond with the results of our study. The highest value of severe erosion according to Reference [92] is greater than 3300 t/km<sup>2</sup> /y, while extreme erosion according to Reference [63] is greater than 8000 t/km<sup>2</sup> /y. The values differ considerably in terms of soil erosion classification. We suppose that the soil erosion classification should depend on the researcher's discretion and the suitability of research results until a set of criteria is developed by relevant credible agencies such as the Lancang–Mekong Cooperation (LMC) or the Mekong River Commission (MRC).

#### *5.2. Estimation of Sediment Yield Using the Modified RUSLE Model*

Soil erosion is the initial process of the sedimentation process of a river channel. The Lancang–Mekong River basin faces the challenge of sediment starvation due to the implementation of water development projects, especially hydropower dams. Most studies confirm that sediments have started to decrease continuously because of sediment trapping by hydropower dams [28,44,47,48,93]. Therefore, sediment yield capacity and sediment deposition should be analysed by relevant organizations and the six riparian country governments when drafting the needed solutions. However, the field measurement of sediment aspects is very difficult due to the limitations of equipment and nations' borders in the Lancang–Mekong River basin. Hence, the modified RUSLE model was used for the estimation of sediment yield. This method was also clearly applied to understand the sediment deposition and erosion.

The developed RUSLE model and new technique used to assess sediment yield capacity and sediment deposition areas were appropriate, and the observed sediment data and the sediment yield results from the RUSLE simulation were well correlated despite the limitation of investigating a large field survey area. The consistency between observed sediment data and the RUSLE results can also improve the accuracy of soil erosion prediction and the analyses of sediment yield capacity and sediment deposition. Nevertheless, the observed sediment data from the 15 stations were insufficient for validation, especially for the upper part of the Mekong River and all the sub-basins. This study could only access two stations (i.e., Jiuzhou and Gajiu) in the upper Mekong River basin. If other sediment data regarding the upper Mekong River can be acquired for analysis, then the effectiveness of the RUSLE model with the abovementioned technique can be effectively assessed. Besides, the results of sediment yield in some of the sub-basins may have been overestimated. Problems may have also been caused by the analysis of the RUSLE input factors which are also likely overestimated values. Additionally, the assessment of sediment deposition and erosion using the modified RUSLE model may have led to overestimated results for the sub-basins. A previous study [20] also showed the same trend for Thailand after applying the modified RUSLE model. Therefore, in the application of the method, the abovementioned limitations should be considered for model enhancement. The method proposed in this study is useful in furthering the research and analysis of sediment load at reservoirs and sediment transport in the Lancang–Mekong River basin. Furthermore, the results can be used as basis to understand the physical process of sedimentation in each sub-basin.

The result in Table 4 shows a quite good comparability of the observed and estimated SSY from the RUSLE model. Almost half of the sub-basins were in approximately 5–10% of the percentage error, while the remaining sub-basins were estimated at more than 10% from the observed values. The sub-basins have a high sediment quantity. The modified RUSLE model can be a well-known simulation. Conversely, if sub-basins have a low sediment quantity, the model shows low performance for sediment yield estimation. These causes may occur from two important factors including the spatial resolution of the RUSLE input factors and the features of the river basins. For the spatial resolution in the analysis, this study chose a rather coarse grid (5 km) resolution despite the limitations of the input data sources. The model can be well-captured in some specific areas from the influence of grid resolution. If this study can be applied to a spatial resolution of 1 km, the sediment yield estimation may be improved efficiently [20]. Meanwhile, the features of the river basins directly affect the sediment yield estimation, especially rainfall from changing climate and land-use change from human activity. Most land in the sub-basins, which have greater values of percentage error (10–29%), have changed from forest areas to agricultural areas (among other types), especially Nam Chi, Nam Mun, Nam Songkhram, and Nam Ngum. This issue created inaccuracies in the analysis of the *C* and *P* factors, because the *C* factor was considered from the MODIS satellite image using the remote sensing techniques, and the *P* factor was also estimated from the *C* factor [86]. Furthermore, sub-basins, which are overestimated values, have features without high slopes when comparing with other sub-basins. Hence, the modified RUSLE model may be able to consider areas with better slopes which is quite

consistent with Reference [20]. Totally, these factors may be the causes of the problem in the study of sediment yield estimation in the Lancang–Mekong River basin.


**Table 4.** Comparative results between observed SSY and estimated SSY from model.

#### *5.3. Soil Erosion Impact on Dams*

Soil erosion can negatively affect hydropower dams in the Lancang–Mekong River basin. Sediments resulting from soil erosion can decrease the storage capacity of dams used to generate electricity and for other purposes. The upper Mekong basin, especially in the northern area, is classified as having extreme and high soil erosion, indicating increased vulnerability to soil erosion rate. Dams under construction and planned dams may also face the risk of increased sedimentation once they become operational. The dams located in the central and south parts of upper Mekong River basin are relatively less risky than those in the north part, because soil erosion in those areas have low and moderate erosion levels. Previous studies [10,53,79] obtained results that coincide with our research for the analysis of soil erosion impact on dams in the upper Mekong River basin. Meanwhile, soil erosion in the lower Mekong River basin, especially from the Khorat Plateau (Thailand) to the Mekong Delta, can also generate sedimentation problems due to the high occurrences of soil loss. The agricultural activities in these areas mainly cause the increase in the soil erosion rate. A dam under construction (Don Sahong) and four planned dams (Ban Koum, Phu Ngoy, Stung Treng and Sambor) may be threatened by soil erosion due to the impact of sub-basins in the Khorat Plateau in Thailand, particularly the confluence with the Lancang–Mekong River's mainstream. In addition, extreme soil erosion occurs in the Mekong Delta. Many studies affirm that the Mekong Delta is the most vulnerable area in terms of risk of soil loss [33–36].

The impact of soil erosion on dams in the Lancang–Mekong River's mainstream can be analysed in two parts based on the water sources of the river, namely, the upper Mekong River (with three river areas from Lancang basin) and the lower Mekong River (composed of the northern highlands, Khorat Plateau, Tonle Sap and Mekong Delta). The upper Mekong River basin covers 180,000 km<sup>2</sup> or approximately 24% of the study area, while the lower Mekong River basin covers 570,000 km<sup>2</sup> or approximately 76%. The soil erosion modulus of the upper Mekong River basin is 235.7 t/km<sup>2</sup> /y. Its percentage relative to total soil erosion is approximately 36%, even if this area is smaller than the lower Mekong River basin. The soil erosion modulus of the lower Mekong River basin is 198.2 t/km<sup>2</sup> /y which represents approximately 64% of the total occurrence of soil erosion. The results of the soil erosion modulus can be explained as that the reservoirs located in the upper Mekong River basin are more vulnerable from soil erosion and increased sediment. Consequently, dams are likely to be at risk

of decreasing storage capacity continually. Our results are consistent with the findings of past studies on the impact of soil erosion on dams and sediment trapping. For instance, Reference [47] reported that the sediment trapping rates of dams under construction and the planned dams in the Lancang–Mekong River basin will increase from 51% to 69% due to the high heterogeneity of specific sediment yield in the different parts of the basin, and much higher trapped sediment load is predicted because of soil erosion resulting from socio-economic development. More than 50% of the sediment load (approximately 140 Mt) in the Lancang–Mekong River basin is expected to be trapped annually. Furthermore, more than 60% of the sediment load originates from China's end of the Lancang–Mekong River's mainstream. Existing dams, dams under construction and planned dams are expected to have the highest impact on storage capacity due to the fact of sediment load. Reference [28] reported that the main dams in the Lancang River, such as Manwan, Gongguoqiao, Dachaoshan and Jinghong, have sediment trapping rates between 30% and 70% because of the high sediment yield in the Lancang–Mekong River's mainstream and sub-basins. The storage capacity of reservoirs will continuously decrease from the sediment load due to the soil erosion in the reservoir upstream.

#### *5.4. Delineating Sediment Form*

This study endeavoured to estimate the sediment yield by considering factors such as soil erosion, gully erosion and rill erosion. These erosions are not the only sources of sediment into the river channel, because sediment yield is fundamentally controlled by climatic conditions, geomorphologic characteristics and anthropogenic forcing [22,48]. Some sediments are formed by erosion in the river channel. Our analysis did not take other factors into account in this study. This study could not consider erosion in the river channel due to the limitation of the modified RUSLE model which solely analyses erosion on land. For the study of channel deformations and changing river morphology, a hydrodynamic model is needed. Besides, sediment data (suspended and bed load sediment) for the Lancang–Mekong River basin are insufficient, because a number of measuring stations continue to be unavailable. This is the main limitation for further study in the basin. The results in this research can be considered together with erosion in the river channel using a hydrodynamic model; it would be able to show the sediment process on both land and in the river.

#### **6. Conclusions**

The RUSLE model was integrated with GIS techniques in this study to assess soil erosion and sediment yield in the Lancang–Mekong River basin. The impact of soil erosion on hydropower dams was also considered. The findings indicate that soil erosion occurs in all areas of the Lancang–Mekong River basin, accounting for 5350 t/km<sup>2</sup> /y of its average soil erosion rate or approximately 45% of the basin. The north part of the upper Mekong River basin and some parts of Thailand have higher terrains than the other areas, and they have good vegetation cover and support practice. Furthermore, the *LS* factor showed that this factor was the strongest influencing factor for soil erosion in the study area. The spatial distribution of soil erosion also indicated that the norther part of the upper Mekong River basin and the central and southern parts of the lower Mekong River basin are the most vulnerable areas in terms of increased soil erosion rates due to the movement of sediments to the river. Hence, the dams in this river are highly threatened by sediment problems.

The value of pursuing research on the sediment capacity of each sub-basin of the Lancang–Mekong River basin are summarized as follows. The size of the sub-basins and their elevation directly affect the sediment capacity of the river. Moreover, the spatial distribution of sediment deposition and erosion indicates that relatively high sediment erosion occurs along the flow direction of the mainstream, from the northern part of upper Mekong River basin to the Mekong Delta. The findings on sediment yield estimation from the modified RUSLE model and the observed sediment data were in good agreement and had high correlation. The proposed technique can be applied in the assessment of sediment yield capacity and sediment deposition in the Lancang–Mekong River basin.

The modified RUSLE method was successfully applied to the assessment of the amount and spatial distribution of soil erosion and sediment deposition in the Lancang–Mekong River basin. The method can be applied not only to this river but also to other important areas. This study can help policymakers and relevant organizations improve their decision making based on the provided valuable information on soil erosion and sedimentation in this region.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4441/12/1/135/s1, Table S1. Lists and location of hydropower dams and reservoirs in Lancang-Mekong River's mainstream, Table S2. Average observed SL and SSY from 1962 to 2010 based on catchment area of the station, Table S3. Data sources for the analysis of the RUSLE factors in this study, Figure S1. Locations of the sediment load observed stations.

**Author Contributions:** Conceptualization, P.C., M.X. and W.T.; Methodology, P.C., M.X. and W.T. Software, M.X. and W.T.; Formal analysis, P.C.; Investigation, M.X. and W.T.; Resources, P.C. and W.T.; Data curation, P.C. and W.T.; Writing—original draft preparation, P.C.; Writing—review and editing, M.X. and W.T.; Visualization, P.C.; Supervision, M.X. and W.T.; Funding acquisition, M.X. and W.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China (Grant Nos. 51579135, 51379104 and 51079070), the State Key Laboratory of Hydroscience and Engineering (Grant Nos. 2013-KY-5 and 2015-KY-5), the Chinese Academy of Sciences (XDA23090401) and the National Key Research and Development Program of China (2016YFC0402407).

**Acknowledgments:** Our sincerest appreciation to the Lancang–Mekong Cooperation and the Mekong River Commission for providing us the observed sediment data. We would also like to thank the APHRODITE Project for allowing us the use of the precipitation product. We are also thankful for the technical recommendations of Prem Rangsiwanichpong and the encouragement from my beloved wife, Usa Chuenchum.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## **Estimating Human Impacts on Soil Erosion Considering Different Hillslope Inclinations and Land Uses in the Coastal Region of Syria**

**Safwan Mohammed 1 , Hazem G. Abdo 2 , Szilard Szabo 3 , Quoc Bao Pham 4,5 , Imre J. Holb 6,7 , Nguyen Thi Thuy Linh 8,9, \*, Duong Tran Anh 10 , Karam Alsafadi 11 , Ali Mokhtar 12,13,14 , Issa Kbibo 15 , Jihad Ibrahim <sup>15</sup> and Jesus Rodrigo-Comino 16,17**


Received: 30 August 2020; Accepted: 4 October 2020; Published: 7 October 2020

**Abstract:** Soils in the coastal region of Syria (CRoS) are one of the most fragile components of natural ecosystems. However, they are adversely affected by water erosion processes after extreme land cover modifications such as wildfires or intensive agricultural activities. The main goal of this research was to clarify the dynamic interaction between erosion processes and different ecosystem components (inclination, land cover/land use, and rainy storms) along with the vulnerable territory of the CRoS. Experiments were carried out in five different locations using a total of 15 erosion plots. Soil loss and runoff were quantified in each experimental plot, considering different inclinations and land uses (agricultural land (AG), burnt forest (BF), forest/control plot (F)). Observed runoff and soil loss varied greatly according to both inclination and land cover after 750 mm of rainfall (26 events). In the cultivated areas, the average soil water erosion ranged between 0.14 <sup>±</sup> 0.07 and 0.74 <sup>±</sup> 0.33 kg/m<sup>2</sup> ; in the BF plots, mean soil erosion ranged between 0.03 <sup>±</sup> 0.01 and 0.24 <sup>±</sup> 0.10 kg/m<sup>2</sup> . The lowest

amount of erosion was recorded in the F plots where the erosion ranged between 0.1 ± 0.001 and 0.07 <sup>±</sup> 0.03 kg/m<sup>2</sup> . Interestingly, the General Linear Model revealed that all factors (i.e., inclination, rainfall and land use) had a significant (p < 0.001) effect on the soil loss. We concluded that human activities greatly influenced soil erosion rates, being higher in the AG lands, followed by BF and F. Therefore, the current study could be very useful to policymakers and planners for proposing immediate conservation or restoration plans in a less studied area which has been shown to be vulnerable to soil erosion processes.

**Keywords:** soil management; land cover changes; Syria; soil erosion; hillslopes

#### **1. Introduction**

Soils are vital components of environmental systems and supply livelihoods, services and goods for humans and natural ecosystems [1,2]. Soils are formed by numerous factors such as parent material, topography, climate, water, organisms, and time; however, it is well-known that this process is slow and endangered by land degradation due to certain human activities [3–5]. Intensification of anthropogenic effects has become a key factor that causes negative structural shifts in the soil matrix and health; thus, there has been an acceleration of the erosional cycle from prehistoric times [6] to today [7]. Current soil erosion rates, caused by water or wind, are high and can be considered one of the most serious ecological threats to land sustainability globally [8], given that more than 75 billion tons per year of soil are lost due to soil erosion [9]. The problem associated with soil erosion by water is the result of spatial integration of physical and human factors, which vary significantly across scales (from pedon to watershed), and make any estimation difficult [10–12]. Soil erosion irremediably reduces the quality of the physico-chemical and biological properties, soil fertility and land productivity, which considerably affect cultivated areas [13,14]. Therefore, nature-based solutions to achieve land degradation neutrality can be a key to conserving ecosystem services [15,16]. However, for any ecological restoration, stakeholders and land managers must be fully motivated and convinced, and this is still a current challenge [17,18].

Soil erosion is progressively limiting the availability of resources, threatening biodiversity, and affecting food production, and is accelerated by specific drivers such as climate change, land use/land cover changes, overgrazing, inappropriate farming procedures, or armed conflicts [19–23]. Consequently, soil erosion is defined as a physical and anthropological challenge [24]. During the 1980s, statistics indicated that about two billion hectares of agricultural land had completely deteriorated since 1000 AD, and currently, the FAO estimates ~75 billion tons of agricultural soil loss, causing an annual cost of USD 400 billion [25]. Consequently, an increasing interest in soil stability and conservation is progressively evolving to deal with this worldwide environmental problem in the context of the landscape changes which have occurred in the current century [9,26,27].

Soil erosion is the outcome of the dynamic interaction between different ecosystem components, e.g., land use, inclination, rainfall intensity, and soil properties [28]. The mechanism of soil erosion by water includes splashing and detachment of soil particles due to the kinetic energy of raindrops, then the transportation of these particles by surface runoff [29]. However, due to its tremendous impact, recent research has been directed more towards erosion control techniques in many parts of the world, for example in Austria [30], Spain [31], China [32], Hungary [33], Germany [34], and France [35], among others.

The components of the Mediterranean environment are considered one of the most fragile around the world, especially as regards soils, which are exposed to severe degradation processes [36,37]. In several cases, rainfall and runoff have induced soil erosion, which is a well-known degradation challenge in terms of ecological mismanagement [8,38]. Rugged and dissected terrains, steep slopes, high rainfall intensities, shallow and skeletal soil thicknesses, receding and sparse vegetation, and

chronic and severe drought stress in summer are among the most important physical factors which drive soil erosion [39,40]. Several authors have reported that the annual rates of soil erosion have reached dangerous levels, exceeding the allowable soil loss tolerance limits (2 to 12 Mg ha−<sup>1</sup> year−<sup>1</sup> ) for agricultural and economic sustainability in the Mediterranean environment [41–45]; nevertheless, these numbers can vary depending on the main goal of the research and the specific area [46]. For example, Kouli et al. [47] determined that more than 1 Mg ha−<sup>1</sup> year <sup>−</sup><sup>1</sup> may be irreversible within a time dimension ranging from 50 to 100 years.

Syria is among the eastern Mediterranean countries which are seriously exposed to the problem of water-related soil erosion, especially in the coastal region of the country (CRoS). This area represents an appropriate terrestrial, structural, climatic, hydrological, and intense anthropological case of the acceleration of soil erosion by water [23,48]. In the CRoS, soil erosion by water is the first threat to agricultural activity, which is the pivot of economic life for 34.8% of the population [49]. Meanwhile, CRoS is considered the first agricultural stability zone in Syria, receiving more than 600 mm of rainfall and being used for rainfed agriculture, with a total agricultural land area of 2.7 million ha [50]. Accordingly, the issue of soil erosion in CRoS has been assessed by many local scholars at the administrative area or catchment area level, using different models such as the Revised Universal Soil Loss Equation (RUSLE) [51,52], the Water Erosion Prediction Project (WEPP) model [48,53], and the Coordination of Information on the Environment (CORINE) model [54]. On the other hand, a limited number of studies have dealt with soil erosion after wildfires. Al-Ali and Kheder [55] stressed the importance of monitoring soil erosion after wildfires, where the soil erosion from burnt forests reached 7.22 Mg ha−<sup>1</sup> year−<sup>1</sup> .

In the CRoS, as well as in the Mediterranean region in general, different anthropogenic activities (i.e., rapid changes in land use driven by intense population pressure or agricultural expansion) and climate change have rapidly exacerbated soil water erosion. However, information about soil erosion on the field-scale in the near-eastern Mediterranean remains limited compared to that in the western and northern Mediterranean. Some representative examples can be found in the territories of border countries, highlighting the importance of assessing land degradation processes from different points of view, e.g., [55–58]. Within this perspective, the main aim of this research was to bridge the gap in the common literature on soil water erosion in the coastal territories of Syria by measuring soil water erosion and runoff under three different land uses (agricultural land (AG), burnt forest (BF), forest/control plot (F)). Our hypotheses were the following: (i) agricultural areas are the main areas at risk of soil loss; (ii) burnt forests are endangered by the increased runoff and severe soil loss; (iii) the effect of inclination on erosion rates has a saturation curve, i.e., above a threshold inclination, the rate of erosion does not increase relevantly; and (iv) slopes and land cover have a significant interactive effect, thus, these two factors determine erosion hazard together.

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

#### *2.1. Study Area*

The study area is located in the western part of Syria (35◦49′ to 36◦31′ E; 34◦49′ to 36◦05′ N) within an area of 5274 km<sup>2</sup> (Figure 1). The elevation of the region ranges from 0 to 1700 m a.s.l. The Syrian coast area consists of three basic geomorphological units: the plain (0–100 m), the plateau (100–400) and the mountains (400–1700) [59]. The study area is characterized by narrow plains near the coast, followed by dissected mountains. The degree of inclination generally ranges from 0◦ to more than 60◦ . The coastal strip was affected by recent tectonics, which caused a fluctuation in the sea level from the Early Pleistocene to the recent "upper" Holocene. This is reflected in the diversity of rock formations such as sandstones, sands and conglomerates, which were laid down as sedimentary deposits with limestone and marls. Interestingly, these rocks were penetrated by basaltic rocks in the southern part of the coastline [60].

**Figure 1.** The coastal region of Syria and locations of the experimental plots.

According to the Köppen climate classification, the study area falls into two categories (Csa and Csb) with the main group being C, which follows the Mediterranean climate. The rainy winter season is mostly concentrated between November and March [61]. In general, the average rainfall ranges from 765 mm (near the coast) to 1250 mm (in the high mountains) [61]. The mean annual temperature in the plain areas is about 19.3 ◦C and in the mountains, it is about 14.8 ◦C [61]. The common soil orders are Inceptisols, Entisols, and Mollisols [62]. The study area includes the governorates of Tartous and Lattakia, with a population of approximately three million [61]. Syria is divided into five agricultural stability zones, according to distributed rainfall and the suitability for rainfed agriculture. The study area is located in the first agricultural stability zone, where rainfall exceeds 600 mm [63].

Traditional agriculture is the most essential economic axis for rural inhabitants, and most fields are cultivated with wheat, and olive and citrus orchards. Between 2010 and 2018, more than 800 wildfires were recorded in the coastal region of Syria. Wildfires usually occur between June and late August (summer season), and are typically induced by human activities. In this research, the experimental burnt sites were selected based on fire time and intensity.

#### *2.2. Experimental Design*

Based on a field survey conducted in the study area, five different locations (SY1, SY2, SY3, SY4, SY5) with different hillslope inclinations (38%, 45%, 15%, 29%, 10%) were chosen as representative sites for measuring soil erosion (Table 1). Three different land uses were selected at each location: (i) agricultural land (AG), where traditional cultivation, sowing, and harvesting occurs, with the absence of any mechanization; the common crops in AG lands are wheat (SY1, SY2), olives (SY3, SY4), and citrus orchards (SY5); (ii) burnt forest (BF), where soil cover varies from 30% to 55% with local natural vegetation; and (iii) forest land (F), which is characterized by mixed forest, and is used as a control plot without recently extensive human disturbance. The soil cover for all treatments was sampled without any disturbance.


**Table 1.** Experimental characteristics of the five locations studied.

Experimental plots of 2 × 1.6 m were installed with metal barriers of 0.5 m height (0.15 m into the soil) to collect runoff and soil loss. This method was previously adopted in Syria by [64], and applied by [61,65] and [61,66]. Nonetheless, the plots designed were similar to [61,67], but of a smaller size.

The amount of rainfall (mm) was measured on-site by placing a metal rainfall gauge at each location. Meanwhile, runoff (L/m<sup>2</sup> ) was recorded at each plot after each rainy event by recording the volume in each sediment collector. Soil loss (kg/m<sup>2</sup> ) was also determined by mixing the collected soil detachment and a representative sample of 5 L each. Finally, the samples were transported to the laboratory. In the laboratory, each sample was placed in a small container and dried in an oven (105 ◦C) for 24 h.

In addition, soil samples were collected at the beginning of the monitoring period from the topsoil (0–0.15 m) in each plot, and soil texture and soil erodibility factors were determined (Table 2). The design and performance of the chosen experimental plot with the sampling strategy were tested following [39] (Figure 2). Data were collected from October 2012 to December 2013 (i.e., the vegetation period). A total of 26 rainy storms were observed during the monitoring period. ′ ′ ′ ′

**Table 2.** Soil texture and K value in the studied locations (SY1-SY5) for three land uses (agricultural land (AG), burnt forest (BF), and forest land (F)). ′ ′ ′ ′


**Figure 2.** Sketch design for the experimental plot.

#### *2.3. Data Analysis*

Average, maximum, minimum, and median values were determined. Soil erosion and runoff data were depicted in boxplots, together with the linear regression among them in each land cover class. Normal distribution was checked by the Shapiro–Wilk test (S-W); as this failed, the non-parametric Kruskal–Wallis test (K-W) [68] was applied as an alternative to the one-way ANOVA. The K–W test aimed to detect the difference between the medians of the treatments with the following hypothesis: *H***<sup>0</sup>** was that the medians of the studied groups were from the same distribution, while *H***<sup>1</sup>** represented the idea that the medians of the studied groups were different. As the K–W test did not show which plot is statistically different from any other, the pairwise comparison among slopes was performed with the Mann–Whitney test with Bonferroni correction. Pairwise analyses in the same hillslope but for different land uses (i.e., AG-F; AG-BF, BF-F) were neglected as we focused on the differences caused by inclination and did not analyze the obvious differences among land use types. Finally, to assess the relationships between the studied variables, a correspondence analysis was carried out. We applied a General Linear Model (GLM) to reveal the importance of rainfall, inclination and land use types. The inclination type was included as ordinal data and land use as a categorical dummy variable. We determined the model parameters, and the effect sizes expressed as partial η <sup>2</sup>p, which expressed the contribution of each variable and the interaction of the factors as a standardized measure [69]. The effect can be very small (η <sup>2</sup>p < 0.01), small (0.01> η <sup>2</sup>p > 0.06), medium (0.06 > η <sup>2</sup>p > 0.14), and large (η <sup>2</sup>p > 0.14). Differences among inclination degrees were analyzed with the t-test and ANOVA using the 1999 Monte-Carlo permutation. Statistical analyses were conducted with SPSS (v24; IBM, Chicago, IL, USA), the EViews software package (v10; [70] New York, NY, USA), and R 3.6.3 [71] with the gamlj package [72].

#### **3. Results**

#### *3.1. Soil Water Erosion and Runo*ff

Observed runoff and soil erosion varied according to both inclination and land use, as can be observed in Appendix A (Figure A1). The total rainfall in the study area exceeded 750 mm, divided into 26 events. The average soil loss ranged between 0.74 <sup>±</sup> 0.33 and 0.14 <sup>±</sup> 0.07 kg/m<sup>2</sup> , while runoff ranged between 42.14 <sup>±</sup> 15.27 and 12.77 <sup>±</sup> 5.84 L/m<sup>2</sup> in the AG (Table 3). Meanwhile, in the BF plots, mean soil loss ranged between 0.24 <sup>±</sup> 0.10 and 0.03 <sup>±</sup> 0.01 kg/m<sup>2</sup> , and runoff from 22.95 <sup>±</sup> 9.33 to 3.77 <sup>±</sup> 1.62 L/m<sup>2</sup> . The lowest amounts of soil loss and runoff were recorded in the F lands, where the ranges were between 0.07 <sup>±</sup> 0.03 and 0.1 <sup>±</sup> 0.001 kg/m<sup>2</sup> , and 11.98 <sup>±</sup> 4.73 and 1.78 <sup>±</sup> 0.78 L/m<sup>2</sup> , respectively.


**Table 3.** Univariate statistics of observed soil loss and runoff in the studied locations (SY1-SY5) under three land uses (AG: agricultural land, BF: burnt forest, F: forest).

Min: Minimum; Max: Maximum; SD (n): Standard deviation (n); ϕ: Standard error of the mean.

The highest soil loss was 1.34 <sup>±</sup> 0.33 kg/m<sup>2</sup> , registered in the AG lands with 38% inclination, and the lowest was in the gentlest slope (10%), reaching 0.29 <sup>±</sup> 0.07 kg/m<sup>2</sup> (Figure 3a; Table 3). Soil loss was the highest in both BF and F with a hillslope inclination of 45%, reaching 0.45 ± 0.10 and 0.13 <sup>±</sup> 0.03 kg/m<sup>2</sup> , respectively (Figure 3b,c; Table 3).

**Location** 

**L/m2**

η

η η

η η

**Figure 3.** Box plots of soil erosion in each ecosystem (with respect to slope): (**a**) agricultural land; (**b**) !burnt forest, and (**c**) forest (median (\_\_\_\_\_); mean (•); median 95% confidence (shaded).

Similarly, a maximum runoff was recorded in the AG lands with 72.5 L/m<sup>2</sup> in the steepest slopes (Figure 4a). In BF, the highest runoff was 41.51 L/m<sup>2</sup> with 45% inclination; meanwhile, the lowest was observed with the gentlest slopes, reaching 1.1 L/m<sup>2</sup> (Figure 4b). In F lands, the highest runoff was 21.50 L/m<sup>2</sup> (SY1) and the lowest was 3.50 L/m<sup>2</sup> (SY5) (Figure 4c).

**Figure 4.** Box plots of runoff in each ecosystem (with respect to slope): (**a**) agricultural land; (**b**) burnt forest, and (**c**) forest (median (\_\_\_\_\_); mean (•); median 95% confidence (shaded).

**Location** 

**Location**

In each studied land-use type, regression analysis detected a high correlation between the generation of runoff and the activation of soil loss: R<sup>2</sup> values were 0.91, 0.87, and 0.89 (p < 0.05) in AG, BF, and F, respectively (Figure 5).

**Figure 5.** Correlation between soil loss (kg/m<sup>2</sup> ) and runoff (L/m<sup>2</sup> ) (regardless slope): (**a**) agricultural land, (**b**) burnt forest, and (**c**) forest.

#### *3.2. Impact of Inclination on Soil Water Erosion*

The Kruskal–Wallis test (K–W) showed that at least one of the studied plots was significantly (p < 0.05) different from other treatments in each slope inclination (SY1, SY2, SY3, SY4, SY5), and land use (AG, BF, F) (Table 4).


**Table 4.** Kruskal–Wallis analysis in each ecosystem for both soil water erosion and runoff (p < 0.05).

The significance level is 0.05.

The pairwise comparison among the inclinations showed that there was a significant difference (p < 0.05) among them in the agricultural lands, both in the case of soil loss and runoff under different inclinations (Table 5). Differences were also significant (p < 0.05) between 15% (SY3) and 45% (SY2), and between 15% (SY3) and 38% (SY1). However, non-significant differences were noticed among the following plots: 10% vs 15%; 29% vs 45%; 29% vs 38%; and 45% vs 38% (Table 5). Just as with the AG, the F plots showed similar values with one exception: in the 29% vs the 38% plots, where the difference was significant regarding the runoff. In BF plots, significant differences were recorded in soil loss data among the following pairs: 10% vs 29%; 10% vs 38%; and 10% vs 45%; meanwhile, the pairwise analysis of runoff data showed identical results in the F plots.

**Table 5.** Pairwise comparisons between slopes for soil loss and runoff for three land uses (AG: agricultural land, BF: burnt forest, F: forest).


Each row tests the null hypothesis that Sample 1 and Sample 2 distributions are the same. Asymptotic significances (2-sided tests) are displayed. The significance level is 0.05. <sup>a</sup> Significance values have been adjusted by the Bonferroni correction for multiple tests. The bold numbers and bold color express the significance (p < 0.05).

The correspondence analysis revealed that erosion and runoff in the F plots can be discriminated and highly differentiated from both AG and BF, as it was located in a position further from the origin (x = 0, y = 0), whilst AG and BF were less distinct (Figure 6a). Similarly, erosion on 45% hillslope inclination, followed by 38%, was differentiated from other hillslope inclinations, while the 10%, followed by 45%, and 38% were differentiated from other slope inclinations in terms of runoff (Figure 6b).

**Figure 6.** Correspondence analysis per plot: (**a**) soil loss and (**b**) runoff.

#### *3.3. Multivariate Analysis of Factors and Covariates*

The GLM revealed that all the factors involved (inclination and land use) and the covariate (rainfall) had a significant (p < 0.001) effect on the soil loss, and the explained variance was 85.1% (based on the adjusted R<sup>2</sup> = 0.851). Furthermore, the statistical interaction also obtained a significant (p < 0.001) effect (Table 6). Regarding the relevance of the predictors, land use registered the largest effect, while the effect of rainfall was 40% smaller, and the inclination effect was about half. The effect of the interaction of inclination and rainfall was similar to the rainfall effect.

**Table 6.** Summary of the General Linear Model (GLM) performed with soil erosion as the target variable (SS: sum of squares, df: the degree of freedom, F: F-statistic, p: significance, η <sup>2</sup>p: effect size; p < 0.001 is highlighted in bold).


Generally, the soil loss rate of the AG lands was the greatest in all hillslopes, while the control areas (F) had the lowest rate. The erosion can be regarded as linear in these areas; locally estimated scatterplot smoothing (LOESS) curves were almost linear in all possible combinations (Figure 7). Visual evaluation of the data showed that inclination degrees can be divided into two different groups based on the soil loss: (i) inclination of 10 and 15%, and (ii) 29, 38 and 45%. In the case of group (i), the erosion rate was below 0.5 kg/m<sup>2</sup> , although the difference between the AG lands was significant (mean difference: 0.096; pM-C < 0.0005). Larger differences were caused by the heaviest rainfalls in the study area with 15% inclination. Erosion rates within group (ii) were similar regarding all the three land cover types, and the soil loss in the 45% inclination area was not significantly different (p > 0.05) from the 38% or 29% inclination degree areas according to the ANOVA test (F = 2.059, df = 2, p = 0.129).

**Figure 7.** Relation between erosion rates and rainfall by slope and land cover type (agricultural land: AG; burnt forest: BF, and forest: F).

#### **4. Discussion**

Soil erosion by water is considered one of the most important agricultural sustainability challenges in the CRoS as a result of the following factors: heavy rainfall, severe inclinations, high erodibility, massive gushes of runoff, land-use changes, and non-sustainable agricultural practices [48]. Therefore, the assessment of water erosion derived from field analysis provides a detailed method of approaching the relationship between erosion, runoff and soil properties.

#### *4.1. Criteria to Assess Current Erosion*

#### 4.1.1. The Role of Physical Features in Erosion

The climate of the study area is characterized by a high-intensity precipitation pattern with the intense kinetic energy of raindrops that hit the hillslopes with different land uses. Land use played an influential role in determining the quantities of eroded material and discharged runoff, which varied according to other physical features such as topography and soil properties [73]. Recently, forest lands in the CRoS were badly affected by severe wildfires, which increase the susceptibility to soil loss in the study area. In this regard, the importance of soil management was clear. In our research, soil loss and runoff were the highest in the AG and BF plots compared to the F plots. Within the study area, the cultivated land (AG plots) and burnt plots remained bare and exposed directly to raindrops, which could explain the high amount of soil erosion and runoff in comparison to the F plots, as other recent investigations in cultivated or abandoned fields have demonstrated [14,74], or in areas after recent wildfires [75].

#### 4.1.2. The Role of Slope Steepness in Erosion

Of the five locations used for measuring soil loss and runoff, three of them were chosen with an inclination higher than 25%, i.e., SY1, SY2, SY4. Our statistical analysis revealed that from 29%, the critical limit was to be found above this value, similarly to a saturation curve, in that a greater slope gradient did not cause a relevant increase in the erosion rate (Figure 8). These results agree with other soil erosion and runoff reports presented by [76–79]. In the light of the high-intensity rainstorms in the study area, inclination was also a driving factor in the occurrence of high-velocity runoff events which enhance soil detachment. Additionally, this inclination could motivate both ponding depth and depressional storage [80–82]. Under the same land use, inclination degree could accelerate the erosion remarkably, as can be revealed from Table 3 and Figure 7. Land use had a relevant effect on the soil erosion rate, with the highest values observed in the AG plots, and the lowest ones in the F plots. In Syria, only a few studies have reported on soil erosion at the plot scale. Barneveld et al. [83] claimed that soil erosion in the NW part of Syria rarely exceeded 5 kg/m<sup>2</sup> in cultivated olive lands with average slopes of 25%. However, this difference in measuring soil erosion could be explained by the physiographic difference between each research location, especially the steepness of the slope, the form and development of terrain, precipitation intensity, soil characteristics, and agricultural practices. Notably, our results are higher than the erosion observed in Mediterranean mountains by [84] (147.3 g m−<sup>2</sup> ) and lower than results reported by [83] (5 kg/m<sup>2</sup> ).

#### 4.1.3. Role of Human Activities in Erosion

If the physical factors are compared to human activities, the latter is the main driver of erosion through poor and unsustainable soil management and tampering with soil structure, and altering its physical, chemical and biological properties, especially its organic matter content [85]. However, these consequences should be considered as serious in fragile and vulnerable soils as in the Mediterranean environment. Intensive tillage and bare soils play a key role in accelerating soil erosion [86,87] by enhancing the separation of macro-aggregates, which negatively affects the soil aggregates' stability [88–92]. Soils in forest plots are protected by more vegetation and we hypothesize that soil aggregates are stronger and are not affected by the negative impacts of the kinetic energy of raindrops.

Some authors have observed that the collapse of soil aggregates can minimize soil porosity by blocking pores by fine particles (silt, clay) and can magnify soil sealing and crusting, and, subsequently, soil erosion can be enhanced [93]. As a consequence, some authors have even reported that the soil erodibility factor (K) is higher in AG plots for this reason, which indirectly indicates the susceptibility of AG plots to soil erosion [94]. In this regard, organic matter (OM) is expected to be higher in the F plots, which significantly enhances aggregate stability against rainy storms, while aggregates in AG and BF plots would be more vulnerable [95–97]. Our results are consistent with [98], who indicated that inappropriate agricultural practices in shallow topsoil can increase the susceptibility of runoff and erosion. This is extensive in various agricultural activities in the Mediterranean belt [99,100]. The relevance of OM content in mitigating erosion has been proved by several authors, i.e., soils with <2% OM are highly susceptible to erosion and runoff [101,102]. In addition, [36,103–105] highlighted the vital role of agricultural activities and the Mediterranean climate in accelerating soil erosion in semi-arid regions, while other studies stressed the importance of ground soil cover for preventing erosion and runoff [99,106–109]. As extensive fieldwork in the CRoS has revealed, in the AG plots there is an absence of most of the agricultural practices that conserve soil, especially crop rotation, maintaining tillage, contour and strip farming, grass water channels, and diversion structures. −

**Figure 8.** Slope gradient and erosion rate (\_\_\_\_\_LOESS fit line with 95% confidence intervals).

#### *4.2. Dimensions of the Current Evaluation*

The CroS constitutes the first agricultural stability zone and the agricultural and economic backbone of the local population, and therefore the protection of its natural resources, especially soils from erosion, is a priority in the framework of agricultural sustainability. Thus, the implementation of some conservation practices (CP) or even the establishment of a national action plan for soil conservation to repair local ecosystems is a high priority. Some authors have recommended CP including soil mulching [110,111], tillage reduction [112,113], buffer strips and minimum cultivation [114] or a correct planification of soil terraces [115]. Nonetheless, field analysis of soil erosion is at the forefront of the measures that will develop strategies for preserving farmland, especially during the ongoing war that has negatively affected the agricultural and food system in the country. In the context of soil erosion, cultivated hillslopes in the CRoS are subject to intensive use pressure which includes poor maintenance technologies, and overuse of fertilizers. Consequently, soil aggregates are more dynamic once there are other agents of erosion, especially high-intensity raindrops. In this regard, the orographic precipitation model imposes high rainfall intensities, and consequently massive runoff which accelerates soil erosion. Unfortunately, in this research, the intensity and duration of rainfall could not be measured. However, further studies should address those elements instead of using the total rainfall amount per event. In

addition, further research should be carried out to address appropriate measures for land conservation, especially with hillslopes of over 29% inclination.

### **5. Conclusions**

In this research, soil loss and runoff were measured in five different locations (hillslopes) with three different land uses (AG, BF, F) in the coastal region of Syria. The main findings of this research are:


Few studies have dealt with soil erosion in Syria, and to our knowledge none of these have measured erosion per rainfall event at different hillslope positions comparing human disturbances under three different land uses. The outcome of this research could play an important role in setting up the first conservation plan in Syria. Moreover, the output of this research will contribute to bridging the gap in the common literature on soil water erosion in the near-eastern Mediterranean, and could be used for the improvement of erosion equations or soil protection policies not only in Syria, but all over the Mediterranean belt.

**Author Contributions:** Conceptualization, S.M., I.K. and J.I.; Data curation and collection, I.K. and J.I.; Resources, I.J.H.; Supervision, S.S. and J.R.-C.; Visualization, K.A., Q.B.P., N.T.T.L., D.T.A., A.M. and I.J.H.; Writing—original draft, S.M., and H.G.A.; Writing—review and editing, All authors. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** This paper is part of a research project of the first author (Safwan Mohammed) funded by the Tempus Public Foundation (Hungary) within the framework of the Stipendium Hungaricum Scholarship Programme. The research was supported by the Thematic Excellence Programme of the Ministry for Innovation and Technology in Hungary (ED\_18-1-2019-0028) projects. Authors would like to thank Ministry of Local Administration and Environment (Syria), Tishreen University (Syria) and Debrecen University (Hungary) for their unlimited support.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
