**The Impact of Vegetation Successional Status on Slope Runo**ff **Erosion in the Loess Plateau of China**

**Enhao Chang 1,2 , Peng Li 1,2, \* , Zhanbin Li 1,3 , Yuanyi Su 1,2 , Yi Zhang 1,2 , Jianwen Zhang 1,2 , Zhan Liu 1,2 and Zhineng Li 1,2**


Received: 4 November 2019; Accepted: 7 December 2019; Published: 11 December 2019

**Abstract:** Slope vegetation restoration is known to influence erosion in the Loess Plateau region in China. The ability of vegetation to mitigate soil erosion under extreme runoff, however, has not been studied in great detail in this region. Here, we examine five typical vegetation communities in the Loess Plateau region that originated from restoration efforts enacted at different times (1, 11, 15, 25, and 40 years). Water scouring experiments were carried out to monitor vegetation community succession and its effects on erosion. These results indicate that the sum of plant importance values increased from 260.72 to 283.06, species density increased from 2.5 to 4.5 per m<sup>2</sup> , and the amount of litter and humus increased from 24.50 to 605.00 g/m<sup>2</sup> during the 1 to 40 years of vegetation community succession. Root biomass and root diameter reached a maximum of approximately 10.80 mg·cm−<sup>3</sup> and 0.65 mm at 40 years of recovery. Slope runoff velocity decreased by 47.89% while runoff resistance increased by 35.30 times. The runoff power decreased by 19.75%, the total runoff volume decreased by 2.52 times, and the total sediment yield decreased by 11.60 times in the vegetation community. Slope runoff velocity and power had the largest correlation with aboveground vegetation (0.76, 0.74), total runoff had the largest correlation with underground roots (0.74), and runoff resistance was most strongly correlated with soil structure (0.71). Studies have shown that the succession of vegetation communities can enhance the aboveground ecological functions of plants, thereby significantly reducing the runoff velocity and power. The development of plant root system significantly reduces the runoff volume; the improved soil structure significantly increased the runoff resistance coefficient.

**Keywords:** vegetation community; vegetation importance value; root system; soil erosion; grey correlation analysis

#### **1. Introduction**

In recent years, there has been a series of studies conducted on soil and water conservation focusing on silt-dam gully engineering, terraced fields of slope engineering, and the Grain for Green Project on the Loess Plateau [1–3]. These efforts play an important role in the ecological restoration of Loess areas. The average annual sediment in the Yellow River has been reduced from an estimated 1530 million tons in the 1950s to 166 million tons in the 2010s [4]. About 40% to 50% of the reduction of the average annual sediment in the Yellow River comes from soil and water conservation measures in the Loess Plateau [5]. Since 1999, the Chinese government has carried out the Grain for Green Project, which has restored slope vegetation along the Loess Plateau. Observations from remote sensing show that the vegetation coverage of the Loess Plateau increased from 31.6% in 1999 to 59.6% in 2013 [6]. Restored forest and grassland areas accounted for about 56.7% of the total area of the Loess Plateau [7]. Some have suggested that vegetation restoration on the Loess Plateau has resulted in a 50% reduction in sedimentation along the Yellow River [8,9]. Ecological restoration of vegetation thus plays an important role in reducing slope soil erosion in the region [10]. Rainfall, however, does not automatically generate runoff. Erosion caused by a few short-duration heavy rainstorms can account for more than 60% to 90% of the total annual erosion [11]. Loess slopes can be damaged by slope runoff when extreme rainstorm events are frequent [12]. The slope is also the pioneer path of sediment production, which has a great contribution to the total amount of sediment observed in the outlet section of the watershed in the Loess hilly region [13]. As such, a better understanding of how vegetation can mitigate slope runoff and sediment under the erosion action of high-intensity slope runoff on the Loess Plateau is urgently required.

The Grain for Green Project facilitates a great change in landscape patterns in a certain sense, with the vegetation changing from annual crops to perennial native plants on the Loess Plateau [1]. The vegetation community, in the process of recovery without intervention, began to take place as species succession because of the extension of the years. Successional dynamics in the Loess Plateau typically follow a pattern of Artemisia plants, giving way to perennial rhizome grasses, and finally to perennial arbuscular herbs [14]. Gramineae, Legume, and Compositae occupy an important position in the natural succession of plant communities [15]. Previous studies have found that leguminous plants, such as *Lespedeza davurica*, can improve the soil organic matter, total nitrogen, total phosphorus, and available nitrogen content [16]. The Legume plants gradually became the dominant species, and soil organic carbon, in the 0–50 cm depth soil, gradually recovered after vegetation succession for about 20 years on Loess Plateau [17]. The Asteraceae and Gramineae plants are mainly represented by *Artemisia capillaris*, *Artemisia sacrorum*, and *Bothriochloa ischaemum*. They can grow naturally in the early stage of vegetation restoration, and their root systems have a strong ability to retain surface soil [18–20]. During succession, plant biomass, ground coverage, root structure, and function will change. The soil physical and chemical properties undergo predictable dynamics as well [21–23]. The number of species and individuals increases rapidly from 1 to 10 years in disturbed vegetation communities on the Loess Plateau [24]. The maximum root length density reached 31.04 mm/cm<sup>3</sup> in the 0–20 cm depth soil at 15 years after abandonment, and maximum root biomass density reached 3.35 mg/cm<sup>3</sup> after 21 years. Likewise, the water absorption capacity and the turnover frequency of root systems gradually increased in the process of vegetation restoration [17]. The slope is a basic unit of erosion, and erosion is primarily driven by the hydrodynamic index of the runoff and the erodibility of the topsoil [25]. The flow rate, flow velocity, and runoff depth are key factors directly affecting the slope separation process, and are also the basic parameters used for calculating other hydrodynamic indicators, such as runoff shear, runoff power, Reynolds number, and Froude number. These, in turn, are affected by factors, such as the underlying surface and vegetation [26,27]. Therefore, studying slope runoff hydrodynamics across vegetation communities in varying degrees of succession can help to reveal the role that the successional status has on mitigating slope erosion.

To study this, we studied slope runoff and erosion in five areas that had vegetation in varying degrees of succession (i.e., from 1 to 40 years after restoration) by the method of artificial simulation. Our goal was to analyze the effects of the vegetation successional status on slope soil erosion under water flushing conditions. We expect that results from this study will provide scientific guidance for future research on water erosion dynamic mechanisms and vegetation regulation principles in the Loess Plateau.

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

#### *2.1. Study Area*

The study area is located in the Xindan watershed on the Loess Plateau in China (E 110◦15′–110◦20′ , N 37◦27′–37◦32′ ; 810–1120 m a.s.l.) (Figure 1). The soil type is to yellow loess soil, and the plough layer is thin (20–30 cm). Xindian watershed is a national soil erosion control area, which has been banned for about 60 years. The average annual temperature and precipitation in the study area is 9.7 ◦C and 486 mm, respectively. Precipitation does, however, vary widely across years and space. In recent years, there has been a notable decrease in precipitation days, an increase in the frequency of heavy rains, and an increase in both droughts and floods. ′ ′ ′ ′

The successional sites included the *Artemisia capillaris* for 1 year since restoration, *A. sacrorum* for 11 years, *Bothriochloa ischaemum* for 15 years, *Lespedeza davurica* for 25 years, and *Ziziphus jujube* for 40 years. These species are dominant because of succession. It is also an inevitable sequence of vegetation succession in the Loess Plateau. This method of selection can be regarded as a method of a spatial sequence equivalent to a vegetation succession time series.

**Figure 1.** (**a**,**b**) Study area in Shaanxi Province, China; (**c**) Digital elevation model; (**d**) Digital image and sampling sites.

#### *2.2. Experiments and Tests*

The five experimental plot areas were 4 m × 0.5 m, and slope steepness ranged from 8 to 9◦ . In order to prevent lateral seepage of the slope flow during the test, plots were separated by a 2-mm thick steel plate. We installed flow-stabilizing devices and jet grooves at either ends of the plots, and dug a circular pit under the catchment groove where the sample collecting barrel was placed (Figure 2). Each scouring experiment lasted for 30 min, and the observation sections were set at 1, 2, 3, and 4 m of the plots to observe the runoff in sections. According to the series of precipitation data of the hydrological station in the past 30 years, the P–III frequency of rainstorms was calculated for each duration (10, 20, and 60 min; 3, 6, 12, and 24 h; and 3 d). In the torrential rains of different durations, the 20-year return period short-duration (60, 20, and 10 min) rainstorm intensities reached 0.9, 1.9, and 2.7 mm/min, respectively (Figure 3). As such, we used a flow rate of 4, 8, and 16 L/min for scouring, which is similar to a rain intensity of 2, 4, and 8 mm/min according to the catchment area of the plot.

Runoff and sediment samples were collected every minute, and runoff data was measured every 2 min. Sediment samples were left to settle first, then dried, and measured for sediment yield. We cut off the vegetation on the ground in each plot, retaining the stem of a certain height (5 cm) and the root system (Figure 2). We used the potassium permanganate stain tracing method to measure the runoff velocity on the slope. Runoff depth and width were measured using an artificial ruler. Runoff depth was used as a reference, and the hydrodynamic calculation uses the formula derivation value. Each experimental plot was scoured three times with different flow rates. A total of 45 experiments were conducted in this study.

**Figure 2.** Device schematic diagram and field photos of the experimental process in the experimental area. (**a**–**e**) Before soil erosion; (**f**–**j**) After soil erosion; (**a**,**f**) *Artemisia capillaris* for 1 year since restoration; (**b**,**g**) *A. sacrorum* for 11 years; (**c**,**h**) *Bothriochloa ischaemum* for 15 years; (**d**,**i**) *Lespedeza davurica* for 25 years; (**e**,**j**) *Ziziphus jujube* for 40 years.

**Figure 3.** Short-term rainstorm frequency P–III curve in the study area.

#### *2.3. Vegetation Community Survey and Sample Collection*

In addition to experimental plots, three 2 × 2 m vegetation survey plots were established at each site. We sampled the vegetation by documenting the number of plants, estimating the percent ground cover, vegetation height, and humus layer thickness (Table 1). We harvested all aboveground plant tissue to measure aboveground vegetation biomass. We collected litter and humus layers for mass calculation. Root samples were collected using a root drill that was 9 cm in diameter with a barrel length of 10 cm from 0–20-, 20–40-, 40–60-, 60–80-, and 80–100-cm soil layers, based on previous reports of the root systems of vegetation on the Loess Plateau [1,28]. We used an iron box with an area of 0.2 × 0.15 m buckled into the soil to obtain undisturbed soil for measuring soil aggregates. A soil wreath knife with a volume of 100 cm<sup>3</sup> was used to obtain undisturbed soil for measuring the saturated hydraulic conductivity.


**Table 1.**Status of the experimental plots.

#### *2.4. Sample Testing*

All root samples were cleaned and scanned using a root scanner (EPSON, TWAIN PRO, Suwa City, Japan). We used the root-system analysis program WinRHIZO (QC., Quebec City, Canada) to analyze output from the root scanner. The program was used to estimate the root length, surface area, root tips, and diameter. Root biomass was measured by weighing the dried roots. Soil aggregate samples were air-dried naturally and then sifted into three grain classes to calculate the percentage, which were 0–0.25, 0.25–2, and >2 mm, respectively. Soils extracted from the drill core were used to measure the soil particle size using a laser particle size analyzer (Malvern, Mastersizer 2000, Birmingham, Britain). We only utilized soil median diameter d50 data in the grey correlation analysis. The saturated hydraulic conductivity of undisturbed soils collected with a wreath knife was measured by the constant head method.

#### *2.5. Data Analysis*

#### 2.5.1. Vegetation Community Index

We measured the root length density (RLD), root weight density (RWD), and root tip density (RTD) for all samples according to Equations (1)–(3):

$$RLD = \frac{L}{V\_{\text{s}}},\tag{1}$$

$$RWD = \frac{M}{V\_s} \,\tag{2}$$

$$RTD = \frac{N}{V\_s} \,\prime \tag{3}$$

where *L* is the sum of all root lengths per unit soil volume (mm); *M* is the dry weight of all roots per unit soil volume (mg); *N* is the sum of all the root tips per unit soil volume; and *V<sup>s</sup>* is the volume per unit of soil.

Next, we measured the soil saturated hydraulic conductivity as an indicator of soil permeability. The higher the saturated hydraulic conductivity, the higher the soil permeability. Increases in soil permeability can increase the infiltration of runoff and play a better role in soil and water conservation. This was calculated according to Equation (4):

$$K = \frac{10QL}{A \Delta H T'} \tag{4}$$

where *K* is the soil saturated hydraulic conductivity (mm/min); *Q* is the outflow (mL) in time *T*; *L* is the linear distance (cm) of the water flow path; A is the cross-sectional area (cm<sup>2</sup> ) through which the water flows; ∆*H* is the total head difference (cm) of the start and end sections of the percolation path; and *T* is the outflow time (min). We only used the soil saturated hydraulic conductivity as one of the factors of the correlation analysis and did not analyze the results in this study.

We used the Shannon–Wiener index (*H*), Margalef index (*R*), and the vegetation importance value (*Z*) as three ecological indicators that can reflect the aboveground structure of vegetation communities. The *Z* reflects the ecological status of a certain plant in a vegetation community and can play a normalization role for the study of more complex vegetation communities. Equations (5)–(7) were as follows:

$$H = -\sum\_{i=1}^{s} \frac{N\_i}{N} \{ \ln \frac{N\_i}{N} \} \tag{5}$$

$$R = \frac{S - 1}{\ln N} \tag{6}$$

$$Z = R\mathbf{D} + R\mathbf{F} + R\mathbf{C}\_r \tag{7}$$

where *N* is the sum of the number of plots; *N<sup>i</sup>* is the number of plants of the *i*-th species; *S* is the total number of species per plot; and *RD* = (density of a species/total density of all species) × 100%; *RF* = (frequency of a species/total frequency of all species) × 100%; *RC* = (cover of a species/total coverage of all species) × 100%.

#### 2.5.2. Calculation of Hydrodynamic Parameters

Hydrodynamic parameters included the Darcy–Weisbach (*f*), runoff shear (τ), runoff power (*P*), Reynolds number (*Re*), and Froude number (*Fr*) following Equations (8)–(15):

$$f = 8 \text{R} \cdot \text{J} \cdot \text{g}/v^2 \text{ }^{\circ}\text{C} \tag{8}$$

$$J = \left[L\cdot\sin\theta - \left(v^2/2\,\text{g}\right)\right]/L,\tag{9}$$

$$h = \frac{Q}{dv'}\tag{10}$$

$$Re = \frac{v\mathcal{R}}{\mathcal{Y}}\,'\,\tag{11}$$

$$\gamma = \frac{0.01775}{1 + 0.0337t + 0.000221t^2},\tag{12}$$

$$\mathbf{Fr} = \upsilon / \sqrt{\mathbf{gh}}.\tag{13}$$

$$
\pi = \rho \mathbb{R} \mathbb{J},\tag{14}
$$

$$
\omega = \tau v\_{\prime} \tag{15}
$$

where *R* is the hydraulic radius in m; *J* is the hydraulic gradient in m/m; g is the gravitational acceleration constant of 9.8 m/s 2 ; *v* is the runoff velocity in m/s; θ is the slope in degrees; *L* is the slope length in m; *h* is the depth of runoff in m; *Q* is the flow rate in m<sup>3</sup> /s; *d* is the runoff width in m; *R* is the hydraulic radius in m, which is approximately equal to the runoff depth, *h*; *t* is the water flow temperature in ◦C; and ρ is the water flow density in kg/m<sup>3</sup> .

#### 2.5.3. Grey Correlation Analysis

The gray correlation analysis formula is as follows, Equations (16) and (17):

$$\xi\_{0i} = \frac{\min\_{i} \min\_{k} \left| \mathbf{x}\_{0}^{\prime}(k) - \mathbf{x}\_{i}^{\prime}(k) \right| + \rho \max\_{i} \max\_{k} \left| \mathbf{x}\_{0}^{\prime}(k) - \mathbf{x}\_{i}^{\prime}(k) \right|}{\mathbf{x}\_{0}^{\prime}(k) - \mathbf{x}\_{i}^{\prime}(k) + \rho \max\_{i} \max\_{k} \left| \mathbf{x}\_{0}^{\prime}(k) - \mathbf{x}\_{i}^{\prime}(k) \right|} \tag{16}$$

$$\gamma(\mathbf{x}\_{0\prime}\mathbf{x}\_{i}) = \frac{1}{n} \sum\_{i=1}^{n} \xi\_{0i\prime} \tag{17}$$

where ξ0*<sup>i</sup>* is the correlation coefficient, γ(*x*0, *xi*) is the correlation degree, *x*<sup>0</sup> ′ (*k*) − *x<sup>i</sup>* ′ (*k*) is the difference sequence, max*i*max*<sup>k</sup> x*0 ′ (*k*) − *x<sup>i</sup>* ′ (*k*)  is the maximum difference, and min*i*min*<sup>k</sup> x*0 ′ (*k*) − *x<sup>i</sup>* ′ (*k*)  is the minimum difference.

#### **3. Results**

#### *3.1. Vegetation Succession Sequence and Structural Characteristics of Ground*/*Underground Parts*

We showed the distribution of *Z* for each species (Figure 4). We also summed the *Z* of each species in the five experimental plots. The rankings of the plant species in terms of *Z* across the experimental plots was: *Artemisia sacrorum* (356.72) > *Artemisia capillaris* (214.36) > *Bothriochloa ischaemum* (189.31) > *Lespedeza davurica* (177.97) > *Artemisia atrovirens* (123.67) > *Ziziphus jujuba* (100.74) > *Artemisia scoparia* (81.39) > *Tripolium vulgare* (58.32) > *Setaria viridis* (36.34) > *Taraxacum mongolicum* (10.27) > *Clerodendrum*

*mandarinorum* (9.66) > *Asparagus cochinchinensis* (5.36). The sum of *Z* of annual plants was 419.42, mainly *Artemisia capillaris*, *Artemisia scoparia*, and *Artemisia atrovirens*. Perennial herbaceous plants had a *Z* of 665.98, mainly *Bothriochloa ischaemum*, *Artemisia sacrorum*, and *Tripolium vulgare*. Small trees and semi-shrub plants had a *Z* of 278.71, which was primarily driven by *Lespedeza davurica* and *Ziziphus jujuba*. Therefore, it can be explained that perennial herbaceous plants are the main biological species for the natural restoration of vegetation in the area.

The *Z* across all species increased from 260.72 to 283.06 during the 40 years of succession that we analyzed. The *Z* of annual plants was 147.11, the perennial herbaceous plants was 93.12, and semi-shrubs and small trees was 20.49 in the 1-year vegetation community. Thus, annual plants were dominant in the 1-year vegetation community. Perennial herbaceous plants were dominant (153.23) in the 11- and 15-year vegetation communities by the same calculation method. Semi-shrubs and small trees were dominant (81.02) in the 25- and 40-year vegetation communities. It can be concluded that annual plants are the dominant species in the early stage of vegetation community succession, perennial grasses are the dominant species in the middle stage, and semi-shrubs and small trees are the dominant species in the later stage.

**Figure 4.** Distribution characteristics of the vegetation importance value.

The species density gradually increased from 2.5 to 4.5/m<sup>2</sup> from 1 to 40 years of succession (Table 2). Plant density increased first and then decreased and then stabilized. Plants density peaked at 135.5/m<sup>2</sup> at 15 years of succession, and then stabilized between 25 and 40 years (Table 2). Plant height increased from 25.25 to 159.32 cm/m<sup>2</sup> over the 40 years. Aboveground biomass increased from 28.83 to 753.33 g/m<sup>2</sup> . The amount of litter and humus increased from 24.50 to 605.00 g/m<sup>2</sup> .


**Table 2.** The ecological structure development of vegetation community succession from 1 to 40 years.

Roots are highly sensitive to the soil environment and occupy an important position in the succession of vegetation communities. The more closely the root system is integrated with the soil, the more obvious its effect on the soil's physical and chemical properties, and the stronger the soil erosion resistance [29]. RWD and root diameter gradually increased from 1 to 40 years of succession (Figure 5). At 40 years, the maximum RWD and root diameter were 10.80 mg/cm<sup>3</sup> and 0.65 mm, respectively. RLD and RTD increased from 1 to 15 years, then decreased from 15 to 25 years before stabilizing. The average RLD and RTD reached a maximum of 7.72 mm/cm<sup>3</sup> and 2.80/cm<sup>3</sup> in the 15-year successional community. According to the results above, it can be explained that slender and thin are the main root morphology of the vegetation community in early succession. Perennial plants, however, increased in dominance with increasing successional age. Specifically, we found an increase in the dominance of semi-shrubs and small trees after 25 years. At this time, the RWD and root diameter was bigger. The RWD, RLD, root diameter, and RTD decreased with soil depth at each successional stage, and there were significant differences in the root index between some soil layers (*p* < 0.05) (Figure 5).

**Figure 5.** Distribution characteristics of root morphology in vegetation communities. Note: The same letter (a,b,c) indicates that there was no significant difference in the indexes of root system among groups, the significant level P = 0.05.

−

#### *3.2. Dynamic Characteristics of Slope Runo*ff *in Di*ff*erent Vegetation Communities*

Soil erosion of vegetation communities on slopes is primarily affected by the hydrodynamic characteristics of runoff and by the composition of the underlying material. Runoff velocity ranged from 0.078 to 0.266 m/s across the whole experiment (Table 3). The runoff velocity range was 0.203 to 0.266 m/s under different flow conditions during the 1 to 15 years of succession. Runoff velocity ranged from 0.078 to 0.180 m·S <sup>−</sup><sup>1</sup> under different flow conditions during the 25 to 40 years of succession. In the later stage of succession, the runoff velocity was 52.37%, 38.77%, and 52.52% lower under the condition that the discharge rate was 4, 8, and 16 L/min, respectively. It can be concluded that vegetation community succession is a vegetation self-restoration process that can effectively reduce runoff velocity.

The Darcy–Weisbach (*f*) metric is used to indicate the resistance of the underpad to runoff. Generally, the larger the resistance coefficient, the more energy that is required for the water to overcome the resistance, and the smaller the sediment yield. The resistance coefficients ranged from 0.462 to 21.792 in the whole experiment (Table 3). The average resistance coefficient under different flow conditions was 0.827 in the 1 to 15 years of vegetation succession. The resistance coefficient was 11.223 during the 25 to 40 years of succession. The resistance coefficient of the slope runoff increased from 0.458 to 16.166 during the 25 to 40 years of vegetation succession.

Generally, the greater the flow shear stress, the greater the effective shear stress acting on the soil surface, and the greater the soil erosion intensity on the slope. The runoff shear stress increases with the increase of the discharge flow. When the scouring flow increased by 2.00 times, the shear stress increased by 1.63 times correspondingly in the whole experiment process. The runoff shear stress, however, also showed an increasing trend with the increase of the vegetation community succession years. The shear stress increased from 3.560 to 8.177 Pa during 1 to 40 years of vegetation community succession.

The runoff power can reflect the comprehensive influence of hydrodynamic characteristics on slope erosion. With the same trend of the shear stress, the runoff power of each vegetation community increased with the increase of the scouring flow (Table 3). This change in runoff power was due to the greater scouring force, faster runoff velocity, and greater shear stress, and the susceptibility to rill erosion. The maximum runoff power was 0.788 to 1.327 N/(m/s) in the 1-year vegetation community, and the minimum was 0.589 to 1.108 N/(m/s) in the 40-year vegetation community (Table 3).


**Table 3.**Characteristics of the dynamic parameters of slope runoffin different vegetation communities.

#### *3.3. Runo*ff *and Sediment Yield under Di*ff*erent Vegetation Communities*

The runoff volume and sediment yield on the slopes of the different vegetation communities were significantly different. Both the runoff volume and sediment yield decreased significantly with increasing successional age at a scouring flow of 4, 8, and 16 L/min (Table 4). Compared with the early successional community, the total runoff volume at 40 years decreased 3.52, 2.74, and 2.29 times, respectively, under the scouring flow of 4, 8, and 16 L·min−<sup>1</sup> . The total sediment yield decreased 16.83, 9.31 and 11.65 times on average. It can be seen from the multiple of reducing runoff and sediment that the effect of vegetation restoration on reducing runoff and sediment decreases with the increase of the erosion discharge. In addition, the total runoff volume and sediment yield under different scouring flows were averaged and the following results were calculated: During 1 to 40 years of vegetation succession, the runoff volume decreased by an average of 2.52 times and the sediment yield decreased by an average of 11.60 times. Therefore, it can be concluded that the contribution of Loess slope vegetation succession to sediment reduction during water erosion is much greater than that for runoff reduction.


**Table 4.** Total runoff volume and sediment yield under different vegetation communities.

The soil erosion rate of the 1-year vegetation community reached a maximum value of 1.35 g/(m<sup>2</sup> ·s) when the scouring flow rate was 16 L/min. The minimum soil erosion rate was 0.01 g/(m<sup>2</sup> ·s) in the 25-year vegetation community under a scouring flow rate of 4 L/min, followed by the 40-year vegetation community at 0.02 g/(m<sup>2</sup> ·s) (Figure 6). The soil erosion rate decreased by 25.36 times from 1 to 40 years of succession under the same scouring flow rate. This indicates that the successional status of vegetation in the Loess Plateau has a significant effect on reducing the erosion rate of runoff.

**Figure 6.** Average soil erosion and runoff rate under different vegetation communities.

#### *3.4. Grey Correlation Analysis between Hydrodynamic Parameters and Ecological Factors of Vegetation Communities*

The correlation between these hydrodynamic parameters and soil erosion rate can be illustrated (Table 5). The erosion rate was positively correlated with runoff velocity and power (*p* < 0.01). In contrast, the erosion rate was negatively correlated with resistance (*p* < 0.01), and was not correlated with shear force (*p* > 0.05).

The hydrodynamic parameters and the soil erosion rate showed a significant Pearson correlation, but the correlation coefficient does not indicate which hydrodynamic factors are most relevant to the erosion rate. Therefore, we identified the hydrodynamic factors most closely related to the erosion rate by the grey correlation analysis method: Runoff power (Table 6).


**Table 5.** Correlations between the soil erosion rate and hydrodynamic factors.

\* Indicates a significant correlation at *p* < 0.05; \*\* Indicates a significant correlation at *p* < 0.01.

**Table 6.** Correlation coefficient between the soil erosion rate and hydrodynamic factors and grey correlation degree.


In order to establish the relationship between the hydrodynamic parameters with the vegetation characteristics, we chose the hydrodynamic parameters as the characteristic indicators, and selected the vegetation community as the sequence index. The correlation degree of the total slope runoff

ξ

volume and ecological factors of the vegetation community were ranked as follows: Underground part (0.74) > aboveground part (0.69) = soil structure (0.69) (Table 7). From the grey correlation analysis of the slope runoff power, it can be concluded that the aboveground part (0.74) > soil texture (0.73) > underground part (0.69) (Table 8). The grey correlation analysis of the slope runoff velocity shows that the aboveground part (0.76) > soil texture (0.70) > underground part (0.68) (Table 9). The grey correlation analysis of the runoff resistance shows that soil texture (0.71) > underground part (0.70) > aboveground part (0.65) (Table 10).


**Table 7.** Grey correlation between the total slope runoff volume and vegetation communities.






#### **4. Discussion**

#### *4.1. E*ff*ects of Vegetation Community Restoration on Soil Structure and Erosion*

Vegetation construction is an important measure for soil erosion control on the Loess Plateau. Specifically, it plays an important role in controlling soil erosion and reducing sediment along the Yellow River. Vegetation restoration can effectively improve soil properties, affecting the relationship between the plant–soil interface through various factors, such as stems, leaves, roots, and root exudates [30]. Improving soil properties takes time, and there are differences in the degree and efficiency of soil improvement between different succession development directions and different vegetation types [31]. Therefore, the growth and decline of dominant species is the basic unit of the succession and development of a vegetation community. The composition of the community and the spatial distribution of individuals constitute the structural characteristics of the vegetation community. Vegetation community structure includes diversity, species composition, community floristic composition, and community age structure [32]. The natural restoration process of abandoned farmland vegetation on the Loess Plateau can be roughly divided into the rapid recovery period, primary succession period, advanced succession period, and stable period [33].

Vegetation communities have the ability to control soil erosion, in part because they increase soil anti-erodibility [34]. The soil erodibility is affected by vegetation as they can increase soil organic matter and soil aggregate stability. In this study, the shear stress increased from 3.560 to 8.177 Pa from 1 to 40 years of vegetation succession (Table 3). At the same time, the erosion rate decreased from 1.35 to 0.02 g/(m<sup>2</sup> ·s) (Figure 5). This demonstrates that communities at later stages of succession have better anti-erosion properties than younger communities. The effects of vegetation on soil erosion reduction can be divided into three parts, namely, the aboveground interception of vegetation, the fixation of underground roots, and the resistance of the soil interface [35]. The vegetation types and cover on the Loess Plateau have changed dramatically over the past 30 years [36], and the regional vegetation ecosystem has been significantly improved. These improvements have effectively alleviated the serious effects of soil erosion in this area. Jiao et al. showed that the average soil erosion intensity in the early stages of vegetation succession was between 3087.6 and 4408.4 t/km<sup>2</sup> /a, and the vegetation succession period was between 1245.2 and 1827.8 t/km<sup>2</sup> /a [37]. This is consistent with the finding that vegetation restoration can effectively reduce the rate of soil erosion (Figure 5). However, with the occurrence of extreme rainstorm events or the increase of scouring intensity, places with better vegetation coverage are more prone to small gravity erosion, such as landslides and collapses. [38]. It can be assumed that when the scouring flow rate is great enough, the vegetation will gradually lose its effectiveness in reducing flow and reducing sediment, and can even induce extreme soil erosion events, such as landslides and collapses. Previous studies have shown that when grasslands on slopes scour under 5 L/min of discharge, the effect of sediment reduction on vegetated slopes is significantly

greater than that on bare slopes while when the scouring flow increases to 8 L/min, the difference of sediment reduction between different grasslands is small, which indicates that the ability of vegetation to prevent and control runoff erosion on slopes is weakened with the increased discharge [39]. The total runoff volume decreased by an average of 3.52, 2.74, and 2.29 times at 4, 8, and 16 L/min of the scouring flow and the total sediment yield decreased by 16.83, 9.31, and 11.65 times on average in our data report (Figure 5). It was also just starting from the scouring flow of 8 L/min, and the difference in the control effect of vegetation communities on soil erosion became smaller. Therefore, the predecessors and our research can at least prove that when the scouring flow is large enough, the soil and water conservation of the vegetation gradually decreases. Whether it will aggravate slope erosion or induce gravity erosion remains to be further studied.

#### *4.2. E*ff*ects of Vegetation Community Restoration on the Hydrodynamics of Slope Runo*ff

Vegetation can effectively reduce water erosion on slopes. and some of the more important reasons are the hydrodynamic parameters affecting the runoff, which will change the runoff velocity, flow regime, and erosion energy of the slope, thus affecting the soil erosion process [25,40]. Vegetation can increase the critical conditions of slope erosion by increasing runoff resistance, and reducing runoff velocity and power to improve the critical conditions of slope erosion [41]. Their research shows that among many hydrodynamic factors, the runoff power on the slope is the most closely related to the average sediment transport rate of runoff, and the runoff power is the factor that can best reflect the soil erosion rate. It is easier to analyze and simulate the soil erosion process by using runoff kinetic energy and power theory [42,43]. Previous studies have shown that vegetation type, vegetation coverage, vegetation litter, humus, and roots are important factors affecting soil erosion on slopes. Vegetation communities regulate the runoff of slopes through the interaction of aboveground and underground parts [44]. On the one hand, the ecological structure of aboveground plants and underground roots in vegetation communities can increase runoff resistance and reduce runoff power [45]. On the other hand, the correct succession of vegetation communities can improve the soil properties, greatly enhance the soil anti-erodibility, so that the formation of rills cannot be fully developed, and the runoff hydraulic energy slope always changes little [46,47]. In view of this problem, we continue to discuss the correlation degree of slope runoff volume, velocity, resistance coefficient, power, and other factors that affect the soil erosion rate. The average slope runoff volume decreased by 2.85 times with the succession and development of the vegetation community (1–40 years) (Table 4). This contribution came primarily from the root system of the vegetation community, and the average correlation degree was 0.74 (Table 7). The RLD was the most effective indicator of the runoff volume for roots. Runoff power was reduced by 19.75% from 1 to 40 years of succession (Table 3). The contribution came primarily from aboveground vegetation, with an average correlation degree of 0.74 (Table 8). The litter and humus quality were the most effective factors affecting the aboveground part. Vegetation community succession reduced the slope runoff velocity by 47.89%, and its deceleration effect came primarily from aboveground tissue, with the average correlation degree reaching 0.76 (Table 9). The vegetation importance value (0.87), species density (0.84), and Shannon–Wiener index (0.76) played key roles in the aboveground part. Therefore, the vegetation community complicates the composition of species through the development of succession, which is more effective in controlling soil erosion than the vegetation coverage of a single species [44]. The succession and development of the vegetation community also increased the content of soil aggregates, improved soil structure, and thus increased the resistance of surface soil to slope runoff. The correlation degree of soil structure, a series of the sequence index, to runoff resistance reached 0.71 (Table 10). Among them, the content of soil micro-aggregate was a key factor for the increased runoff resistance, and the correlation degree reached 0.8.

#### *4.3. Implications for the Relationship between Vegetation Community Restoration and Slope Erosion*

Around 50 years ago, the main contradiction in China's Loess Plateau was food production and ecological restoration. Strong soil erosion conditions have not allowed humans to cultivate on slopes [48,49]. Therefore, the Chinese government has adopted a series of eco-economic compensation measures in the hopes of resolving this important problem. However, after seeing significant increases in vegetation cover, the global climate changed, and the Loess Plateau still experiences strong soil erosion under extreme precipitation events [50]. Therefore, the new contradiction points to the benefit and mechanism of vegetation in controlling soil erosion. There were differences in the development direction of vegetation community succession, which leads to differences in the underlying surface. In the future, the difference analysis and quantitative description of slope erosion patterns should be examined more carefully. On the one hand, quantitative discussion of the effects of aboveground parts, underground parts, and soil structure of vegetation communities on slope erosion is required, and on the other hand, it is necessary to deeply analyze the relationship between the anti-erodibility of plant communities and the improved soil anti-erodibility, and to establish an evaluation model. This information will help us to understand the mechanism of a vegetation community regulating runoff and sediment on a slope more comprehensively. At the same time, the coupling and feedback between the slope hydrodynamic process and ecological vegetation processes can also be clarified from a scientific point of view.

#### **5. Conclusions**

We comprehensively analyzed the successional development of vegetation community restoration on the Loess Plateau in China. We found that vegetation communities in later stages of succession have the ability to control runoff erosion on slopes. In the early stages of vegetation community succession, communities were dominated by wormwood plants. Perennial grasses were dominant in the middle stages of succession, and semi-shrub and small trees became dominant in the later stages. The plant aboveground and underground parts, such as plant density, the number of species, root length density, and root biomass, increased gradually with the development of vegetation succession. After 40 years of natural succession of vegetation communities, the slope runoff velocity decreased by 47.89%, the runoff resistance coefficient increased by approximately 35.30 times, the runoff power decreased by approximately 19.75%, the total runoff volume decreased by approximately 2.52 times, and the total sediment yield was reduced by approximately 11.60 times. On the one hand, the role of vegetation in preventing and controlling slope water erosion indicated that the vegetation important value, number of species, vegetation diversity (Shannon–Wiener index), and litter humus layer of the aboveground part significantly reduced the runoff velocity and power. The total amount of runoff was significantly reduced by the development of vegetation roots. On the other hand, vegetation communities improved the soil structure, in which runoff resistance was significantly increased by the restoration of soil micro-aggregate. Our results are important for vegetation restoration in the Loess erosion slope. They provide a scientific basis for the study of the influence of the vegetation community on the resistance and control of soil erosion and can be used to benefit an evaluation of water and soil conservation in the "Grain for Green Project" on the Chinese Loess Plateau.

**Author Contributions:** Conceptualization, E.C. and P.L.; methodology, Z.L. (Zhanbin Li) and Y.S.; formal analysis, Y.Z. and J.Z.; software, Z.L. (Zhan Liu) and Z.L. (Zhineng Li); E.C. wrote the manuscript and all authors contributed to improving the paper.

**Funding:** This study was supported by the National Key Research and Development Program of China (2016YFC0402407), the National Basic Research Program of China (No. 2016YFC0402404), the Natural Science Foundations of China (No. 51779204, 41701603), Shaanxi Province Innovation Talent Promotion Plan Project Technology Innovation Team (2018TD-037) and Shaanxi Provincial Technology Innovation Guidance Project (2017CGZH-HJ-06).

**Acknowledgments:** We thank the reviewers for their useful comments and suggestions. We would like to thank Murphy Stephen at Yale University for his assistance with the English language and grammatical.

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

### *Article* **Evaluation of Rainfall Erosivity Factor Estimation Using Machine and Deep Learning Models**

**Jimin Lee <sup>1</sup> , Seoro Lee 1 , Jiyeong Hong 1 , Dongjun Lee 1 , Joo Hyun Bae 2 , Jae E. Yang 3 , Jonggun Kim <sup>1</sup> and Kyoung Jae Lim 1, \***


**Abstract:** Rainfall erosivity factor (R-factor) is one of the Universal Soil Loss Equation (USLE) input parameters that account for impacts of rainfall intensity in estimating soil loss. Although many studies have calculated the R-factor using various empirical methods or the USLE method, these methods are time-consuming and require specialized knowledge for the user. The purpose of this study is to develop machine learning models to predict the R-factor faster and more accurately than the previous methods. For this, this study calculated R-factor using 1-min interval rainfall data for improved accuracy of the target value. First, the monthly R-factors were calculated using the USLE calculation method to identify the characteristics of monthly rainfall-runoff induced erosion. In turn, machine learning models were developed to predict the R-factor using the monthly R-factors calculated at 50 sites in Korea as target values. The machine learning algorithms used for this study were Decision Tree, K-Nearest Neighbors, Multilayer Perceptron, Random Forest, Gradient Boosting, eXtreme Gradient Boost, and Deep Neural Network. As a result of the validation with 20% randomly selected data, the Deep Neural Network (DNN), among seven models, showed the greatest prediction accuracy results. The DNN developed in this study was tested for six sites in Korea to demonstrate trained model performance with Nash–Sutcliffe Efficiency (*NSE*) and the coefficient of determination (*R* 2 ) of 0.87. This means that our findings show that DNN can be efficiently used to estimate monthly R-factor at the desired site with much less effort and time with total monthly precipitation, maximum daily precipitation, and maximum hourly precipitation data. It will be used not only to calculate soil erosion risk but also to establish soil conservation plans and identify areas at risk of soil disasters by calculating rainfall erosivity factors.

**Keywords:** rainfall erosivity factor; USLE R; machine learning; Deep Neural Network

#### **1. Introduction**

Climate change and global warming have been concerns for hydrologists and environmentalists [1–3]. Hydrologic change is expected to be more aggressive as a result of rising global temperature, that consequently results in a change in the current rainfall patterns [4]. Moreover, the Intergovernmental Panel on Climate Change (IPCC) [5] report showed that increasing rainfall events and rainfall intensity are expected to occur in the coming years [6]. Due to the frequent occurrence of greater intensity rainfall events, rainfall erosivity will increase, thus topsoil will become more vulnerable to soil erosion [7]. Soil erosion by extreme intensive rainfall is a significant issue from agricultural and environmental perspectives [8]. A decrease in soil fertility, the inflow of sediment into the river ecosystem, reduction of crop yields, etc., will occur due to soil erosion [9,10]. Therefore, effective best management practices should be implemented for better sustainable management

**Citation:** Lee, J.; Lee, S.; Hong, J.; Lee, D.; Bae, J.H.; Yang, J.E.; Kim, J.; Lim, K.J. Evaluation of Rainfall Erosivity Factor Estimation Using Machine and Deep Learning Models. *Water* **2021**, *13*, 382. https://doi.org/ 10.3390/w13030382

Academic Editor: Su-Chin Chen Received: 17 December 2020 Accepted: 28 January 2021 Published: 1 February 2021

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

**Copyright:** © 2021 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 (https:// creativecommons.org/licenses/by/ 4.0/).

of soil erosion. Furthermore, there is a need for a regional estimate of soil loss to proper decision-making related to appropriate control practice, since erosion occurs diversely over space and time [11].

During the last few decades, various empirical, physically based, and conceptual computer models [12] such as Soil and Water Assessment Tool (SWAT) [13], European Soil Erosion Model (EUROSEM) [14], Water Erosion Prediction Project (WEPP) [15], Sediment Assessment Tool for Effective Erosion Control (SATEEC) [16], Agricultural Non-Point Source Pollution Model (AGNPS) [17], Universal Soil Loss Equation (USLE) [18], Revised Universal Soil Loss Equation (RUSLE) [19] have been developed. Among the models, the USLE model is one of the most popular and widely used empirical erosion models to predict soil erosion because of its easy application and simple structures [20,21]. The USLE model [18] calculates the annual average amount of soil erosion by taking into account soil erosion factors, such as rainfall erosivity factors, soil erodibility factor, slope and length, crop and cover management factor, and conservation practice factor.

The Ministry of Environment of Korea has supported for use of USLE in planning and managing sustainable land management in Korea. To these ends, the USLE has been extensively used to predict soil erosion and evaluate various soil erosion best management practices (BMPs) in Korea. Various efforts have been made for the development of sitespecific USLE parameters over the years [22]. Yu et al. [23] suggested monthly soil loss prediction at Daecheong Dam basin in order to improve the limitation of annual soil loss prediction. They found that over 50% of the annual soil loss occurs during July and August. The rainfall erosivity factor (R-factor) is one of the factors to be parameterized in the evaluation of soil loss in the USLE. The R-factor values are affected by the distribution of rainfall amount and its intensity over time and space.

Rainfall erosivity has been widely investigated due to its impact on soil erosion studies worldwide. Rainfall data at intervals of less than 30 min are required to calculate USLE rainfall erosivity factors. The empirical equations related to R-factor based on rainfall data, such as daily, monthly, or yearly, available in various spatial and temporal extents, have been developed using numerous data [24,25].

Sholagberu et al. [26] proposed a regression equation based on annual precipitation because it is difficult to collect sub-hourly rainfall data to calculate maximum 30-min rainfall intensity. Risal et al. [27] proposed a regression equation that can calculate monthly rainfall erosivity factors from 10-min interval rainfall data. In addition, the Web ERosivity Module (WERM), web-based software that can calculate rainfall erosivity factor, was developed and made available at http://www.envsys.co.kr/~werm. In the study by Risal et al. [27] on the R-factor calculation for South Korea, 10-min interval rainfall data, which cannot give the exact estimate of maximum 30-min rainfall intensity, was used. The Korea Meteorological Administration (KMA) provides 1-min rainfall data for over 50 weather stations in Korea. Estimation of R-factor values for South Korea using a recent rainfall dataset is needed for present and future uses because climate change causes changes in precipitation pattern and intensity to some degrees. However, the process of calculation of R-factor from rainfall data is time-consuming, although the Web ERosivity Module (WERM) software can calculate rainfall erosivity factor [27]. Furthermore, the radar rainfall dataset can be used to calculate spatial USLE R raster values using Web Erosivity Model-Spatial (WERM-S) [28]. These days, Machine Learning/Deep Learning (ML/DL) has been suggested as an alternative to predict and simulate natural phenomena [29]. Thus, ML/DL has been used for the prediction of flow, water quality, and ecosystem services [30–34]. These studies have implied that ML/DL is an efficient and effective way to calculate R factor values using recent rainfall time-series data provided by the KMA.

The objective of this study is to develop machine learning models to predict the monthly R-factor values, which are comparable with those calculated by the USLE method. For this aim, we calculated R-factor using 1-min interval rainfall data to estimate the maximum 30-min rainfall intensity of the target values, which is monthly R factor values at 50 stations in S. Korea. In the previous study by Risal et al. [27], the R-factor values

for South Korea were calculated using 10-min interval rainfall data, which cannot give an exact estimate of maximum 30-min rainfall intensity. The procedure used in this study is shown in Figure 1.

**Figure 1.** Study procedures.

#### **2. Methods**

#### *2.1. Study Area*

— — Figure 2 shows the location of weather stations where 1-min rainfall data have been observed over the years in South Korea. The fifty points marked in circles are observational stations that provide data used for training and validation to create machine learning models predicting rainfall erosivity factors, while six stations marked in green on the right map—Chuncheon, Gangneung, Suwon, Jeonju, Busan, and Namhae—represent the stations for the final evaluation of the results predicted by machine learning models selected through validation. Thiessen network presented using red lines of the map on the right shows a range of the weather environment around the weather stations.

**Figure 2.** Weather stations in the study area.

or more as specified in the USLE users'

−1

−1 −1 −1 −1 −1

USLE users' manual in order to calculate the R

≤ →

−1 −1

→

−1

–

Σ

#### *2.2. Monthly Rainfall Erosivity Calculation*

Monthly rainfall erosivity (R-factor) was calculated for each of the 50 weather stations in South Korea from 2013 to 2019. It was calculated based on the equation given in the USLE users' manual in order to calculate the R-factor value [18]. According to Wischmeier and Smith [18], a rainfall interval of fewer than six hours is considered a single rainfall event. In addition, the least amount of rainfall that could cause soil loss is at least 12.7 mm or more as specified in the USLE users' manual [35].

However, if the rainfall is 6.25 mm during 15 min, it is defined as a rainfall event that can cause soil loss. The calculations for each rainfall event are as follows.

$$\text{IF I} \le 76 \text{ mm/hr} \to \text{e} = 0.119 + 0.0873 \log\_{10} \text{I} \tag{1}$$

$$\text{IF I} > 76 \text{ mm/hr} \to \text{e} = 0.283 \tag{2}$$

$$\mathbf{E} = \boldsymbol{\Sigma} \text{ ( $\mathbf{e} \times \mathbf{P}$ )}\tag{3}$$

$$\mathbf{R} = \mathbf{E} \times \mathbf{I}\_{\text{30max}} \tag{4}$$

where I (mm h−<sup>1</sup> ) is the intensity of rainfall, e (MJ mm ha−<sup>1</sup> ) is unit rainfall energy, P (mm) is the rainfall volume during a given time period, E (MJ ha−<sup>1</sup> ) is the total storm kinetic energy, I30max (mm h−<sup>1</sup> ) is the maximum 30-min intensity in the erosive event, and R (MJ mm ha−<sup>1</sup> h −1 ) is the rainfall erosivity factor. In this study, the monthly R-factor (MJ mm ha−<sup>1</sup> h <sup>−</sup><sup>1</sup> month−<sup>1</sup> ) was estimated by calculating monthly E and multiplying it by I30max. In addition, the monthly rainfall erosivity factor was calculated using Equations (1)–(4) [18] using the 1-min precipitation data provided on the Meteorological Data Open Portal site of the KMA (Korea Meteorological Administration).

#### *2.3. Machine Learning Models*

Machine learning can be largely divided into supervised learning, unsupervised learning, and reinforcement learning [36,37]. In this study, supervised learning algorithms were used. A total of seven methods (Table 1) were used to build models to estimate R-factor. Table 1 shows the information on machine learning models utilized in this study.


**Table 1.** Description of machine learning models.

Decision Tree, Random Forest, K-Nearest Neighbors, Gradient Boosting, and Multilayer Perceptron imported and used the related functions from the Scikit-learn module (Version: 0.21.3), while eXtreme gradient boost was taken from the XGboost library (License: Apache-2.0) and used the regression functions. Deep Neural Network is trained by taking Dense and Dropout functions from "Keras.models.Sequential" module of TensorFlow (Version: 2.0.0) and Keras (Version: 2.3.1) framework. In this study, the standardization method was used during the pre-process for raw data. Moreover, the "StandardScaler" function, a preprocessing library of Scikit-learn, was used.

#### 2.3.1. Decision Tree

The Decision Tree (DT) model uses hierarchical structures to find structural patterns in data for constructing decision-making rules to estimate both dependent and independent variables [38]. It first learns by continuing the yes/no question to reach a decision [39]. In this study, the DT model in the Scikit-learn supports only the pre-pruning. Entropy was based on classification and 2 for min\_samples\_split was given in Table 2.

**Table 2.** Critical hyperparameters in machine learning models.


A model hyperparameter is a value that is set directly by the user when modeling. Table 2 shows the hyperparameter settings of the regressors used in this study.

#### 2.3.2. Random Forest

Random Forest (RF) is a decision tree algorithm developed by Breiman [40] that applies the Bagging algorithm among the Classification and Registration Tree (CART) algorithm and the ensemble technique. RF creates multiple training data from a single dataset and performs multiple training. It generates several decision trees and improves predictability by integrating the decision trees [41]. Detailed tuning of the hyperparameter in RF is easier than an artificial neural network and support vector regression [42].

In this study, the hyperparameters in the RF are the following: 52 for n\_estimators, and 1 for min\_samples\_leaf.

#### 2.3.3. K-Nearest Neighbors

K-Nearest Neighbors (KNN) is a non-parametric method which can be used for regression and classification [43]. In this study, KNN was used for regression. KNN is an algorithm that finds the nearest "K" neighborhood from the new data in training data and uses the most frequent class of these neighbors as a predicted value [44]. In this study, the number of the nearest neighbors in KNN's hyperparameter was set as 3. The weights were calculated using a simple mean, and the distance was calculated by the Minkowski method [45].

#### 2.3.4. Gradient Boosting and eXtreme Gradient Boost

Gradient Boosting (GB) is an ensemble algorithm belonging to the boosting family that can perform classification and regression analysis [46,47]. In GB, the gradient reveals the weaknesses of the model that have been learned so far, whereas other machine learning models (e.g., DT and RF) focus on it to boost performance [48]. In other words, the advantage of gradient boosting is that the other loss functions can be used as much as possible. Therefore, the parameters that minimize the loss function that quantifies errors in the predictive model can found for better R-factor prediction. In this study, the hyperparameters in the GB are the following: 0.01 for learning\_rate, 4 for min\_samples\_split.

The eXtreme Gradient Boost (XGB) model is faster in training and classifying data than GB using parallel processing. It also has a regulatory function that prevents overfitting, which results in better predictive performance [49]. XGB is trained only by important features so that it calculates faster and performs better when compared to other algorithms [50,51]. The hyperparameters in the XGB are the following: gbtree for booster, and 10 for max\_depth.

#### 2.3.5. Multilayer Perceptron

Multilayer Perceptron (MLP) is a neural network that uses a back-propagation algorithm to learn weights [52]. MLP network consists of an input layer, a hidden layer, and an output layer (the R-factor). In this study, the hidden layer consisted of 50 nodes.

The hidden layers receive the signals from the nodes of the input layer and transform them into signals that are sent to all output nodes, transforming them into the last layer of outputs [53]. The output is used as input units in the subsequent layer. The connection between units in subsequent layers has a weight. MLP learns its weights by using the backpropagation algorithm [52].

#### 2.3.6. Deep Neural Network

Deep Neural Network (DNN) is a predictive model that uses multiple layers of computational nodes for extracting features of existing data and depending on patterns learn to predict the outcome of some future input data [54]. The invention of the new optimizers enables us to train a large number of hyperparameters more quickly. In addition, the regularization and dropout allow us to avoid overfitting. The package used to build DNN in this study was TensorFlow developed by Google. In this study, the DNN model structure consisted of 7 dense layers and 1 dropout (Figure 3). Additional details about DNN can be found in Hinton et al. [55].

**Figure 3.** Illustration of the proposed Deep Neural N **Figure 3.** Illustration of the proposed Deep Neural Network (DNN) for rainfall erosivity (R-factor) prediction.

**Count** 

–

'

#### *2.4. Input Data and Validation Method*

Input data were compiled as shown in Table 3 to develop machine learning models to assess the R factor. The corresponding month from Jan. to Dec. was altered to numerical values, because rainfall patterns and their intensity may vary every month over space. Total monthly precipitation, maximum daily precipitation, and maximum hourly precipitation were calculated monthly and selected as the independent variables. The data can be easily downloaded in the form of monthly and hourly data among the Automated Synoptic Obstruction System (ASOS) data from the Korea Meteorological Administration (KMA)'s weather data opening portal site and organized as input data.


**Table 3.** The input data for machine learning models.

The monthly R-factors data in the manner presented in the USLE for the 50 selected sites from 2013 to 2019 were designated as target values, and as the features are given in Table 4, month (1–12), total monthly precipitation, maximum daily precipitation, and maximum hourly precipitation were designated as the features. Among the data, 80% of randomly selected data were trained, the model was created, and then the remaining 20% of data were used for the validation of the trained model.

To assess the performance of each machine learning model, Nash–Sutcliffe efficiency (*NSE*), Root Mean Squared Errors (*RMSE*), the Mean Absolute Error (*MAE*), and coefficient of determination (*R 2* ) was used. Numerous studies indicated the appropriateness of these measures to assess the accuracy of hydrological models [56–58]. *NSE*, *RMSE*, *MAE*, and *R 2* for evaluation of the model accuracy can be calculated from Equations (5)–(8).

$$NSE = 1 - \frac{\sum \left(O\_t - M\_t\right)^2}{\sum \left(O\_t - \overline{O\_t}\right)^2} \tag{5}$$

$$RMSE = \sqrt{\frac{\sum \left(O\_t - M\_t\right)^2}{n}} \tag{6}$$

$$MAE = \frac{1}{n} \sum \left| M\_t - O\_t \right| \tag{7}$$

$$R^2 = \frac{\left[\sum (O\_l - \overline{O\_l}) \left(M\_l - \overline{M\_l}\right)\right]^2}{\sum \left(O\_l - \overline{O\_l}\right)^2 \sum \left(M\_l - \overline{M\_l}\right)^2} \tag{8}$$

where *O<sup>t</sup>* is the actual value of *t*, *O<sup>t</sup>* is the mean of the actual value, *M<sup>t</sup>* is the estimated value of *t, M<sup>t</sup>* is the mean of the estimated value, and *n* is the total number of data.


#### **Table 4.** Monthly R-factor calculated by the Universal Soil Loss Equation (USLE).

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

#### *3.1. USLE R-Factor*

For the selected 50 sites, monthly rainfall erosivity factors for each year from 2013 to 2019 were calculated, and the average monthly rainfall erosivity factors for seven years were obtained. Then, the seven-year average monthly rainfall erosivity, R-factor, was generated and shown in Table 1. Moreover, to give a comprehensive look at the degree of rainfall patterns by site, the average annual rainfall erosivity factor for each site is also presented in Table 4.

In this study, rainfall erosivity factor maps were generated to examine patterns of monthly R-factor calculated by USLE using rainfall data from 50 selected sites for evaluation. The R-factor distributions were mapped reflecting the geographical characteristics in South Korea (Figure 4). The high R-factor distribution in all regions during the summer months of July and August can be confirmed.

– **Figure 4.** Spatial distribution of monthly R-factor calculated by USLE, using rainfall data from 50 weather stations for the period 2013–2019.

− − − − − − The monthly R-factors for two months from July to August contribute more than 50% of the total average annual R factor value of Korea. The rainfall occurs mainly in the wet season and the likelihood of erosion is very high compared to the dry season. In such a case, using the average annual R-factor value can give a misleading amount of soil erosion. For these reasons, the monthly R-factor would be helpful in analyzing the impact of the rainfall on soil erosion rather than the average annual R-factor.

#### *3.2. Validation of Machine Learning Models*

Table 5 shows the prediction accuracy results (*NSE*, *RMSE*, *MAE*, *R 2* ) of seven machine learning models, by comparing the predicted R-factor. The results from the Deep Neural Network (DNN) showed the highest prediction accuracy with *NSE* 0.823, *RMSE* 398.623 MJ mm ha <sup>−</sup><sup>1</sup> h <sup>−</sup><sup>1</sup> month −1 , *MAE* 144.442 MJ mm ha <sup>−</sup><sup>1</sup> h <sup>−</sup><sup>1</sup> month −1 , and *R <sup>2</sup>* 0.840.


**Table 5.** Prediction accuracy results of seven machine learning models.

− − − − − − When comparing the results of DNN and the other machine learning models (Decision Tree, Random Forest, K-Nearest Neighbors, Multilayer Perceptron, Gradient Boosting, and eXtreme Gradient Boost), we can see that DNN provided more accurate prediction results over other machine learning algorithms. Moreover, the highest value of *NSE*, *RMSE*, *MAE*, and *R <sup>2</sup>* was found when the DNN was employed for the prediction R-factor values.

DNN had been proven for its good performance in a number of studies about the environment. In the study conducted by Liu et al. [59], the DNN showed better results, compared with results obtained by other machine learning algorithms, in predicting streamflow at Yangtze River. Nhu et al. [60] reported the DNN has the most impactful method in machine learning for the prediction of landslide susceptibility compared to other machine learning such as decision trees and logistic regression. In the study by Lee et al. [61], a DNN-based model showed good performance as a result of evaluating the heavy rain damage prediction compared to the recurrent neural network (RNN) in deep learning. Sit et al. [62] reported the DNN can be helpful in time-series forecasting for flood and support improving existing models. For these reasons, it has been shown that DNN performs better in various studies. − − − − −

In this study, the second best-predicted model is the K-Nearest Neighbors (KNN). The result from the KNN model showed prediction accuracy with *NSE* 0.817, *RMSE* 405.327 MJ mm ha <sup>−</sup><sup>1</sup> h <sup>−</sup><sup>1</sup> month −1 , *MAE* 149.923 MJ mm ha <sup>−</sup><sup>1</sup> h <sup>−</sup><sup>1</sup> month −1 , and *R <sup>2</sup>* 0.794 which indicates that the KNN is the most effective, aside from DNN, in predicting R-factor. According to Kim et al. [63], KNN has good performance results in predicting the influent flow rate and four water qualities like chemical oxygen demand (COD), suspended solids (SS), total nitrogen (TN), and total phosphorus (TP) at a wastewater treatment plant. **− − − − − −**

On the other hand, Decision Tree has prediction accuracy, with *NSE* 0.518, *RMSE* 657.672 MJ mm ha <sup>−</sup><sup>1</sup> h <sup>−</sup><sup>1</sup> month −1 , *MAE* 217.408, MJ mm ha <sup>−</sup><sup>1</sup> month- −1 , and *R <sup>2</sup>* 0.626. This means that Decision Tree is less predictable than other machine learning models (Random Forest, K-Nearest Neighbors, Multilayer Perceptron, Gradient Boosting, eXtreme Gradient Boost, and Deep Neural Network). Hong et al. [37] also reported Decision Tree has less accuracy for the prediction of dam inflow compared to other machine learning models (Decision tree, Multilayer perceptron, Random forest, Gradient boosting, Convolutional neural network, and Recurrent neural network-long short-term memory).

Figure 5 shows the scattering graphs of the R-factors predicted by the seven machine learning models and calculated by the USLE method. All machine learning results represent a rather distracting correlation with less agreement. However, in Figure 5h, the Deep Neural Network algorithms predicted USLE R values calculated using the method suggested by USLE users' manual with higher accuracy, *NSE* value of 0.823. gested by USLE users' manual with higher accuracy,

**Figure 5.** *Cont.*

**Figure 5.** Comparison of (**a**) Decision Tree, (**b**) Multilayer Perceptron, (**c**) K-Nearest Neighbor, (**d**) Random Forest, (**e**) Gradient Boosting, (**f**) eXtreme Gradient Boost, and (**g**) Deep Neural Network calculated R-factor with validation data, and (**h**) comparison of machine learning accuracy.

Among the data, 80% of randomly selected data were trained, the model was created, and then the remaining 20% of data were used for the validation of the trained model. To prevent overfitting, the K-fold cross-validation was implemented for *R <sup>2</sup>* as shown in Table 6. As a result of the five attempts of K-fold cross-validation, the DNN showed the best results with an average *R <sup>2</sup>* of 0.783.


**Table 6.** K-fold cross validation results of seven machine learning models.

Figure 6 shows the results of the prediction of the five machine learning models (i.e., Multilayer Perceptron, K-Nearest Neighbor, Random Forest, eXtreme Gradient Boost, and Deep Neural Network) at six sites for the testing of the selected models, as well as the time series comparison graph for 2013–2019 of the monthly R-factor values calculated on the USLE basis. At most sites, it showed that the time series trend fits well with a pattern similar to the USLE calculation value. In particular, looking at the distribution in Figure 6b Gangneung, the value of 9303 MJ mm ha <sup>−</sup><sup>1</sup> h <sup>−</sup><sup>1</sup> month −1 in October 2019, which represented the peak value of the rainfall erosivity factor, was generally well predicted by all machine learning models. Among the models, the result of the Random Forest model estimated a similar value with 8133 MJ mm ha <sup>−</sup><sup>1</sup> h <sup>−</sup><sup>1</sup> month −1 .

On the other hand, among six sites, the time series distribution values of the model prediction result in Busan showed a slightly different pattern from the USLE calculation R-factor. In particular, the result was overestimated as the values of 8241 MJ mm ha <sup>−</sup><sup>1</sup> h −1 month −1 in August 2014, and Multilayer Perceptron was almost twice overestimated at 16,725 MJ mm ha <sup>−</sup><sup>1</sup> h <sup>−</sup><sup>1</sup> month −1 . – −1 −1 −1

However, the Random Forest (8188 MJ mm ha <sup>−</sup><sup>1</sup> h <sup>−</sup><sup>1</sup> month −1 ) and eXtreme Gradient boost (8395 MJ mm ha <sup>−</sup><sup>1</sup> h <sup>−</sup><sup>1</sup> month −1 ) algorithms were predicting very similar values. Therefore, the machine learning results could be seen as good at predicting the peak value. −1 −1 −1

A comparison of the machine learning model accuracies of *NSE* and *R <sup>2</sup>* of the test (validation) results at the six sites is shown in Figure 7. All five models had a coefficient of determination of 0.69 or higher, and the simulated values of the USLE method calculation and machine learning models showed high accuracy prediction. However, compared to Deep Neural Network, the *NSE* results of the four models (Multilayer Perceptron, K-Nearest Neighbor, Random Forest, eXtreme Gradient Boost) were less than 0.58, and the Deep Neural Network model showed 0.87 in both *NSE* and *R 2* . Therefore, the monthly average value of the R-factor, predicted by the DNN would be a good candidate algorithm for USLE R factor prediction (Table 5 and Figure 7). −1 −1 −1 −1 −1 −1 −1 −1 −1 −1 −1 −1

**Figure 6.** *Cont.*

**Figure 6.** The comparisons of forecasting results of R-factor using machine learning in (**a**) Chuncheon, (**b**) Gangneung, (**c**) Suwon, (**d**) Jeonju, (**e**) Busan, and (**f**) Namhae.

**Figure 7.** Comparison of prediction accuracy results by machine learning models in test sites.

− − −

Table 7 shows average monthly rainfall erosivity factor values at the six sites for testing, Chuncheon, Gangneung, Suwon, Jeonju, Busan, and Namhae, along with the USLE calculation and Deep Neural Network (DNN) prediction.

− − −

− −

−


**Table 7.**Monthly R-factor calculated (C) by the previous method and predicted (M) by Deep Neural Network.

Among average annual vales, the results for Busan showed a good performance with the Deep Neural Network (DNN) resulting in the average annual value of the rainfall erosivity factor of 257 MJ mm ha−<sup>1</sup> h <sup>−</sup><sup>1</sup> year−<sup>1</sup> difference over the USLE calculation result. In the case of Chuncheon, DNN also showed a good performance with an average annual rainfall erosivity factor difference of 298 MJ mm ha−<sup>1</sup> h <sup>−</sup><sup>1</sup> year−<sup>1</sup> difference over the USLE calculation result. On the other hand, the USLE calculation results for Namhae showed an average annual value of the rainfall erosivity factor difference of 3361 MJ mm ha−<sup>1</sup> h <sup>−</sup><sup>1</sup> year−<sup>1</sup> greater than the DNN result.

This is because, in the case of Namhae, the rainfall tendency lasted for a long period in the dry season from February to June compared to the other testing sites like Chuncheon, Gangneung, Suwon, Jeonju, and Busan. Moreover, the monthly R-factor calculation of Namhae in dry seasons was two to four times more than other testing sites. In particular, the monthly R-factor for February in Namhae figure being about five times higher than the monthly R-factor in Busan. This means that if the single set of learning data has a huge deviation or variation from other sets, it may result in the uncertainty of the entire result data. Therefore, the monthly R-factor of Namhae in the dry season from is containing uncertainty. In the future study, when predicting the R-factor of the Namhae, DNN model analysis will be implemented in consideration of rainfall trends by supplement the historical rainfall data.

R-factor can be calculated by machine learning algorithms with high accuracy and time benefit. The spatio-temporal calculation of the rainfall erosivity factor using machine learning techniques can be utilized for the estimation of the soil erosion due to rainfall at the target value. The DNN will be incorporated into the WERM website in the near future after further validation.

#### **4. Conclusions**

The main objective of this study is to develop machine learning models to predict monthly R-factor values which are comparable with those calculated by the USLE method. For this, we calculated R-factor using 1-min interval rainfall data for improved accuracy of the target value. The machine learning and deep learning models used in this study were Decision Tree, K-Nearest Neighbors, Multilayer Perceptron, Random forest, Gradient boosting, eXtreme Gradient boost, and Deep Neural Network. All of the models except Decision Tress showed *NSE* and *R <sup>2</sup>* values of 0.7 or more, which means that most of the machine learning models showed high accuracy for predicting the R-factor. Among these, the Deep Neural Network (DNN) showed the best performance. As a result of the validation with 20% randomly selected data, DNN, among the seven models, showed the greatest prediction accuracy results with *NSE* 0.823, *RMSE* 398.623 MJ mm ha−<sup>1</sup> h −1 month−<sup>1</sup> , *MAE* 144.442 MJ mm ha−<sup>1</sup> h <sup>−</sup><sup>1</sup> month−<sup>1</sup> , and *R <sup>2</sup>* 0.840. Furthermore, the DNN developed in this study was tested for six sites (Chuncheon, Gangneung, Suwon, Jeonju, Busan, and Namhae) in S. Korea to demonstrate a trained model performance with *NSE* and *R <sup>2</sup>* of both 0.87. As a result of the comparative analysis of R-factor prediction through various models, the DNN was proven to be the best model for R-factor prediction in S. Korea with readily available rainfall data. The model accuracy and simplicity of machine learning and deep learning models insist that the models could replace traditional ways of calculating/estimating USLE R-factor values.

We found that the maximum 30 min intensity derived from 1-min interval rainfall data in this study is more accurate than that estimated from previous research. These methods can provide more accurate monthly, yearly, and event-based USLE R-factor for the entire period. Moreover, if the user has input data (month, the total amount of monthly precipitation, maximum daily precipitation, maximum hourly precipitation) as described in Table 3, the monthly R-factor can be easily calculated for the 50 specific stations in S. Korea by using the machine and deep learning models. Since the updated R-factor in this study reflected the recent rainfall data, which have high variability, it can improve the accuracy of the usage of the previous R-factor proposed by the Korean Ministry of

Environment [64] for future study. The results from this study can help the policymakers to update their guideline (Korean Ministry of Environment) [64] regarding the updated version of R-factors values for S. Korea.

It is expected that it will be used not only to calculate soil erosion risk but also to establish soil conservation plans and identify areas at risk of soil disasters by calculating rainfall erosivity factors at the desired temporal-spatial areas more easily and quickly.

However, this study evaluated the R-factor using machine learning models in S. Korean territory, under the monsoon region. Although deep learning models such as Deep Neural Network's applicability in S. Korea has been confirmed in this study, few studies have investigated and benchmarked the performances of a Deep Neural Network model-based USLE R-factor prediction trained. Therefore, future studies should be carried out for the diverse conditions of the other countries such as European countries, the United States, and African countries to broaden the applicability of machine learning technology in USLE R-factor (erosivity factor) analysis.

**Author Contributions:** Conceptualization and methodology: J.L., J.E.Y., J.K., and K.J.L.; formal analysis: J.H.B. and D.L.; data curation: J.L., J.H., and S.L.; writing—original draft preparation: J.L.; writing—review and editing: J.K.; supervision: K.J.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Ministry of Environment of Korea as The SS (Surface Soil conservation and management) projects [2019002820003].

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

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

#### **References**


#### *Article*

## **Impacts of the Degraded Alpine Swamp Meadow on Tensile Strength of Riverbank: A Case Study of the Upper Yellow River**

**Haili Zhu 1 , Peng Gao 2, \*, Zhiwei Li 3 , Jiangtao Fu 4 , Guorong Li 1 , Yabin Liu 1 , Xilai Li <sup>5</sup> and Xiasong Hu 1**


Received: 1 July 2020; Accepted: 18 August 2020; Published: 21 August 2020

**Abstract:** In the meandering riverbank of the Upper Yellow River (UYR), the native alpine swamp meadow (AS) has continuously degenerated into an alpine meadow (AM) due to climate change and intensified grazing. Its implication on river morphology is still not well known. This study examined this effect by in situ measurings of (1) physical properties of roots and their distribution in the soil-root mixture of the upper bank layer, and (2) the tensile strength in terms of excavating tests for triggering cantilever collapses of AS and AM riverbanks. The results showed that the root number in AS was significantly greater than that in AM, though the root distribution in both was similar. Also, the average tensile strength of individual roots in AS was 31,310 kPa, while that in AM was only 16,155 kPa. For the soil-root mixture, it decreased from 67.39 to 21.96 kPa. The weakened mechanical property was mainly ascribed to the lessened root number and the simpler root structure in the soil-root mixture of AM that reduces its ability to resist the external force. These findings confirmed that healthy AS can enhance bank stability and delay the development of tensile cracks in the riverbank of the meandering rivers in the UYR.

**Keywords:** alpine swamp meadow; alpine meadow; degradation of riparian vegetation; root distribution; tensile strength; tensile crack

#### **1. Introduction**

Banks of meandering rivers are often composed of silts and sand that have significantly higher compressive strength than tensile strength or cohesion [1–4]. There has been a consensus that riparian vegetation can reinforce bank strength [5–9] and stability through the interaction between soil and roots [10–16]. Mechanisms of riverbank failure are quite different from those found on hillslopes because steeper and shorter riverbanks tend to have a more variable profile and relatively small size of the failed block [17]. Accordingly, vegetation types and their root distributions throughout the bank profile play a critical role in resisting riverbank failure. In addition, the lower soil layer of the composite riverbank is subject to fluvial erosion, often resulting in a cantilever upper layer that includes the mixture of cohesive soil and vegetation. Once the gravity moment generated by the upper cantilever

layer exceeds its tensile moment, vertical cracks are developed and continuously extend through the upper layer, causing cantilever bank failure [10,18]. The ability of the soil-root mixture to resist external forces is a key factor in evaluating the bank stability.

Many studies have analyzed the tensile strength of fiber-reinforced soils and explored its relationship with soil physical indices using laboratory experiments [19–23]. These studies found that the benefits of natural or synthetic fiber reinforcement include (1) improved ductility in tension compared with pure earth blocks and (2) inhibition of tensile crack propagation [21–24]. Since the root system is intertwined in the soils, artificial soil blocks remolded in these laboratory experiments cannot reflect the actual root branching structure and the complex interaction between roots and soils. It is thus imperative to perform in situ tensile tests using naturally rooted soils for revealing the tensile properties of riparian vegetation.

In recent decades, a series of studies have been carried out for examining the shear strength of rooted soils using either in situ or indoor shear-test experiments. These studies helped better understand the enhancement of the vertically extending root system through soils by revealing that coarse roots tend to stabilize the soil-root mixture and fine roots may enhance its mechanical strength [25–27]. Other studies have confirmed that the number of roots passing through the potential shear plane, root system distribution and its strength [27,28], the initial water content of the soils [12,29], and the soil-root friction [30], contribute to soil reinforcement. However, there is still a lack of studies for quantifying the tensile strength of the rooted soil and the reinforced traction in the soil-root mixture.

The Upper Yellow River (UYR) watershed is located in the hinterland of the Qinghai-Tibet Plateau and includes numerous rivers, lakes, and wetlands. It is an important source of fresh water but a fragile ecosystem in the western region of China. When the UYR flows from the source into the open terrain with low-lying hills and valleys on the eastern edge of the Qinghai-Tibet Plateau, it forms a unique curved shape in plane form, commonly called the first bend of the UYR [31]. Meandering rivers are widely developed in this area, accompanied by the vegetation cover of the alpine swamp meadow and its degraded type, the alpine meadow. Because of global climate change and anthropogenic disturbances, the alpine swamp meadow on the eastern Qinghai-Tibet Plateau has been undergoing severe degradation [32–35]. The degradation succession of the alpine swamp meadow caused changes in the composition of the plant community. Dicotyledonous plants with a straight root system replaced sedges and gramineous plants with a dense clump root system, leading to changes in the underground biomass of the plant community, reduction of the spatial distribution of the root system, and a significant decrease of root activity and bulk density [35,36]. The biological degradation affects not only the root distribution of riparian vegetation but also the tensile and shear strength of the vegetated riverbanks. This means that degraded riparian vegetation could decrease bank stability and affect lateral evolution trends of meandering rivers in the UYR. Therefore, quantifying the influence of alpine swamp meadow degradation on the strength of riverbank soils in the UYR is an important issue that needs to be addressed. Most previous studies have examined the shear strength of the vegetated soils [8–13,16,25,26,37]. Little has been done on quantifying the effect of degraded riparian vegetation on the tensile strength of the soil-vegetation mixture.

In this study, we addressed this issue by focusing on the riverbanks of a meandering reach that has initially formed a cantilever arm under fluvial erosion. The vegetation community is featured by a healthy alpine swamp meadow (AS) and moderately degraded alpine meadow (AM). We measured root numbers and their vertical distributions in the soil-root mixture within a depth range of 0.3 m. Furthermore, by excavating the sandy layer below the top soil-vegetation layer of the bank in situ for artificially initiating cantilever bank collapse, we measured the tensile strength of individual roots using in situ pullout tests and subsequently calculated that for the soil-root mixture of the collapsed bank blocks. By comparing the differences of these results between AS and AM mixtures, we revealed the effects of meadow degradation on bank tensile strength and the role of roots in slowing bank crack development.

#### **2. Experiences at Field Scale and Procedure for Determining the Cantilanver Bank Collapse**

#### *2.1. Study Sites and Degradation of Alpine Swamp Meadow*

The Upper Yellow River (UYR) watershed is within the Qinghai-Tibet Plateau located in western China (Figure 1a). In its downstream reach, the main channel is joined by a relatively small tributary, the Lanmucuo River (Figure 1b). Our study sites are in the upstream reach of this tributary (34◦26′ N–35◦02′ N; 101◦29′ E–101◦35′ E) (Figure 1c). This reach has elevations ranging between 3400 and 4200 m a.s.l. with a mean channel gradient of 0.19%. The UYR region is subject to the alpine monsoon climate that features a long, cold dry season from October to early May, and a short, warm wet season from middle May to September. The annual mean precipitation is 560.5 mm [38], most of which is concentrated in the period from June to September, accounting for more than 83% of the total. The annual mean evaporation and temperature are 1278 mm and −0.16 ◦C, respectively [39]. Under this climate, the ground consists of seasonally frozen soils. ′ ′ ′ ′ −

**Figure 1.** (**a**) The geography of the Qinghai-Tibet Plateau and Upper Yellow River; (**b**) the location of the study area in the Lanmucuo river; (**c**) the specific locations of the two selected sites; (**d**) the spatially distributed AS and AM along the riparian zone of the Lanmucuo River; (**e**) illustration of *Blysmus sinocompressus*, the dominant species of AS; (**f**) demonstration of the AM that is rich in species composition (some degraded species are discernable).

Grassland is the major type of land cover in the region. It is dominated by alpine swamp meadow (AS), which provides important ecosystem services to the regional environment. Affected by global warming and intensified grazing, AS in the region has gradually degraded to alpine meadows (AM) and alpine steppe meadows (ASM) [40]. Spatially, AS and AM are often seen along riverbanks, the Lanmucuo River, while ASM is mostly distributed on the upper sloping parts of the piedmonts (Figure 1d). Because the area covered by ASM is far away from river banks, it is not a concern in this study. Based on a set of qualitative and semi-qualitative indicators used for pasture degradation classification [33,35], a survey for riparian vegetation along a 22 km reach of the upper Lanmucuo River was conducted in the 2017–2019 period [16,41]. During the survey, the coverage and number of species, dominant species, and underground biomass of AS and AM were determined within each sampling plot with the size of 1 × 1 m. The underground biomass was represented by the root mass within the soil layer that is 0.3 m in depth. The results (Table 1) showed that the values of all measured metrics were different with statistical significance between AS and AM.



\* standard deviation; AS: alpine swamp meadow; AM: alpine meadow; Different superscripts of a and b denote significant differences (*p* < 0.05) between different vegetation types.

In general, distributions of surficial plant species, coverage, and rooting depth and lateral root spread in underground root systems are determined by the prevailing physical habitats [42]. In particular, AS tends to appear in areas with topographic lows and occupied by seasonally saturated water (i.e., in swales), whereas AM is typically developed around the apex of river bends (Figure 1d). Vegetation communities of AS are mainly composed of cold-tolerant hygrophytes and hydromesophytes, which have a simple community structure. *Blysmus sinocompressus* is the dominant and healthy specie (Figure 1e and Table 1) and its coverage reaches nearly 98%. The remaining 2% is contributed from other species, such as *Ranunculus nephelogenes* and *Pedicularis longiflora*. The dominant species of AM are *Kobresia pygmaea*, *Elymus nutans*, and *Potentilla saundersiana*, accounting for 30% of the coverage (Table 1). Other herbaceous species, such as *Poa annua*, *Nardostachys jatamansi*, *Saxifraga montana*, *Aconitum tanguticum*, and degraded species, such as *Leontopodium pusillum*, *Oxytropis ochrocephala*, take up to 30% of the coverage. This means that AM only takes 60% of the habitat (Figure 1f). Compared with AS, the number of species and the mesophytes in AM are significantly high (Table 1). These surficial ecological differences between AS and AM must affect their underground properties and the associated mechanical characteristics, which were investigated in this study.

#### *2.2. Field Experiments and Measurements*

Field experiments were conducted at the two selected sites that are about 200 m apart from each other in the study area (Figure 1c). These sites have natural vegetation communities and are not disturbed by livestock and human activities. Therefore, they are representative of the general vegetation distribution in the study area. Site 1 is covered by AS, while site 2 is topped by AM. The height of the riverbank to the water level of the channel flow at site 1 is about 0.2–0.4 m lower than that at site 2, suggesting that the groundwater level at site 1 is higher than that at site 2. This ·difference is the main factor that determines the spatial distribution of AS and AM in the study area. The riverbank at both sites had been undercut by river flows, forming a cantilever arm with a width of about 0.2 m and thickness of about 0.3 m at site 1, while 0.3 and 0.35 m, respectively at site 2 (Figure 2). Both cantilever arms were stable with no tensile cracks developed on the top. They were underlined by a layer of silt and sand with a thickness of approximately 0.3 to 0.6 m, whose bottoms were gradually intergraded into fine gravels deposited at the toes of both banks (Figure 2). Field experiments were conducted during an October storm event in 2018, such that the entire bank profiles were exposed, facilitating subsequent bank excavation and measurements.

**Figure 2.** Riverbanks and their vertical profiles at (**a**) site 1, which is covered by AS, and (**b**) site 2, which is covered by AM (not to scale).

At each site, experiments were performed at the face of a bank profile stretched longitudinally for 8 m (Figure 3). An area of 1 m long on each side (i.e., E1 and E2 in Figure 3) was selected for extending the existing cantilever arm by excavating the lower part of the bank profile. The two areas were selected apart from each other by 3 m to ensure that the excavation of the first one would not affect the stability of the second. During each experiment, excavation was executed gradually, such that the development of tensile cracks on top of the cantilever arm could be observed. Their emergence indicated the status of arm stability. Excavation ended when the cantilever arm reached the threshold that triggered the failure of the cantilever arm. Three types of measurements were subsequently performed after bank collapse, determining tensile strengths of individual roots, measuring root diameter, number, and distribution, and sampling soil for both in situ and laboratory analyses.

The tensile strength of a single root for a given plant species (*T<sup>r</sup>* , MPa) may be determined by [2],

$$T\_r = 4F/(\pi d^2) \tag{1}$$

− − where *F* is the maximum pullout force of measured individual roots (N) and *d* is the diameter of the corresponding single root (10−<sup>3</sup> m). The value of *F* was measured using a HP-500 digital push-pull meter with a maximum load of 500 N and an indication error of ±0.5% (Leqing Aidebao Instruments Co. Ltd., Leqing, China). This measurement was taken at three plots within the experimental area; two (i.e., R1 and R3) were on the new face of the upper soil-root layer after cantilever failure and the other (i.e., R2) was at the face of the original upper bank (Figure 3). This design assured that the measured tensile strength accounted for spatial variability. At each plot, vegetation roots were separated from their surrounding soils by brushing soil particles away. Then, every single root was connected to the tension meter by a clamp. The value of *F* was recorded after applying a horizontal tensile force at a uniform speed until the root is broken or pulled out (Figure 4). The diameter of the

same root was also measured using a Vernier caliper with an accuracy of 2 <sup>×</sup> <sup>10</sup>−<sup>5</sup> m. If the broken position of the root was close to the clamp, then the recorded *F* value was biased and not counted. There were more than 30 roots measured in each plot.

**Figure 4.** Measuring the maximum pullout force of individual roots using a tension meter.

− Root diameters and numbers were measured within the front faces of the soil-root layer (0.3 m in thickness) in the collapsed bank block at R1 and R3, and the face of the same layer in the original bank at R2 (Figure 3). The surface of the face was washed by water to expose the root system and then divided into 30 small zones by laying a grid whose cells had the size of 0.10 × 0.10 m on the top. Within each of the 30 cells, root diameters were measured using the same Vernier caliper and the number of roots was recorded. All the roots that passed through the soil-root layer were measured. These values were classified into three groups based on the root diameter (i.e., <sup>&</sup>lt;0.5, 0.5–1.0, and <sup>&</sup>gt;1.0 <sup>×</sup> <sup>10</sup>−<sup>3</sup> m), and the depths with the soil-root layer (i.e., 0–0.10, 0.11–0.20, and 0.21–0.30 m). The root area ratio (RAR), defined as the ratio of the total area of all roots passing through the collapsed soil-root layer to the area of the plot (i.e., 1 × 0.3 m), was then calculated using the measured values. The measurements were repeated three times at each plot (i.e., R1, R2, and R3).

−− Soil samples were taken from two small plots in the soil-root layer around the depth of 0.10–0.20 m beneath the ground surface (i.e., S1 and S3 in Figure 3) and from two plots in the lower layer between 0.3 and 0.5 m from the ground surface (i.e., S2 and S4 in Figure 3) at both sites. These samples were carefully stored and transported back to the soil mechanics laboratory of Geological Engineering Department, Qinghai University for further analysis. The oven-drying method [43] was used to determine the soil moisture content, and the sieving method [43] was applied to determine the grain size distribution. Each measurement was repeated three times for accuracy. Additionally, a cutting

−

ring with a capacity volume of 6.0 <sup>×</sup> <sup>10</sup>−<sup>5</sup> <sup>m</sup><sup>3</sup> was used in situ to determine the soil bulk density at each plot (Figure 3).

#### *2.3. Calculation*

*σ*

−

In the process of cantilever bank collapse, the weight of the cantilever body is originally balanced by the tensile strength from the soil-root mixture of the body. Extension of the cantilever arm (by artificial excavation in this case or by fluvial erosion under the natural condition) will increase its weight. Once the tensile strength is insufficient to balance the weight, tensile cracks develop from the top of the cantilever body. At this time, both tensile stress and compressive stress exhibit a triangle distribution, and the center axis of the failing cantilever block is located in the stress center of the cantilever body below the crack [44] (Figure 5). Further extension of the cantilever arm catalyzes the development of the tensile crack until the failure of the cantilever body occurs. Under the critical condition of the failure, the external moment of the cantilever body is balanced by the resistance moment of the soil-root mixture, which may lead to [44],

$$\mathcal{W}b\_{\mathfrak{c}}/2 = \frac{l(d\_1 - d\_t)^2}{\mathfrak{Z}(1+a)^2}\sigma\_l + \frac{a^2l(d\_1 - d\_t)^2}{\mathfrak{Z}(1+a)^2}\sigma\_{\mathfrak{c}} \tag{2}$$

2 2

1 1

where *W* = ρ*gbcd*1*l* (N) is the weight of the cantilever body; ρ is the bulk density of the soil-root mixture (kg m−<sup>3</sup> ); *b<sup>c</sup>* is the critical width of the cantilever arm (m); *d*<sup>1</sup> is the thickness of the cantilever layer (m), *d<sup>t</sup>* is the depth of the crack from the top of the bank (m); *l* is the unit length of the cantilever layer (m), that is 1 m; *a* = σ*<sup>t</sup>* / σ*c*, σ*<sup>t</sup>* and σ*<sup>c</sup>* are tensile and compressive stress of the soil-root mixture (N m−<sup>2</sup> ), following Ajaz's results [45], *a* = 0.1. Substituting *W* = ρ*gbcd*1*l* into Equation (2), the tensile strength σ*<sup>t</sup>* of the soil-root mixture of the cantilever body may be expressed as 3(1 ) 3(1 ) <sup>1</sup> *ρ* − *σ σ σ σ* <sup>−</sup> 1

$$
\sigma\_t = 3(1+a)\rho g b\_c^2 d\_1 / 2(d\_1 - d\_t)^2 \tag{3}
$$

**Figure 5.** Stress analysis of the upper root-soil layer of a cantilever riverbank (modified from Figure 7 in Xia et al. [44], (not to scale)).

*σ* Using the measured *d*<sup>1</sup> and *d<sup>t</sup>* in our tested AS and AM bank blocks, we calculated their σ*<sup>t</sup>* values. The product of *T<sup>r</sup>* and RAR also reflects the tensile strength of a soil-root mixture. We compared these two types of calculations for the tensile strength.

− − −

#### **3. Results and Analysis**

#### *3.1. Characteristics of Soils and Root Distribution*

Our results demonstrated that the composition and physical properties of the two layers in the vertical profile of the riverbank are significantly different between the two sites (Table 2). Generally, the upper soil-root layer was composed of silt, while the lower layer consisted of silty sand with some poorly graded fine gravel. In the upper layer, the average bulk density and moisture content at the two sites were 1260 kg·m−<sup>3</sup> and 1560 kg·m−<sup>3</sup> ; 39.36% and 41.72%, respectively. They were about 180 kg·m−<sup>3</sup> and 405 kg·m−<sup>3</sup> ; 17.12% and 31.42% higher than those in the lower layer. These two soil properties in the soil-root layer of AS were 21.96% and 3.57% higher than those in that of AM. For the lower layer, the average bulk density and moisture content at site 1 were greater than those at site 2 by 118.67% and 3.27%, respectively. It is obvious that the physical parameters (i.e., bulk density and moisture content) of the soil-root mixture of AS and AM are disparate. These differences could have different influences on the tensile strength of the root system in the soil [27], which lend the support for our analyses of tensile strength for both single roots and the soil-root mixture in this study.


**Table 2.** The physical properties of the two vertical layers at both sites.

Although the distribution of roots over different classes of root diameters is highly variable for different species [46], it may still be characterized by vertical patterns of root number and diameters along the depth of a riverbank. In this study, AS is mainly composed of *Blysmus sinocompressus*, belonging to the *Cyperaceous* family. This species has a typical dense and fibrous root system that may be up to 0.8 m long, and its rhizomes are typically about 0.25–0.60 m long. These roots mix with the surrounding soil, forming a soil-root layer that extends from the bank surface to the depth of 0.30 m. Our results showed that the total root number within the experimental block of this layer, which was 1.0 m long and 0.3 m deep (Figure 3), was 4345, with 2505, 1160, and 680 in the depth ranges of 0–0.10, 0.11–0.20, and 0.21–0.30 m, respectively. Among these roots, those with the root diameter <0.5, 0.5–1.0, and <sup>&</sup>gt;1.0 <sup>×</sup> <sup>10</sup>−<sup>3</sup> m took about 66%, 25%, and 9%, respectively. AS is dominated by finer roots, which is evidenced by the fact that within each of the three depth ranges, they took 66%, 67%, and 76%, respectively (Figure 6a). Roots with medium diameters (0.5–1.0 <sup>×</sup> <sup>10</sup>−<sup>3</sup> m) in all three depth ranges took about 25%, while those with the diameters <sup>&</sup>gt;1.0 <sup>×</sup> <sup>10</sup>−<sup>3</sup> m only existed in the depths of 0–0.10 and 0.11–0.20 m, taking merely 9% and 8%, respectively. More fine roots (<0.5 <sup>×</sup> <sup>10</sup>−<sup>3</sup> m) developed in the shallow depth range rather than in the deep depth range, supported by their distributions of 56%, 26%, and 18% in the depth ranges of 0–0.10, 0.11–0.20, and 0.21–0.30 m, respectively (Figure 6a). The roots with the mean diameter (0.5–1.0 <sup>×</sup> <sup>10</sup>−<sup>3</sup> m) followed a similar vertical distribution, featured by 58%, 27%, and 15% in the three depth ranges, respectively. Roots with the diameters <sup>&</sup>gt;1.0 <sup>×</sup> <sup>10</sup>−<sup>3</sup> m only existed in the first two depth ranges with 71% and 29%, respectively.

−

−

*ρ* **<sup>−</sup>** *ω*

−

−

−

−

**Figure 6.** Root distributions in the three depth ranges of the experimental blocks for (**a**) AS and (**b**) AM. The number within each color of the three columns represents the percentage of roots in each sub-layer of the soil-root layer (i.e., the sum of the numbers in each color is 100%). The number outside of the columns represents the percentages of three root diameters in each sub-layer (i.e., the sum of these numbers along each column is 100%). The scaling of the y-axis is different between (**a**,**b**).

The AM is dominated by *Kobresia pygmaea, Elymus nutans,* and *Potentilla saundersiana*. Their root number took about 93% of the total in AM. The first two plants have dense and fibrous roots that are about 0.15–0.55 m long. Besides the three dominant plants, AM contains degraded species such as *Nardostachys chinensis*, *Oxytropis ochrocephala*, *Saxifraga montana*, and *Leontopodium pusillum*, whose roots are generally sparse and short. They have a typical tap root system with the length ranging between 0.04 and 0.13 m. The total number of roots in AS was only 2355 with 1420, 580, and 355 in the three depth ranges downward, respectively. Similarly to those in the AS, roots with the diameters <sup>&</sup>lt;0.5 <sup>×</sup> <sup>10</sup>−<sup>3</sup> <sup>m</sup> dominated in each of the three depth ranges, taking 71%, 71%, and 76%, respectively (Figure 6b). The remaining roots at each depth were those with the mean diameters (0.5–1.0 <sup>×</sup> <sup>10</sup>−<sup>3</sup> m). They took roughly the same percentage of the total roots within each depth range, which was 20%, 22%, and 24%, respectively. Again, the coarse roots (i.e., diameter <sup>&</sup>gt;1.0 <sup>×</sup> <sup>10</sup>−<sup>3</sup> m) only occupied the first two depth ranges, taking a small portion of the total roots in each (i.e., 9% and 7%, respectively). Along the vertical direction, roots with each diameter class demonstrated a similar distribution to their AS peers. For example, the fine roots (diameter <sup>&</sup>lt;0.5 <sup>×</sup> <sup>10</sup>−<sup>3</sup> m) took 60%, 24%, and 16% from the top to the bottom depth ranges (Figure 6b). A similar vertical distribution appeared for the roots with the mean diameter (i.e., diameters were between 0.1 and 1.0 <sup>×</sup> <sup>10</sup>−<sup>3</sup> m) with 57%, 26%, and 17% in the three depth ranges, respectively. Also, the coarse roots only stayed in the first two depth ranges. Most of them (75%) occupy in the 0–0.10 m depth range (Figure 6b).

Although roots in AS and AM showed similar spatial and diameter distributions, which led to their similar average root diameters (i.e., 0.46 <sup>±</sup> 0.31 <sup>×</sup>10−<sup>3</sup> m and 0.41 <sup>±</sup> 0.39 <sup>×</sup>10−<sup>3</sup> m, respectively) in the entire mixed layer of 0.30 m, their root numbers were greatly different, which can be proved by a two-sample difference test (*p* < 0.05). The root number of AS was 46% higher than that of AM. Consequently, the RAR of the AS experimental block was 0.225% on average, which was about twice that of the AM block (Table 2). These results showed that the root number, which plays an important role in resisting tensile force imposed to the experimental block, is reduced greatly when AS is degraded to AM.

#### *3.2. Tensile Trength of Individual Roots and Soil-Root Mixture*

The effect of roots on soil strength of the riverbank does not only depend on the root number, but also the tensile strength of a single root (*Tr*) (i.e., Equation (1)). The average *T<sup>r</sup>* of the dominant plant *Blysmus sinocompressus*, whose root diameters ranged between 0.20 and 1.76 <sup>×</sup> <sup>10</sup>−<sup>3</sup> m at site 1, was 31,310 kPa, and the maximum value can reach up to 128,000 kPa. The *T<sup>r</sup>* for roots with the diameters ranging between 0.24 and 1.88 <sup>×</sup> <sup>10</sup>−<sup>3</sup> m at site 2 (i.e., AM) was 16,160 kPa. The average

*T<sup>r</sup>* of the AS block was 48.4% higher than that of the AM block. There existed a strong relationship between the root diameter and *T<sup>r</sup>* (Figure 7), which indicated that *T<sup>r</sup>* decreases as the root diameter increases for both AS and AM. This finding is consistent with those reported in earlier studies [11,46–49]. This relationship may be described by a power function, whose exponent was different between the two sites (Figure 7). For the AS block, when the root diameter was less than 0.70 <sup>×</sup> <sup>10</sup>−<sup>3</sup> m, *<sup>T</sup><sup>r</sup>* decreased sharply as the root diameter increased, nearly following a linear relationship. Specifically, the value of *T<sup>r</sup>* decreased drastically from 127,930 to 30,810 kPa, as the root diameter only increased from 0.20 to 0.70 <sup>×</sup> <sup>10</sup>−<sup>3</sup> m. As the root diameter continuously increased from 0.71 to 1.76 <sup>×</sup> <sup>10</sup>−<sup>3</sup> m, *<sup>T</sup><sup>r</sup>* only decreased from 26,480 to 7050 kPa. The differences of *T<sup>r</sup>* were 97,120 kPa and 19,430 kPa before and after the threshold root diameter (i.e., 0.7 <sup>×</sup> <sup>10</sup>−<sup>3</sup> m). For the AM block, when the diameter was less than 0.4 <sup>×</sup> <sup>10</sup>−<sup>3</sup> m, *<sup>T</sup><sup>r</sup>* decreased greatly as the root diameter increased, and its value decreased from 73,380 to 23,590 kPa with the average of 39,050 kPa. However, the difference of *T<sup>r</sup>* was only 25,860 kPa when the root diameter increased from 0.4 to 1.88 <sup>×</sup> <sup>10</sup>−<sup>3</sup> m. This change was much less than that for root diameters less than 0.4 <sup>×</sup> <sup>10</sup>−<sup>3</sup> m (i.e., 49,790 kPa). Therefore, the root diameter of 0.7 and 0.4 <sup>×</sup> <sup>10</sup>−<sup>3</sup> <sup>m</sup> can be taken as a threshold of the AS and AM blocks, respectively. Values of *<sup>T</sup><sup>r</sup>* were more sensitive to the changes of root diameters less than the threshold while remaining less variable for coarse roots whose diameters are greater than the threshold.

**Figure 7.** Root tensile strength vs. root diameters for the alpine swamp meadow (AS) and alpine meadow (AM).

− Because there were about 93% and 90% fine roots with diameters less than 1 <sup>×</sup> <sup>10</sup>−<sup>3</sup> m in the mixed layer of AS and AM blocks, respectively, the difference in the root number should not be the main cause for the different *T<sup>r</sup>* values between AS and AM blocks. Rather, their difference represented the true mechanical discrepancy between the two. In other words, the tensile strength of individual roots in healthy meadow plants (i.e., AS) is generally higher than that in the degraded meadow plants (i.e., AM).

*σ* σ *σ* The tensile strength of the soil-root mixture should be relevant to the physical properties of the mixture. At site 1, the moisture content and volume of the two tested bank blocks were different. They were 40.92% and 0.175 m<sup>3</sup> for block 1, and 41.71% and 0.187 m<sup>3</sup> fo block 2, respectively (Tables 2 and 3). The tensile strength of the soil-root mixture (σ*t*) between the two tested sites was also slightly different. It was 66.86 kPa for Site 1 and 67.93 kPa for Site 2 (Table 3). At site 2, the moisture content and volume of the two tested experimental blocks were different by 1.06% and 0.004 m<sup>3</sup> (Tables 2 and 3), respectively. However, the variation of σ*<sup>t</sup>* between the two was still minor, merely 1.34 kPa (Table 3). This shows that for either AS or AM, repetitive measurements of σ*<sup>t</sup>* for different experimental blocks were consistent with each other, though the moisture content and volume of the experimental blocks might be slightly different.

σ

σ

*σ*


**Table 3.** Physical and mechanical properties of the soil-root mixture at the two sites.

These results suggest that variation of the moisture content and volume of the experimental blocks for the same type of plants has a negligible impact on σ*<sup>t</sup>* . Nonetheless, between the AS and AM, the mean tensile strength of the rooted soil was significantly different, which was 67.39 and 21.96 kPa, respectively (Table 3). The value of σ*<sup>t</sup>* for the AS at site 1 was about 3.07 times higher than that of the AM at site 2 (Table 3). The dramatic difference clearly reflected the differences in the mechanical characteristics of the root systems between the AS and AM blocks. It follows that the characteristics of root distribution and tensile strength of individual roots are the important factors that influence the tensile strength of the soil-root mixture.

The product of RAR and *T<sup>r</sup>* is often used to evaluate the contribution of the root system to the soil tensile strength [9,27,49–54]. For a soil-vegetation mixture, the tensile strength is often viewed as that from the soil and root system, and the latter is much greater than the former [9,51]. In this study, the product of RAR and *T<sup>r</sup>* for the AS and AM was 70.45 and 18.58 kPa, respectively. Their ratio was 3.79, which was greater than that of σ*<sup>t</sup>* for the AS and AM blocks (i.e., 3.07) (Table 3). This indicates that conditions of the root system in a soil-root mixture are critical for determining the mechanical characteristic of the mixture.

In the two excavation tests for AS, the depths of the developed cracks in the collapsed block were 0.035 and 0.037 m, respectively, with the average of 0.036 m (Table 3). However, for AM, these depths were 0.048 and 0.050 m, respectively, giving rise to the average of 0.049 m (Table 3). The crack developed in the AM collapsed block was deeper by about 27% than that in the AS block.

#### **4. Discussion**

#### *4.1. E*ff*ect of Degraded Riparian Vegetation on Tensile Strength of Individual Roots and the Soil-Root Mixture*

Continuous vegetation degradation has forced alpine meadow, dominated by members of the Cyperaceae, to transform into bare land and subsequently Heitutan on the Qinghai-Tibet Plateau in western China. Areas of severely degraded alpine meadow on the Qinghai-Tibet Plateau are referred to as Heitutan and are characterized by increased proportions of bare land, reduced edible herbage, and commensurate increases in the dominance of less palatable species [55]. Specifically, when AS degenerates into AM, the distribution of dominant plants reduced significantly, from an original coverage of 98% to only about 30%. This change of the surface vegetation communities has further led to changes in the plant's underground biomass. Our results showed that the root number of AM block was 46% less than that of the AS block within the depth of 0–0.30 m (Figure 6). For plants of AS, many roots were much longer than this depth, some of which may have a length of up to 0.80 m (Figure 2a). For the plants of AM, however, the lengths of most of their roots are less than 0.30 m (Figure 2b). Thus, transforming from AS to AM means that deep-rooted plants with dense root systems are gradually replaced by plants with short and sparse root systems [35,55–58]. This change indicates that the shallow-rooted plants of AM are more vulnerable to high evaporation, livestock trampling, and human activities in the study area, resulting in further degradation [35].

In addition to their different densities and lengths, the two types of plants also have different root structures. The dominant plant of AS (i.e., *Blysmus sinocompressus*) has developed rhizomes. Roots derived from them often grow laterally (Figure 8a). Moreover, some roots have a wave shape (Figure 8a). In AM, however, most roots of the dominant plant (i.e., *Kobresia pygmaea*) grow

vertically with no rhizome (Figure 8b). The other two dominant plants (i.e., *Elymus nutans* and *Potentilla saundersiana*) share a similar root structure. Compared with that of plant roots in AM, the structure of the dominant plant in AS effectively enhances the tensile strength of the soil-root mixture because (1) the laterally distributed roots may increase root density and (2) the wave-shaped roots have greater contact area with the surrounding soils and may resist higher external force by straightening its shape [59].

*σ*

**Figure 8.** (**a**) The laterally distributed and wave-shaped root system of *Blysmus sinocompressus* in AS; (**b**) the distribution of straight roots in *Kobresia pygmaea* of AM.

The above-mentioned characteristics of root diameters and branching, and the tortuosity of plant roots, mainly affected the mechanical behavior of individual roots [28]. Our study showed that the influence of vegetation degradation is on not only the root number and structure but also on the mechanical characteristics of roots. In the study area, effective roots (i.e., roots with diameters less than <sup>1</sup> <sup>×</sup> <sup>10</sup>−<sup>3</sup> m) account for about 93% in the layer of 0–0.30 m (Figure 6) below the surface, indicating that most root systems mainly play the role of reinforcement [60]. When the cantilever arm of the meandering riverbank is formed, the soil and root system is subjected to the external load from the weight of the arm, giving rise to the deformation of the arm. The root system can convert part of the external load into the tensile stress and dissipate it to the surrounding soil through the soil-root interface. In this way, the root system and the surrounding soil particles can work together to balance the load and enhance the soil-root tensile resistance [61].

The tensile strength of the root system is a critical factor that directly reflects the effect of rooted soil consolidation. In our study, the average tensile strength of the root for the dominant plants in the AS reached 31,310 kPa, which was 48% higher than that of the AM (16,155 kPa). Mattia et al. [47] and Li et al. [62] measured the root tensile strength of Gramineae plants of *Lygeum spartum*, *Stipa purpurea* and sedge plant of *Kobresia pygmaea*. They found that their tensile strength ranged from 36,260 to 45,670 kPa, close to that of the sedge plant, *Blysmus sinocompressus* (i.e., 31,310 kPa) in our study. The tensile strength of the root for degraded plants of *Potentilla bifurca*, *Ajania tenuifolia* and *Saussurea salsa* decreases significantly, from 5110 to 25,610 kPa [62]. In our study, this value for the degraded meadow plants (i.e., the AM) was within this range. Evidently, the mechanical characteristics of roots changed because of degeneration. The root tensile strength of healthy meadow plants (i.e., the AS) is much higher than that of the degraded plants (i.e., the AM).

The contribution of the root system to the tensile strength of the rooted soil may be appropriately evaluated using the product of RAR and *T<sup>r</sup>* [9,27,49,51,52]. To further analyze the reduction of the tensile strength in the soil-root mixture due to degradation, we compared our calculated values with those from earlier studies [12,47,62]. All of the products of RAR and root tensile strength of the four degraded herbaceous plants (i.e., 1, 2, 3, 5 in Figure 9) and the AM (i.e., 4) fell in the range of 2.04–38.4 kPa, which are lower than those of four healthy herbaceous plants (i.e., 7, 8, 9, 10 in Figure 9) and the AS (i.e., 6), which are within the range of 72.01–118.74 kPa. Clearly, there is a discrepancy of this product between healthy and degenerated plants, separated by the threshold of 50 kPa (Figure 9). −

**Figure 9.** The product of *T<sup>r</sup>* and RAR for healthy and degraded grassland plants. 1. *Saussurea salsa* [62]; 2. *Ajania tenuifolia* [62]; 3.*Potentilla bifurca* [62]; 4. *AM*; 5.*eontopodium nanum*[62]; 6. *AS*; 7. *Stipa purpurea* [62]; 8. *Lygeum spartum* [47]; 9. *Kobresia pygmaea* [62]; 10. *Helictotrichon filifolium* [46]. The hollow circles represent cited results, and the solid ones refer to the results in this study.

Because not all of the tensile strength of the roots is mobilized instantaneously at the moment of bank failure [9,50,51,63], the product of RAR and *T<sup>r</sup>* overestimates the true tensile strength (Table 3). As more and more fiber materials have been used in engineering practices, their ability to improve the tensile strength of the composite has been tested for natural fibers [24,64–67] and synthetic fibers [18,24,58,68]. It is generally believed that fiber materials can reinforce tensile strength, and the degree of reinforcement varies for different fiber materials. Tang et al. [3] proposed that the tensile strength of soil-fiber composite (σ*composite*) includes two parts, which are the tensile strength of natural soils (σ*soil*) and the increase of the tensile strength due to the fiber (∆σ*fiber*):

$$
\sigma\_{\text{composite}} = \sigma\_{\text{soil}} + \Delta \sigma\_{f\text{iber}} \tag{4}
$$

According to Zhu et al. [69], the tensile strength of unsaturated clays is 0.7~0.8 times that of the shear strength. In the current study, the ratio of tensile strength to the shear strength of soil without root is taken as 0.75, and the shear strength is 9.29 kPa [70]. Using these values and Equation (4), the reinforced tensile strengths of plant roots in AS and AM may be calculated as 60.43 and 15.0 kPa, respectively. The reinforcement of the healthy plant root system of AS is about 4 times higher than that of the AM. Also, they are less than the products of RAR and *T<sup>r</sup>* by 11.58 and 4.36 kPa, respectively.

The relationship between the value of ∆σ*fiber* for short (0.03 m) and long (0.05 m) fibers and different fiber contents [24] (Figure 10) was established by setting the root content of AS and AM as 0.48% and 0.18%, respectively [70]. It demonstrated a positive correlation between ∆σ*fiber* and the fiber content [1,20,21,67], and showed that the longer the fiber length, the greater the strengthening effect [3,24]. In Figure 10, the points representing AS and AM are in line with the established relationship for the fiber lengths of 0.05 and 0.03 m, respectively. This consistency indicates that both the root length and root content contribute to ∆σ*fiber*, and both variables are higher for plants of AS than those for plants of AM.

Furthermore, the tensile strength may also be affected by the initial water content and dry bulk density. In the fiber-reinforced soils, it is negatively related to water content and positively related to dry bulk density [3,24]. The average density of the natural soil-root mixture for AS is 18% higher than that of AM (Table 2). This may be attributed to the decreased root number and the dramatically increased pore volume of surface soil (i.e, soil in the 0–0.05 m layer) [58]. Higher bulk density should increase not only bonding forces between particles, which enhances the tensile strength, but also the interfacial contact area of the fiber-matrix structure, which improves the interfacial shear strength and associated friction [24]. In theory, the bonding forces between particles and friction of the soil-root

mixture for AS is stronger than that of AM. Because the difference in the moisture content of the upper layer of soils between AS and AM is merely 1.43% (Table 2), its influence on the tensile strength of the soil-root mixture is relatively small. In this study, the influence of plant degradation on the soil tensile strength is mainly analyzed from a mechanical point of view, and the influences of soil properties (e.g., bulk density, moisture content, porosity, and particle size) on tensile strength and the impacts of degradation on hydraulic properties need to be further studied. It has been well known that the root system can change the hydraulic characteristics of soil [71–74], because it occupies the pore spaces of soil, thus reducing the porosity and increasing the water-holding capacity of soil [75]. In this study, the root number of AM decreased by about 46% compared with that of AS in the depth of 0–30 cm on the upper part of the riverbank, suggesting that the proportion of pores in the soil mass occupied by roots was reduced, which is evidenced by the fact that the density of soil mass was relatively lower than that of AS. Therefore, the water holding capacity of soil mass should be lower in AM than that in AS, implying that its suction of soil mass should also be lower [76]. *σ σ* Δ*σ* Δ*σ*

According to Equation (4), the contributions of the root systems of AS and AM to the tensile strength of the soil-root mixture are 89.67% and 68.3%, respectively. This means that the influence of degradation on the tensile strength of the riverbank is mainly derived from the roots. In our study, when the AS degenerated to the AM, the rooted soil tensile strength decreased by 67.42% (Table 3). Li et al. [58] showed that the shear strength of the root-soil mixture decreases by 36.0% and 52.3% from severely degraded alpine meadows to moderately and slightly degraded alpine meadows, respectively. Thus, degradation of alpine swamp meadows obviously reduces the mechanical properties of the soils and weakens their ability to resist bank failure, wind and water erosion, and other external forces. Δ*σ* Δ*σ*

**Figure 10.** The relationship between fiber (root) content and reinforced tensile strength. The letter l in the diagram stands for the fiber length. The hollow circles and triangles are the data from Meriem et al. [24], while the solid circle and triangle refer to the reinforced tensile strength of AS and AM in this study, respectively.

#### *4.2. The Role of Root System in Preventing Development of Riverbank Cracks*

The development of cracks greatly destroys the integrity of the soil structure, weakens the mechanical properties of soils, reduces stability, increases permeability, intensifies evaporation, and increases soil erosion, resulting in a series of subsequent adverse effects on geotechnical engineering and the environment [77–80]. The tensile strength of the soil is an important mechanical parameter that controls the initiation and propagation of tensile cracks [3]. Therefore, enhancing the tensile strength of the soil and preventing the occurrence or slowing the expansion of cracks are critical for riverbank protection. Because plants in the AS have dense roots and strong single-root pullout resistance, they may enhance the tensile ductility of the rooted soils and inhibit the initial formation of tensile cracks. In our study, the expansion rate of tension cracks in cantilever arms due to weights

of the AS blocks and the duration of block failure are generally slower than those for the AM blocks. This indicates that types of riparian vegetation, fiber root number within the soils, and the root tensile strength may considerably constrain the number of cracks and their propagation rates [81].

When a tested block of cantilever arms is about to fail, a penetrating crack forms from the surface of AS or AM by following the path of least resistance (Figure 11). For the tested block of AS, the crack had a zigzag shape, and its length was about 1.24 m, which is longer than the crack length of AM by 0.14 m. The crack propagated with a preference angle until it was interrupted by roots. This angle seemed greater in the block with a higher root content. This property is consistent with the result of Meriem et al. [24]. At the depth of 0–0.3 m, the root number of AS that passing through the collapsed profile was 1.84 times that of the root number of AM (Figure 6). As a result, a crack in AS propagated in a way that avoided most roots. Given that the length of the crack was longer than that of AM, the time needed to penetrate through the soil-root mixture is longer. Nonetheless, the crack pattern is smoother on the surface of the AM block because of its lower number of roots.

**Figure 11.** The planform crack pattern at the surface of the block in (**a**) AS and (**b**) AM.

Considering the impact of the degraded alpine swamp meadow on tensile strength of riverbank in the UYR region, the grassland management should practice regional rotation grazing, rational distribution of herds, and balanced use of the riparian alpine swamp meadow. In addition, anthropogenic disturbances, in particular engineering construction, should be reduced as much as possible.

#### **5. Conclusions**

Under the influence of climate change and increased anthropogenic disturbances, the alpine swamp meadow (AS), which is the main type of vegetation cover in the meandering riverbank of the Upper Yellow River (UYR), is subject to severe degradation and typically has transformed into the alpine meadow (AM). However, little is known about the tensile strength of the roots in soils and the influence of riparian vegetation degradation on the tensile strength of riverbanks. To reflect the actual interaction between soils and the roots with the natural root branching structure for AS and AM and its effect on bank strength, we measured properties of root vertical distribution and number and performed in situ root pullout and artificial excavation tests. Our results led to the following conclusions. First, though spatial and size distributions of roots in the soil-root mixture of AS and AM were similar, AS was characterized by a higher number of roots than AM. Second, the tensile strength of individual roots decreased with the diameter of roots in both AS and AM. Yet, the former always had higher tensile strength than the latter for any given root diameter. Similarly, the tensile strength of the soil-root mixture in AS was about three times higher than that in AM. This difference is mainly caused by the fact that the lateral extended and wave-shaped root structure in AS is more effective in enhancing the resistance of the soil-root mixture to the external force than the simple vertically distributed root structure in AM. Third, the tensile crack developed in the collapsed block for AM was deeper than that for AS, indicating the reduced resistance to the external force as AS gradually degraded to AM. The impact of the roots for the degraded vegetation (i.e., AM) was also reflected by the relatively smooth and shorter crack route on the surface of the collapsed block for AM. These findings call for better ecological management for preventing AS from being degraded to AM.

Although grassland degradation is a worldwide problem [82,83], how such degradation affects the mechanical properties of riparian riverbanks in the alpine environment has not been fully understood. Our study provides firsthand evidence of weakened bank strength in the soil-root mixture due to degradation of vegetation and will serve as a benchmark for future investigation of riverbank strength, not only in the UYR, but also in other alpine regions in the world.

**Author Contributions:** Conceptualization, H.Z. and P.G.; Data curation, H.Z., G.L.; Formal analysis, H.Z., P.G.; Funding acquisition, H.Z., Z.L., G.L. and X.L.; Investigation, H.Z., Z.L., J.F., G.L., Y.L. and X.H.; Methodology, H.Z., P.G. and Z.L.; Project administration, H.Z.; Resources, X.L., G.L. and J.F.; Writing—original draft, H.Z.; Writing—review & editing, P.G. 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 No. 41302258, 41762023, 41662023, and 51979012), the Project of Qinghai Science & Technology Department (Grant No. 2017-ZJ-776), the Natural Science Foundation Innovation Team Project of Qinghai Provincial Department of Science and Technology (Grant No. 2020-ZJ-904).

**Acknowledgments:** The authors thank undergraduate students B.L., B.X., R.Z. and Y.G. for data collection in the field. Thanks are extended to the anonymous reviewers and Editor for useful comments that improved the manuscript.

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

#### **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/).

### *Article* **Static Liquefaction Capacity of Saturated Undisturbed Loess in South Jingyang Platform**

**Rui-Xin Yan 1,2 , Jian-Bing Peng 2, \*, Jin-Yuan Zhang 3, \* and Shao-kai Wang 2**


Received: 15 June 2020; Accepted: 11 August 2020; Published: 16 August 2020

**Abstract:** According to a previous geological investigation, high-speed and long-distance loess landslides in the South Jingyang platform in Shaanxi Province are closely related to the static liquefaction of loess. Considering the typical loess landslides in this area, isotropic consolidated undrained (ICU) triaxial tests and scanning electron microscopy analyses were conducted in this study. The main conclusions are as follows: (1) The stress-strain curves indicate strong strain softening under different confining pressures. The pore water pressure increases significantly and then remains at a high level; (2) The liquefaction potential index (LPI) shows an increasing trend followed by stabilization; the larger the LPI is, the smaller the state parameter (ψ) is. The steady-state points of the loess are in the instability region; however, the steady-state strength is not zero; (3) Based on the ICU test results, the average pore diameter decreases; the shape ratio remains essentially unchanged; and the fractal dimension and roundness show different trends. The proportions of the macropore and mesopore decrease; that of the small pore increases slightly; and that of the micropore increases significantly; (4) The compression deformation of the highly spaced pores causes rapid strain hardening. A rapid strain softening results from the pore throat blockage at the beginning of particle rearrangement and reorganization. A stable strain softening is related to the agglomeration blocking of the reconstructed pore throat in the gradually stable stage of particle rearrangement and reorganization.

**Keywords:** loess; ICU; static liquefaction; mechanical behavior; pore structure

#### **1. Introduction**

The phenomenon of static liquefaction was first discovered by Terzaghi and Peck [1] when they performed an experiment on saturated silty fine sand. The phenomenon was described as follows: saturated silty fine sand with uniform viscous liquid properties appears under a very small disturbance action; this was termed "spontaneous liquefaction". Subsequently, many scholars [2–4] have reported this phenomenon. They proposed that this phenomenon should be distinguished from the liquefaction of saturated sand under a dynamic load, and they studied it separately. Therefore, a more accurate conceptual description of "static liquefaction" was gradually formed. In other words, the stress-strain curve of the material element or the sample shows an obvious strain softening characteristic in the static loading process. This indicates that the deviator stress can only maintain a very low shear

strength after reaching the peak value. When the external load continues to act, the soil shows strong instability. With the further cognition of the phenomenon of static liquefaction, many scholars [5–7] began to understand the influence of static liquefaction on the loss of the soil resistance, and performed detailed analysis of engineering disasters resulting from static liquefaction. Several laboratory tests were also conducted, and their results have provided a clear understanding of the mechanism of static liquefaction.

Due to the characteristics of loess, such as being porous, having weak cementation, and water sensitivity, it is prone to cause slope instability due to static liquefaction under the saturated condition. According to analyses of the characteristics of the loess landslides in the South Jingyang platform, Shaanxi Province, the loess landslides are closely related to static liquefaction of the saturated loess. Relevant scholars have carried out systematic research on the static liquefaction typical characteristics of loess landslides in this region. Leng et al. [8] speculated that the movement mechanism of the landslide was related to the liquefaction of loess, considering the water accumulation at the toe of the slope and the unique properties of the slip soil after the occurrence of the landslide in this area. In the process of shearing, the loess liquefaction led to a sharp decrease in the shear strength of the displaced landslide materials. Through a field investigation of the Western Miaodian and Zhaitou landslides and using a digital elevation model, Peng et al. [9] found that the static liquefaction of the sliding surface was caused by the loss of shear resistance and bond strength. Xu [10] conducted scores of field investigations and measurements of the loess landslides in the South Jingyang platform. It was determined that the shear opening of the flow landslides was low, and most of the landslides were in the loess saturated zone on the edge of the plateau and were caused by a rise in the groundwater level. The mechanism for these kinds of loess landslides is related to the static liquefaction of the saturated loess at the bottom. In addition, the slope of the flow slip (about 12◦ ) is much lower than the internal friction angle of the saturated loess, which also preliminarily confirms loess liquefaction in the process of flowing. According to Li et al. [11], the pore water pressure produced by the undrained shear at the bottom of the sliding body in the loess tableland area is much higher than the pore water pressure inside the sliding body. Therefore, shearing occurred at the bottom. The high pore water pressure liquefied the bottom material, which resulted in a large difference between the driving force and resistance.

The phenomenon and mechanism of loess static liquefaction depends on the systematic experimental analysis. Therefore, relevant scholars have carried out laboratory experimental research and analyzed the mechanism of loess static liquefaction. For example, Ma et al. [12] determined that surface water infiltration led to water accumulation at the bottom of the loess layer through isotropic consolidated undrained (ICU) and ring shear tests. The loess in this layer was highly liquefied: a high pore water pressure was generated rapidly, and a deviatoric stress reached the maximum value under low deformation. Zhuang et al. [13] conducted ring shear tests on saturated loess in this region and determined that the pore water pressure and shear resistance immediately rose to the maximum value after shearing. Thereafter, the pore water pressure was maintained near the maximum value; however, half of the shear resistance was lost after a small shear displacement. The effective stress path showed that the saturated loess was prone to liquefaction under the undrained shear action. This indicates that the loess landslide had liquefied during the sliding process. Yan et al. [14] conducted ICU tests on the loess in this area and determined that it was easy to form a saturated softening zone on the top of the paleosol layer (relative aquiclude) at the toe of the slope. The relative sliding between the loess particles at this location can increase the pore water pressure sharply and reduce the effective stress. After reaching the steady state, the saturated loess was in a low-confining-pressure unstable state and was prone to plastic flow. Li and Jin [15] used the global digital systems (GDS) triaxial apparatus to conduct isotropic/anisotropic consolidated undrained (ICU/ACU) compression tests and constant-shear-drained (CQD) triaxial tests of the undisturbed loess in the Dongfeng landslide for the South Jingyang platform. The results showed that the stress-strain mode of the soil displayed a strong strain softening type, and the starting friction angle of the soil failure was far less than the steady-state

friction angle. When the stress path passed through the instability area, incomplete drainage shear shrinkage failure occurred, which is a typical static liquefaction phenomenon. Liu et al. [16] carried out undrained monotonic triaxial shear tests on saturated undisturbed loess samples taken from Heping Town, Yuzhong County, Lanzhou City. They found that loess samples with confining pressure of 100 kPa and 150 kPa exhibited complete static liquefaction characterized by that the pore water pressure reached the level of the initial confining pressure and the effective confining pressure decreased to zero. Wang et al. [17] also determined that the final pore pressures in the undrained triaxial compression tests were as great as 65% of the initial confining pressures. Besides, through effective stress path passing through flow liquefaction line (FLL), it can be concluded that the loess is susceptible to flow liquefaction failure. Pei et al. [18] conducted a series of consolidated undrained triaxial compression tests on undisturbed loess obtained from the source area of Shibeiyuan landslide. The tests revealed that although the confining pressures were different (≤200 kPa), the effective stress decreased with the accumulation of excess pore pressure and the steady state strength decreased continuously to zero. Xu et al. [19,20] carried out ICU tests on undisturbed loess retrieved from the backwall of a loess landslide in Heifangtai plateau, Lanzhou City, Gansu Province. The loess reached peak shear strength with an axial strain of <2%, accompanied by a sharp increase in pore pressure and followed by monotonic strain softening. The results showed that the saturated loess had strong liquefaction potential. Based on the study of the mechanical behavior of the static liquefaction of the saturated loess, the mechanism of the static liquefaction of loess has been clearly summarized as follows: O<sup>1</sup> the explosion of the pore water pressure; O<sup>2</sup> the effective stress drops sharply, even to zero; O<sup>3</sup> and the occurrence of flow plastic deformation. However, the internal mechanism is still unclear, such as: O<sup>1</sup> an unclear understanding of the mechanism of the explosion of the pore water pressure in the process of static liquefaction; O<sup>2</sup> incomplete cognition of the relationship between the occurrence of static liquefaction and the change in the internal structure of the saturated loess; and O<sup>3</sup> the essence of the sudden drop in the strength of the saturated loess due to static liquefaction, which can result in a sharp drop in the effective stress. However, the mechanism of the sudden drop in strength that is caused by the loess structure is not clear. Therefore, analyzing the internal mechanism of loess liquefaction and scientifically understanding the three typical characteristics of static liquefaction from multiple perspectives are the key work in the next stage.

In view of the above shortcomings, we carried out the ICU tests under different confining pressures for the saturated undisturbed loess. In addition, scanning electron microscopy was performed before and after the ICU tests. The stress-strain characteristics of the saturated undisturbed loess were studied. The variation trend of the pore water pressure and of the corresponding mechanical mechanism was also analyzed. Finally, the evolution process and internal mechanism of static liquefaction were determined.

#### **2. Development Characteristics of the Landslide and the Correlation Analysis with Static Liquefaction**

The South Jingyang platform is located in the South Bank of the Jinghe River, Jingyang County, Shaanxi Province. The platform is 30–90 m higher than the Jinghe River, and the altitude is 450–500 m. It is almost parallel to the flow direction of the Jinghe River and is distributed in a belt. The South Jingyang platform presents the east–west extension trend as a whole; the total extension is 28 km; and it mainly spans the Taiping, Jiangliu, and Gaozhuang administrative townships [21]. Under the influence of the buried fault of the Jinghe River, the quaternary loess was uplifted and was deposited on the south of the Jinghe River, and the loess platform was formed [22]. It has been affected by large-scale agricultural irrigation since 1980. As of April 2018, there have been 92 loess landslides on the South Jingyang platform. Among them, sliding had occurred multiple times for 17 landslides, which have caused serious casualties and property losses [23]. Based on the triggering factors and landslide movement characteristics, loess landslides in the study area can be classified into four types [22]: irrigation flow landslides, irrigation slide landslides, erosion slide landslides, and engineering-induced landslides (Figure 1). In previous studies, 92 landslides have been investigated, which include 35

irrigation flow landslides, 37 irrigation slide landslides, 15 erosion slide landslides, and 5 engineering induced landslides (artificial loading, cutting slope toe, engineering disturbance, etc.). The statistical results show that the number of irrigation flow landslides and the irrigation slide landslides in the study area is equivalent, which accounts for 38.04% and 40.22% respectively, which is far more than erosion slide landslides and engineering induced landslides (only 21.73% in total). The irrigation flow landslides largely slide out from the Middle Pleistocene (Q2) loess layer at the slope toe. Since the water content of the slope toe and the silt (or sand gravel) layer at the terrace is large, the soil near the sliding surface can produce high pore water pressures during movement. This makes the soil in the sliding zone liquefy. Even if the terrace slide bed is gentle, it still presents significant high-speed and long-distance characteristics. The plane shape of the sliding mass is primarily semicircular or circular, and the sliding distance is generally 200–300 m. The volume is mostly between hundreds of thousands of square meters to one million square meters, which indicates that these landslides are a medium-sized landslide. The shear opening of irrigation slide landslides is relatively high and is predominantly in the unsaturated loess layer. The irrigation slide landslides are associated with the reduction of suction (shear strength) of the loess matrix caused by the infiltration of surface water such as farmland irrigation. From the perspective of landslide evolution, sliding occurs mainly in the second and third stages of the landslide, while some of the sliding, which results from the deep buried groundwater level, occurs in the first stage.

**Figure 1.** The loess landslides distribution in the South Jingyang platform.

The geological structure of the loess slope exhibits different types of structural planes, including joint fissures and weak planes. In the South Jingyang platform, the joint fissures are mainly the fissures on the edge of the platform, and the weak planes are mainly the loess-paleosol interface. According to a field geological structure survey, fissures within 7 m from the edge of the platform are intensively developed; most of these are unloading fissures. The fissures appear along the top of the platform or are near the steep slope, and they continue to extend and penetrate under the unloading effect. However, due to the cutting action of the steep cliff, the extended length of these fissures is often small. Fissures more than 15 m away from the edge of the platform are primarily collapsible fractures and are predominantly at the top of the platform. Under the action of uneven collapsibility, they extend

deeply [12,14]. Although the number of collapsible fissures is less than unloading fissures, the plane extension length and joint opening are larger. Collapsible fissures are the dominant channel for the surface water to infiltrate the loess slope [12,24,25], which can cause the surface water to quickly invade the inside of the landslide body and to accelerate the rise of the groundwater level [21]. The South Jingyang platform presents typical interbedded sequence characteristics of the loess-paleosol [14], and the loess-paleosol contact interface is significant (Figure 2). However, because it is influenced by the tension of the fissures on the edge of the platform, the S1–S5 paleosol layer often breaks, which can easily cause water infiltration. Paleosol has a lower infiltration rate and water-holding capacity than loess does [26] and can be regarded as a natural aquiclude in comparison with loess [14]. The accumulation of surface water near the interface of a certain loess-paleosol is enhanced by the infiltration-promoting effect of the fissures on the edge of the platform. This causes the bottom of the Lishi loess to be saturated, which produces the saturated softening layer. Therefore, liquefaction may occur at the loess-paleosol interface, after which the saturated softening layer becomes the bottom sliding surface of the landslide [27].

**Figure 2.** The profile of loess-paleosol contact interface in the South Jingyang platform.

At about 10:00 AM on 26 July 2014, a large-scale loess landslide occurred in Hetan village, Western Miaodian. This can be regarded as the first large-scale landslide because of the long-time interval from the Western Miaodian landslide in 1987. In fact, four landslides occurred on 26–28 May 2015 and 7–8 August 2015. However, with the increase in the sliding sequence, the volume in the landslide decreased gradually [12]. This indicates that the landslide that occurred on July 26, 2014 had a triggering effect. The main sliding direction of the landslide was NE61◦ , which is slightly different from that of the 1987 landslide (NE48◦ ). A giant sliding body with a sliding distance of 278 m and a total volume of about 50 <sup>×</sup> 104 m<sup>3</sup> was formed (Figure 3). The posterior wall of the landslide was about 40 m high and 222 m wide, and it was characterized by a typical concave ring chair. The top 10–15 m was vertical and was controlled by the Late Pleistocene (Q3) and Q2 loess vertical joints. There were five paleosol layers that were exposed, and the landslide started from the vicinity of the fifth paleosol layer; moreover, the water retention effect of the paleosol layer was significant. In addition, the sliding surface of the posterior wall was smooth, and the water content was high.

**Figure 3.** Development characteristics of the "7.26" high-speed and long-distance landslide in Hetan village, Western Miaodian.

#### **3. Experimental Samples and Process**

≤ Based on the above investigation and analysis, the occurrence of landslides in this region is closely related to static liquefaction. In order to further analyze the evolution process and the typical characteristics of the static liquefaction of loess in this region, a static liquefaction test was conducted. In addition, because the thickness of the Holocene (Q4) loess (generally ≤ 5 m), Late Pleistocene (Q3) Malan loess (generally < 10 m), and other new loess were not large, these were mainly located at the top of the loess slope. Their failure mode was mainly a tensile fracture, and shear liquefaction was unlikely to occur. Although the Middle Pleistocene (Q2) Lishi loess (10–70 m) was the main component of the loess slope in the study area, shear failure and sliding were the main failure modes in this layer. Therefore, the Q2 loess was selected for the ICU tests.

#### *3.1. Experimental Samples*

The undisturbed loess samples were acquired from the Zhaitou (ZT) and Shutangwang (STW) landslide. Among them, the ZT loess samples were obtained from the adit of the posterior wall of the landslide. The sampling location is about 21 m away from the top of the platform, and the corresponding stratum number is L2. The ZT loess samples were yellow-brown, slightly wet, and have slightly dense uniform soil. The samples were mainly composed of silt particles. In addition, needle pores, a small amount of snail shell fragments, and concretions were found in the samples. The STW loess samples were obtained from the surface of the posterior wall of the landslide. The sampling location is about 17 m away from the top of the platform, and the corresponding stratum number is L2. The STW loess samples were light yellow-brown and dry, and its other characteristics were like those of the ZT loess samples. The ZT and STW loess consisted of Q2 Lishi loess.

The physical properties of the samples were determined by referring to the relevant provisions of the standard for the soil test method (GB/T 50123-1999) [28]. The water content was measured using the drying method; the density was measured using the round knife method; and the dry density value was obtained via conversion. The specific gravity was measured using the pycnometer method, and the liquid limit and plastic limit water content were determined using the combined liquid plastic limit method. The basic physical properties of the two loesses are shown in Table 1.


**Table 1.** Basic physical properties of the Zhaitou (ZT) and Shutangwang (STW) loess.

It can be observed from Table 1 that the physical properties of the ZT and STW Q2 loess were close except for the water content. ZT loess samples were procured from the adit, while the STW loess samples were procured from the surface of the slope. Due to the water evaporation at the surface loess, the measured water content of the STW loess samples was relatively low. The liquid index of the two loess was less than zero and the plasticity index was low. This indicates that their plasticity in the natural state was low.

In addition, a laser particle size analyzer was used to analyze the grading of two Q2 loess samples, and a semi quantitative analysis of their mineral compositions was performed using an X-ray diffractometer. The selected diffraction angle was 0.5–30◦ . The specific results are presented in Table 2.

**Table 2.** The particle and mineral composition of the ZT and STW loess.


As shown in Table 2, the particle composition of the ZT and STW Q2 loess were quite similar, with the proportion of the silt being the most, clay being the second, and fine sand being the least. Meanwhile, the corresponding mineral composition and content were highly similar. The highest content was quartz and carbonate, which constituted the main components of the cement and the skeleton of the loess. Clay minerals mainly consisted of illite, chlorite, and a relatively small amount of kaolinite. The hydrophilic characteristics of illite and chlorite were obvious, which makes it easy for the clay minerals to agglomerate and to form tuberculosis under water saturation conditions. The above mineral composition objectively revealed the typical structural characteristics of the loess in the area.

#### *3.2. Experimental Process*

In view of the typical characteristics of the static liquefaction phenomenon of the South Jingyang platform landslide: the pore water of the saturated loess caused by the groundwater immersion cannot be discharged in time during the shear process, which caused the pore water pressure to continue rising and to maintain a high level. Therefore, to accurately reflect the stress state of the saturated loess during the sliding of the loess slope, isotropic consolidated undrained (ICU) triaxial tests were conducted on ZT and STW loess. In addition, to further explore the internal mechanism of the static liquefaction of loess in this area and finding the internal relationship between the static liquefaction and the microstructure characteristics directly, it is necessary to observe and analyze the microstructure of the ZT and STW loess before and after the ICU tests. The specific experimental process is as follows:

(1) The undisturbed loess samples were molded into <sup>Φ</sup>50 mm <sup>×</sup> 100 mm cylinders according to the requirements of the standard for the soil test method (GB/T 50123-1999) [28]; and these cylinders were placed in the load cell of the Wykeham Farrance (WF, Wykeham Farrance, Milan, Italy) stress path triaxial apparatus (Figure 4). After connecting the pressure sensor, pore water pressure sensor, cell pressure controller, and the back pressure controller with the corresponding port of the load cell, the WF stress path triaxial apparatus was used to carry out the back pressure saturation treatment on the undisturbed loess. O<sup>1</sup> A 50 kPa cell pressure was applied to the

samples, and the back pressure was set at 10 kPa lower than the cell pressure to prevent structural damage of the samples caused by the pore water pressure, which is higher than the confining pressure. After the back pressure inflow was stable, the back pressure was kept constant for 6–8 h to make the water gradually infiltrate the samples. O<sup>2</sup> Under the premise of ensuring that the cell pressure was 10 kPa higher than the back pressure, the cell pressure and back pressure were applied to the samples step by step. Each stage was increased by 50 kPa and it remained for 6–8 h before conducting the next stage. O<sup>3</sup> When the B value was higher than 0.95, the samples were considered to be saturated. Multiple tests showed that the undisturbed loess samples generally required an applied back pressure that ranged from 190 to 240 kPa, and in 1–2 days to ensure that the soil sample reached the saturation state.


**Figure 4.** Wykeham Farrance (WF) stress path triaxial apparatus.

#### **4. Results**

#### *4.1. Macroscopic Characteristics Analysis of the Shear Failure and Static Liquefaction*

At first, the stress-strain characteristics under the different confining pressures were studied to obtain the liquefaction potential index (LPI) evolution process of the undisturbed loess. Moreover, the variation law of the pore water pressure and the influence of the confining pressures on the pore water pressure were analyzed. Finally, the evolution process of the static liquefaction was determined.

#### 4.1.1. Stress-Strain Characteristics Analysis

Based on the experimental data from the ICU tests, the stress-strain curves of the ZT and STW undisturbed loess under the confining pressures of 150 kPa, 250 kPa, 350 kPa, and 450 kPa were plotted (Figure 5).

**Figure 5.** Stress-strain curves of the ZT and STW loess. (**a**) ZT. (**b**) STW.

It can be observed from Figure 5 that the stress-strain curves of the ZT and STW loess under the different confining pressures had a similar change in the trends, and all of them present obvious strain softening characteristics [12]. The stress-strain curves can be divided into three stages.


*ε*

*ε*

of the loess structure tends to be stable, and its initial structure is completely destroyed under the action of the shear stress. In addition, the new orientation arrangement has statistically reached a stable state [31]. Therefore, in the ICU tests of the saturated undisturbed loess, when the axial strain reaches 10–15%, the loess is considered to undergo steady-state deformation. When the loess reaches a steady state, its shear resistance is known as the steady-state strength [32]. Compared with its peak strength, the steady-state strength is about 0.4–0.6. Therefore, there is the possibility of static liquefaction. The above structural changes can also find the corresponding intuitive evidence in the later microstructure morphology analysis.

*ε*

*ε*

#### 4.1.2. Variation Analysis of the Pore Water Pressure

Based on the ICU tests data, the pore water pressure-strain curves and effective confining pressure-strain curves of the ZT and STW loess were plotted, as shown in Figures 6 and 7.

**Figure 6.** Pore water pressure-strain curves of the ZT and STW loess. (**a**) ZT. (**b**) STW.

**Figure 7.** Effective confining pressure-strain curves of the ZT and STW loess. (**a**) ZT. (**b**) STW.

According to Figure 6, the pore water pressure curves of the ZT and STW loess exhibited almost the same change trend and were mainly divided into two stages: (1) Sharp rise stage: at the beginning of loading (the axial strain ε*<sup>a</sup>* was less than 3%), the pore water pressure rose rapidly, close to the confining pressure level. (2) Stable stage: after reaching the curve inflection point, the increment of the pore pressure was very small with the increase in the strain, and it tended to be gradually stable. The inflection point of the curve can be considered as the critical point of the pore water pressure change and is closely associated with the structural yield stress being reached [33]. Under the

undrained condition, the pore water pressure rose rapidly. Accompanied by the particle rearrangement, the spaced pore structure collapses, structural failure of the loess structure begins, and the effective stress gradually transfers to the pore water pressure [16]. When the collapse of the structure tended to end, the clay aggregates slipped and squeezed into the pore throat, which seriously hindered the discharge of the pore water. As a result, the pore water pressure remains at a relatively high level. The variation trend is accurately captured through the microstructure morphology analysis.

It can be observed from Figure 7 that under the different confining pressures, the effective confining pressure of the ZT and STW loess decreased sharply when the axial strain was less than 3%, and it tended to be gradually stable with the increase in the axial strain. The effective confining pressure of the ZT loess in the steady state increased with the increase in the initial confining pressure. Under the confining pressure of 250 kPa, the pore pressure of the STW loess increased more; hence, the effective confining pressure was the lowest in the steady state. Under the other confining pressure conditions, the effective confining pressure of the STW loess in the steady state also increased with the increase in the initial confining pressure. Combined with Figures 6 and 7, it is not difficult to determine that when the axial strain was less than 3%, the effective confining pressure decreased with the increase in the pore water pressure, and the pore water pressure and the effective confining pressure tend to be stable with the increase in the axial strain. Although the effective confining pressure at the end of the shear was a lower value (20–100 kPa), it did not reach zero. This means that the loess had not reached complete static liquefaction; however, it had decreased significantly in comparison to the initial confining pressure. At this point, the loess was in a very unstable state and had a strong flow plastic characteristic [14].

#### *4.2. Analysis of the Microstructure Characteristics before and after the ICU Tests*

Based on the aforementioned analysis of the ICU test results for the loess under the different confining pressures, the ZT and STW loess samples show prominent strain softening characteristics. After reaching the peak value, the partial stress decreased significantly; the pore water pressure increased sharply and remained high; and the effective confining pressure decreased rapidly. As a result, the ZT and STW loess have the possibility of static liquefaction. In order to find the relationship between the static liquefaction and the microstructure characteristics and to understand the internal mechanism of the static liquefaction in this area more clearly, the scanning electron microscope (SEM) images captured before and after the ICU tests were analyzed in this section.

In Figures 8 and 9(a1,a2) are SEM images of loess after the back pressure saturation treatment on the triaxial apparatus. Therefore, according to the experimental process, the SEM images of (a1) and (a2) reflect the micromorphology of loess before consolidated undrained triaxial tests. Moreover, (b)–(e) in Figures 8 and 9 are the SEM images of loess after consolidated undrained triaxial tests under different confining pressures. Through a comparative analysis of the microstructure changes before and after the ICU tests under the different confining pressures (Figures 8 and 9), the following results were determined. (1) Before the ICU tests, there were a large number of scattered spaced pores between the loess particles, which was confirmed to the microstructure characteristics of the loess [34]. After the ICU tests, the original pore throat channel of the spaced pores between the loess coarse particles was filled by fine particles, which can be also observed in the previous study [9]. Therefore, the pore throat channel was blocked, the boundary shape of the pore throat channel tended to be more complicated, and the intergranular pores were gradually dominated by mosaic pores. (2) Before the ICU tests, the contact between the coarse particles was mainly point-point and point-face. After the ICU tests, the appearance of fine particles, such as clay aggregates, in the contact position of the particles destroyed the contact state between the coarse particles of the loess. This resulted in the contact state between the coarse particles gradually transforming to face-face. The face-face contact relation mentioned in here includes direct face-face contact and indirect face-face contact [35,36]. (3) Under the confining pressures of 350 kPa and 450 kPa, the edges of the coarse particles peeled off from the fine particles. The same phenomenon was found in another study [37]. This indicates that the

spaced framework formed by overlapping of the coarse particles becomes fully contacted under the compression deformation; however, the framework structure cannot bear a greater load. At the contact point of the coarse particles, the cracks would appear due to the stress concentration and they develop towards the inside of the particles, which gradually produces fragmentation.

(**a1**) (**a2**)

(**b1**) (**b2**)

(**c1**) (**c2**)

**Figure 8.** *Cont*.

(**e1**) (**e2**)

**Figure 8.** SEM images (2000×) of the ZT and STW loess before and after the isotropic consolidated undrained (ICU) tests. (**a1**) Saturated (ZT). (**a2**) Saturated (STW). (**b1**) 150 kPa (ZT). (**b2**) 150 kPa (STW). (**c1**) 250 kPa (ZT). (**c2**) 250 kPa (STW). (**d1**) 350 kPa (ZT). (**d2**) 350 kPa (STW). (**e1**) 450 kPa (ZT). (**e2**) 450 kPa (STW).

**Figure 9.** SEM images (500×) of the ZT and STW loess before and after the ICU tests. (**a1**) Saturated (ZT). (**a2**) Saturated (STW). (**b1**) 150 kPa (ZT). (**b2**) 150 kPa (STW). (**c1**) 250 kPa (ZT). (**c2**) 250 kPa (STW). (**d1**) 350 kPa (ZT). (**d2**) 350 kPa (STW). (**e1**) 450 kPa (ZT). (**e2**) 450 kPa (STW).

#### **5. Discussion**

#### *5.1. Analysis of Capacity of Static Liquefaction*

The stress-strain curves, pore water pressure-strain curves, and the effective confining pressure-strain curves of the ZT and STW loess under the different confining pressures were analyzed. The loess under the different confining pressures shows prominent strain softening characteristics, but the deviator stress at the end of the shear was not zero. In addition, the pore water pressure rose rapidly under a small deformation (ε*<sup>a</sup>* < 3%), and the effective confining pressure decreased with the increase in the pore water pressure and it reached a lower value (20–100 kPa) at the end of the shear. Therefore, the saturated undisturbed loess may undergo static liquefaction under the undrained shear action; however, the possibility of static liquefaction still needs further discussion.

#### 5.1.1. Analysis Change Characteristics of the LPI

To further understand the effect of the confining pressure on the liquefaction resistance of the saturated undisturbed loess, the previous concept of the liquefaction potential index (LPI) was used for the analysis. Ng. C [38] defined the LPI as follows: *ψ ψ*

$$LPI = \frac{q\_{\text{max}} - q\_{\text{min}}}{q\_{\text{max}}} \tag{1}$$

where *q*max is the peak shear stress, and *q*min is the quasi-steady shear stress. This parameter is similar to the brittleness index (*IB*) proposed by Bishop (1967) [39] for the strain softening materials. It can be used to characterize the degree of reduction in the deviatoric stress of the saturated undisturbed loess under the different confining pressures; thus, it reflects more accurately the ability of the loess to resist the static liquefaction. When the *LPI* = 1, this indicates that there is complete static liquefaction of the soil. When the *LPI* = 0, this indicates that there is no static liquefaction of the soil. When 0 < *LPI* < 1, this indicates that the soil has a certain static liquefaction potential and the ability to resist static liquefaction weakens with the increase in the *LPI*. According to Formula (1), the *LPI* of the ZT and STW loess under the confining pressures of 150 kPa, 250 kPa, 350 kPa, and 450 kPa is calculated, and the results are shown in Figure 10a.

*ψ* **Figure 10.** Measurement index of the liquefaction capacity. (**a**) liquefaction potential index (LPI). (**b**) ψ*.* Been [40] defined the state parameter (ψ) as follows:

$$
\psi = e\_0 - e\_{\rm ss} \tag{2}
$$

′ *σ* ′ *σ* ′

′ *σ* ′-*σ* ′

where *e*<sup>0</sup> and *ess* are the void ratios in the initial and steady states, respectively, for the same effective mean normal stress. When ψ > 0, the soil is in a "loose" state, and the undrained shear can cause shear shrinkage. The larger the state parameter ψ is, the more obvious the shear shrinkage behavior of soil is and the stronger the liquefaction capacity is. According to Formula (2), the ψ of the ZT and STW loess for the confining pressures of 150 kPa, 250 kPa, 350 kPa, and 450 kPa were calculated, and the results are shown in Figure 10b.

It can be observed from Figure 10a that the *LPI*s of the ZT and STW loess were almost equal under the different confining pressures and range between 0.4 and 0.6. In addition, with the increase in the confining pressure, no clear monotonic change rule; this is similar to the results of previous ICU tests on saturated undisturbed loess in the South Jingyang platform [9,12,30]. At first, with the increase in the confining pressure from 150 to 250 kPa, the *LPI* shows an upward trend. At this point, the *LPI*s of the ZT and STW loess were 0.61 and 0.59, respectively. Then, with the increase in the confining pressure from 250 to 350 kPa, the *LPI* decreased to some extent. Finally, with the increase in the confining pressure from 350 to 450 kPa, the *LPI* remained unchanged. It can be determined from Figure 10b that ψ of the ZT and STW loess increased with the increase in the confining pressure from 250 to 450 kPa; however, it remained unchanged with the increase in the confining pressure from 150 to 250 kPa. Moreover, for a confining pressure of 150 kPa and 250 kPa, the LPI was larger; yet ψ was smaller. When the confining pressure was 350 kPa and 450 kPa, the LPI was smaller, but ψ was larger.

Furthermore, the possibility of static liquefaction existed in different layers and the buried depth of the undisturbed loess. In addition, there was no obvious relationship between static liquefaction and confining pressure. In fact, according to the research on the buried depth of the sliding surface of the loess landslides in the South Jingyang platform, it is determined that there is no clear loess stratigraphic unit for the sliding of the high-speed and long-distance loess landslide; however, there is a prominent sliding face for the conventional sand landslides [41,42]. Although no prominent sliding face was observed in the loess landslides of the South Jingyang platform, the loess-paleosol of the South Jingyang platform presents typical interbedded sequence characteristics. Therefore, an evident water retention effect was easily produced at the paleosol interface. This is because the paleosol layer was equivalent to a relatively dense aquiclude. Some landslides in this area do have the possibility of sliding along the paleosol layer [12,24], but the specific formation was not regular.

#### 5.1.2. Analysis of the Stress Path Evolution

The stress-strain curves and the pore pressure-strain curves of the ZT and STW loess indicate that the loess exhibited strong strain softening, and the pore water pressure remained relatively high. For the post peak state description and the steady-state evaluation of this kind of strain softening soil, the potential liquefaction state can be well recognized. Based on previous studies [43–45], the stability of the ZT and STW loess was evaluated by establishing a steady state line and an instability line. The stress paths of the ZT and STW loess were plotted using the experimental data from the ICU tests; the abscissa is the effective average stress *p* ′ = 1/3(σ<sup>1</sup> ′ + 2σ<sup>3</sup> ′ ) (MPa), and the ordinate is the deviator stress *q* ′ = σ<sup>1</sup> ′ − σ<sup>3</sup> ′ (MPa), as shown in Figure 11.

Figure 11 shows that the ZT and STW loess under the different confining pressures all have shear shrinkage phenomenon. With the increase in the confining pressure, the peak point of the effective stress path gradually increased. In addition, after the effective stress path of the ZT and STW loess reached the peak point, they all reached the steady state with an increase in the strain. The stress path was divided into the stable region and instability regions by the instability and steady-state lines. If the soil stress path entered the instability region between the instability line and the steady-state line, the soil was very likely to undergo static liquefaction. If the soil stress path is in the stable region below the instability line, static liquefaction of the soil is unlikely to occur [43,44,46]. Accordingly, it can be concluded that static liquefaction of the ZT and STW loess is more likely to occur under different confining pressures.

**Figure 11.** The stress paths of the ZT and STW loess. (**a**) ZT. (**b**) STW.

For the purpose of making the steady-state line pass through the steady-state points of the loess as much as possible, after trying to draw the steady-state lines, it was determined that this effect can be achieved when the steady-state line is a straight line without passing through the origin. Based on the stress path that from the ICU tests of the loess, some scholars [17,47] observed that the steady-state line is a straight line that does not pass through the origin. Although the uniqueness of the steady-state line of sand is controversial [48,49], most of it is a straight line that passes through the origin [43–45]. Loess still has cohesion in the steady-state stage [50], while the cohesion of sand is usually negligible. In addition, the macroscopic mechanical behavior of loess is the result of the coupling effect of the particle cementation strength and the friction strength. After overcoming the cementation strength between the loess particles, the greater the friction strength between the particles, the weaker the strain softening phenomenon [50]. The friction strength between the particles depends on the shape, size, and the relative position of the particles after structural reorganization [51]. The difference in the friction strength between the particles results in a different static liquefaction phenomenon after the peak value for the loess and sand; however, both display shear shrinkage characteristics.

#### *5.2. Correlation Analysis of the Macro and Micro Characteristics of Static Liquefaction*

Through the analysis of the mechanical behavior of ICU tests and the structural evolution, the following conclusions can be obtained. (1) The most obvious feature of the static liquefaction is that the pore water pressure rises sharply and it remains high. (2) The stable-state points fall into the instability region; however, the steady-state strength is not zero, which reflects that there is still a certain residual strength, and it should belong to incomplete static liquefaction. Furthermore, the typical characteristics of static liquefaction are all derived from the change in the loess internal structure [52,53]. As demonstrated from the comparison of the SEM images of the loess before and after the ICU tests (Figure 8), significant changes had taken place in the internal morphology; however, a quantitative analysis of the changes in the pore structure is still needed to accurately understand the relationship between the changes in the macroscopic mechanical behavior and the changes in the microstructure.

To further analyze the changes in the pore structure of the loess before and after the ICU tests, image analysis software (Image-pro Plus (IPP, Media Cybernetics Inc. Rockville, MD, USA)) was used to process the SEM images at 500 times magnification (Figure 9). This software can effectively distinguish the image from the background by threshold segmentation of the original SEM image, and then it obtains the determined value by performing a statistical analysis of the corresponding quantitative parameters. The specific operation flow is as follows. (1) Each of the original SEM images is uniformly corrected (including improving the image brightness and contrast, using the median filter to denoise the image to accurately define the pore boundary, using image binary processing technology to accurately identify the pore structure, etc.). Among them, the median filter is a noise reduction method that the middle value of the gray of the neighbor pixels replaces the gray value of the central

pixel, and is mainly used to eliminate mutational noise points and improve image quality. The above uniform corrections of SEM images are shown in Figure 12. (2) The quantitative parameters were selected for the statistical analysis, among which, the quantitative parameters involved in this study include the average pore diameter, shape ratio, roundness, and the pore contour fractal dimension. The concepts related to the quantitative parameters are described as follows, and the conceptual graph (Figure 12e) is shown for understanding adequately. (3) The quantitative characterization of the pore structure is achieved through the above parameters, and the detailed results are illustrated in Figure 13.


**Figure 12.** SEM images correction process. (**a**) Original SEM images. (**b**) Improving brightness and contrast. (**c**) Median filter. (**d**) Binary processing. (**e**) Quantitative parameters of different types pore structure.

μ

*π*

−

(**c**) (**d**) **Figure 13.** The quantitative parameters of the ZT and STW loess before and after the ICU tests. (**a**) Average pore diameter. (**b**) Shape ratio. (**c**) Roundness. (**d**) Pore contour fractal dimension.

It can be observed from Figure 13 that after the ICU tests, the average pore diameter decreased and the shape ratio was basically unchanged; however, the roundness and pore contour fractal dimension shows two different trends: the roundness and the pore contour fractal dimension of the ZT loess increased after the ICU tests. After the ICU tests, the roundness of the STW loess decreased, and the pore contour fractal dimension increased or decreased. It can be inferred that a large deformation (e.g., crushing, rotating, sliding, etc.) occurred in the aggregate spaced pore structure formed by overlapping of the coarse particles inside the loess; this resulted in the sudden and interlocking initial structural damage. The deformation of the pore structure essentially involves particle arrangement and reorganization [46,58]. In addition, the differences in the shape ratio, roundness, and pore contour fractal dimension before and after the ICU tests reflect the uncertainty rearrangement and reorganization of the particles. However, the difference in the average pore diameter before and after the ICU tests shows that the fine particle aggregates can enter the original pore throat in the process of rearrangement and reorganization of the particles. This can be observed intuitively in Figure 8b–e and is also verified by other studies [9,17,57]. At the beginning of the particles rearrangement and reorganization, the water transport channel formed by spaced pores is blocked to a certain extent. Hence, the macroscopic mechanical behavior is characterized by a rapid decrease in the deviatoric stress after the peak stress

*ε*

and a slow increase in the pore water pressure (rapid strain softening stage). Subsequently, the process of particle rearrangement and reorganization tends to be stable, and the migration and agglomeration of fine particles under the action of pore water pressure continue [26]. The synergistic effect of particle rearrangement and reorganization and the pore water pressure causes the agglomeration blocking of the reconstructed pore throat. Therefore, the macroscopic mechanical behavior shows that the deviatoric stress gradually tended to shift to the steady-state stage, and the pore water pressure was maintained at a high level (stable strain softening stage). In addition, the macroscopic mechanical behavior of the loess before reaching the peak strength displayed a rapid increase in the deviatoric stress and the pore water pressure (rapid strain hardening stage). It is reasonable to speculate that the widespread spaced pore structure in loess often requires a certain compression deformation to ensure full contact and the formation of a force chain [59]. This deformation process can also compress the pore water transport channel.

Based on the different pore diameter ranges of the loess (macropore > 32 µm, mesopore 8–32 µm, small pore 2–8 µm, and micropore < 2 µm) proposed by Lei [60], and the IPP was used to perform a statistical analysis on the different pore proportion, the detailed results are presented in Figure 14. As illustrated in Figure 14, after the ICU tests, the different pore proportion change characteristics of the ZT and STW loess were as follows: the macropore and mesopore decreased, the small pore increased slightly, and the micropore increased significantly. The micropores are mainly cement pores; the small pores are mainly mosaic pores and a small amount of cement pores; the mesopores are mainly spaced pores and some mosaic pores; and the macropores mainly includes root holes, wormholes, and fractures [60]. This shows that during the particles reorganization, the intergranular pores are gradually filled by fine particles such as clay aggregates; thus, resulting in an increase in the proportion of the micropore and small pore, which is also illustrated by the decrease in the average pore diameter and Figure 8b–e after the ICU tests. After the ICU tests, most of the above pore structure quantitative parameters (e.g., average pore diameter, shape ratio, roundness, pore contour fractal dimension, and the proportion of the different pores) did not show an obvious monotonic change law with the increase in the confining pressure. The static liquefaction characteristics of the loess come from the change in the internal structure of the loess, which again reflects that there was no clear relationship between the confining pressure and the static liquefaction of the loess in the South Jingyang platform.

**Figure 14.** Different pore proportion of the ZT and STW loess before and after the ICU tests. (**a**) ZT. (**b**) STW.

The characteristics of the two grain phases (the gain composition was mainly silt and clay) of the ZT and STW loess shown in Table 2 and the pore proportions shown in Figure 14 indicate that the spaced structure formed by the accumulation and overlapping of silt constituted the entire structural framework of the loess. As the particle size of clay was sufficiently small compared with that of silt, the van der Waals force and repulsive force between the clay counteracted the filling effect due to gravity in the relatively loose silt skeleton structure. Therefore, the clay forms relatively closed agglomerates through agglomeration [46]; however, the high content of clay rendered the clay aggregates inlaid or resulted in the filling of the pores or contact areas in the spaced structure (Figures 8 and 9). Therefore, the structure of the loess, which is composed of silt and is filled with clay aggregates, was similar to that of sand, which is composed of coarse grains and is filled with fine grains such as silt [44]. According to a hypothetical model of the interaction between fine grains and coarse grains proposed by Lade and Yamamuro [44], when the fine grains occupy the pores between the coarse grains, they only increase the material density and have negligible influence on the behavior of the soil. When the fine grains occupy the location near the contact of the coarse particles, under the isotropic shear action, the fine grains tend to slide into the pores between the coarse grains. This promotes the rearrangement of the coarse particle structure, it increases the volume shrinkage of the soil, and it produces a greater liquefaction potential under the undrained condition. Therefore, the structure characteristics of the loess provides better structural conditions for the occurrence of static liquefaction phenomenon, which makes it easy to have a structural reorganization under the action of an external load. Clay aggregates would enter the spaced structure that is composed of overlapping silt and hinders the dissipation of pore water [9,17,57].

#### **6. Conclusions**

According to the stress-strain curves, pore water pressure-strain curves, and the regions of the stress path, the static liquefaction of the loess in the South Jingyang platform is incomplete. Through a comparative analysis of the microstructure characteristics before and after the ICU tests, the evolution process and internal mechanism of the static liquefaction were revealed. The main conclusions are as follows:


pore to be blocked to a certain extent. The macroscopic mechanical behavior shows rapid strain softening. When ε*<sup>a</sup>* = 10–15%, the gradually stable rearrangement and reorganization of the particles and the migration and agglomeration of finer grains under the effect of the pore water pressure caused agglomeration blocking of the reconstructed pore throat. The macroscopic mechanical behavior indicates stable strain softening.

**Author Contributions:** Writing—Original Draft Preparation, R.-X.Y.; Funding Acquisition, J.-B.P.; Writing—Review and Editing, R.-X.Y. and J.-Y.Z.; Data Curation, S.-k.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the Major Program of National Natural Science Foundation of China (Grant No. 41790441), the National Natural Science Foundation of Shaanxi Province, China (Grant No. 2018JQ5124), and the Foundation of Key Laboratory of Western Mineral Resources and Geological Engineering of Ministry of Education, Chang' an University (Grant No. 300102268503).

**Acknowledgments:** We thank the entire team for their efforts to improve the quality of the article. At the same time, we would like to thank editor for his timely handling of the manuscripts.

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

#### **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/).
