# Effect of Vegetation Roots on the Threshold of Slope Instability Induced by Rainfall and Runoff

Authors: Gang Huang; Mingxin Zheng; Jing Peng
Journal: Geofluids
 (2021)

---

## Abstract

Although vegetation is increasingly used to mitigate landslide risks, how vegetation roots affect the landslide threshold of slope has rarely been explored, particularly in the case of lateral runoff. In this study, we established a two-dimensional saturated-unsaturated infiltration equation considering the hydraulic effects of vegetation roots. The analytical solution for the shallow unsaturated two-dimensional coupled infiltration of vegetated slope (VS) was obtained by a Fourier transform technique. The numerical method was used to evaluate the stability of VS caused by four root architectures, the rainfall amount, and the rainfall duration. Subsequently, the transformation law in runoff, vegetation evaporation, and landslide threshold was analyzed. The results indicate that the factor of safety (FOS) increases with increasing drying time and decreases with increasing depth; the minimum FOS is at the junction of the root-rootless zone. Runoff and vegetation evaporation are favorable for the shallow stability of VS. The time of the safe area is 35 h for rainfall amount 500 m in the uniformly root clay slope. Moreover, four landslide threshold curves that reflected the root architecture, rainfall amount, and rainfall duration are developed, which are more realistic than those created using one-dimensional instability modeling.


---

## Body

## 1. Introduction

Rainfall-induced landslides have occurred frequently in mountainous environments worldwide [1]. Among the types of rainfall-induced landslides, shallow landslides of the vegetated slope (VS) have resulted in damage to infrastructure and ecological environment, economic losses, and harm to human lives. They are usually due to the infiltration of rainfall, which reduces the shear strength of the slope soil and eventually causes the slope to fail [2–4]. Heavy rainfall for a short time or moderate rainfall for a long time can lead to a large number of landslides, most of which occur in shallow locations above ground level [5]. According to statistics, VS accounted for 38% of the 133 shallow landslides in northern Switzerland since 2015 [6]. Therefore, it is necessary to carry out the study of rainfall-induced landslides of VS.In the current study, the prediction of rainfall-induced landslides has attracted extensive attention [7, 8]. Landslide threshold is a common method for predicting and early warning of landslides, which can effectively predict the catastrophic consequences of the slope in heavy rainfall. The main influence of the landslide threshold was normalized to the intensity and duration of rainfall [9]. A wide range of landslide thresholds have been developed [10, 11]. Recent studies had shown that vegetation could be a considerable spatial predictor of landslide occurrence [12]. The presence of vegetation could affect the occurrence of shallow instability [13, 14]. Moreover, vegetation has been proven to be available for shallow slope stability, which is due to a series of hydrodynamic effects of roots and leaves [15]. A large number of studies have shown that the main effects of vegetation on slope stability are mechanical reinforcement and hydromechanical of roots [16]. Mechanical reinforcement refers to increasing the shear strength of the slope through the mechanical properties of roots, which improve the stability of slope soil and reduce the occurrence of the shallow landslide [17, 18]. The hydrological effect refers to the absorption and transpiration of roots to increase the substrate suction of soil, thereby improving the stability of shallow slopes [19]. At present, research on the effects of hydrological on the stability of shallow slopes has become a hot topic [20, 21]. Relevant scholars in Hong Kong have conducted experiments to study the effect of roots on pore water pressure and the factor of safety (FOS) of a shallow slope in the atmospheric environment, moving towards the fundamental understanding of atmosphere-plant-soil interaction [22–24]. A mathematical model for the rate of root water uptake has been developed, including the root growth rate considering topographic factors, the type of vegetation, and climatic parameters [25]. Moreover, a one-dimensional hydromechanical model for the effect of root architecture on the mechanical and hydrological behaviour of soil has been developed [26]. The one-dimensional model can consider the water absorption of vegetation [27], but there is no way to consider the effect of slope scale, especially the rainfall-induced lateral runoff.Meanwhile, increasing studies have focused on the impact of runoff on rainfall thresholds [28, 29]. The influence of slope runoff at the fluid-solid interface should be considered [30]. The lateral runoff affects the cumulative seepage flow and head pressure of the slope [31]. A landslide analysis procedure combining hydrological model and landslide susceptibility model is established [32]. The influence of lateral runoff is demonstrated by an in-house physical-process-based shallow landslide model [33]. A GIS-based approach for analyzing and assessing the simplified landslide susceptibility in the Entella River catchment is presented [34]. Lateral runoff is the susceptibility factor that triggers the shallow landslide [35]. Understanding the rainfall threshold under lateral runoff and hydrological behaviour could be especially useful, but these links are rarely explored due to the difficulty of modeling the two-dimensional coupled infiltration [36].In the light of the above discussion, this study is aimed at evaluating the shallow stability of VS under rainfall and lateral runoff. In this paper, the Richards equation and limit equilibrium method are used to establish a two-dimensional hydromechanical model and simulate the shallow instability of VS. Analytical solutions for the shallow unsaturated two-dimensional coupled infiltration of finite extent were obtained by a Fourier transform technique. The effects of the root architectures, rainfall amount, and rainfall duration in shallow landslides were studied, and the landslide thresholds of different root architectures were discussed.


## 2. Two-Dimensional Coupled Hydromechanical Model

### 2.1. Infiltration Equation

Rainfall infiltration is the major inducement of shallow instability of VS, which is a typical unsaturated fluid-structure coupling phenomenon. As shown in Figure1, the transpiration of VS can reduce the shallow soil water content, increase matrix suction, and improve slope stability; rainfall first infiltrates the shallow root soil and then enters the unsaturated zone under the action of water pressure and unsaturated infiltration. Changes in the seepage field of the soil cause the pore water pressure in the unsaturated zone to increase, forming a temporarily saturated zone and resulting in the instability of VS.Figure 1
Schematic diagram of the hydraulic action of VS.The water movement in the unsaturated zone of the VS can be expressed by Richards’ equation, which satisfies the conservation of mass and Darcy’s law. Considering the suction effect of the vegetation roots, the sink termSRx,z,t is incorporated into the two-dimensional seepage Equation (26):
(1)n∂Sr∂ψ+ηSrγw1+2νxEx+1Ez∂ψ∂t=∂∂xkx∂ψ∂x+∂∂zkz∂ψ∂z+cosα∂kzψ∂z−SRx,zΓx,z−h2,where n is the porosity of the soil; Sr is the saturation of the soil; η is the void factor, η=1−n2/1−no, no is the initial porosity; γw is the water severity; Ex and Ez are the elastic moduli of the soil in the x and z directions, respectively; νx is Poisson’s ratio of the soil; ψ is the pore pressure head; kx and kz are the permeability coefficients of the soil in the x and z directions, respectively, assuming that the soil is isotropic, kx=kz; SRx,z is the water absorption rate of root. Γx,z−h2 is defined as
(2)Γx,z−h2=00≤x≤L,0≤z<h2,10≤x≤L,h2≤x<h1+h2,Rootlesszone,Rootzone,where h1 is the vertical distance between the lower edge of the root zone and the slope surface, h2 is the vertical distance between the lower edge of the root zone and the groundwater level, and L is length of the zone.According to Fatahi et al.’s [25] suggestion, the water absorption of the roots should consider the three aspects of soil water absorption, root distribution, and vegetation transpiration. The water absorption function SRx,z of the roots can be expressed as
(3)SRx,z=α⋅PT⋅Gz,where α is the function of soil transpiration reduction. According to Feddes [37] and Chiu et al. [32], the value of α is 1, Gz is the density function of the vegetation roots, and PT is the transpiration of vegetation.The vegetated transpirationPT is related to the leaf index, and the expression can be as follows:
(4)PT=PET⋅C,where PET is the potential evaporation of the vegetation; C is the canopy coverage factor, C=1−exp−b⋅LAI; b is the radiation parameter of vegetation leaves; and LAI represents the leaf index of vegetation.The roots of vegetation are mainly composed of the primary root and secondary root. The primary root plays the role of skeleton and the fibrous root absorbs water. According to Nyambayo and Potts [23] and Figure 1, the expression of four different root architectures is as follows:
(5)Gz=1h1,Uniformly,expz−H−h1−1exph1−h1−1,Exponential,2h1z−Hh1+1,Triangular,6h13z2+2H−h1z+h1H−H2,Oval.For the solution of the vegetation seepage equation, it is necessary to obtain the water absorption characteristics of the vegetation roots and the hydraulic characteristics of the root soil. Void ratio caused by vegetation roots is proposed by Ni et al. [24], as follows:
(6)e=e0−Rv1+e01−Rv1+e0,where e0 is the porosity of the rootless zone and Rv is the volume ratio of the roots.The relationship between the saturated permeability coefficient and porosity of the VS is proposed by Yin [38], as follows:
(7)ks=λdecexpρ⋅eexpρ⋅e0ks0,where λdec is the increasing porosity ratio due to the increase of the roots, ranging from 1.3 to 6.5; ρ is the fitting parameter; and ks0 is the saturation permeability coefficient.Gardner’s model [39] is adopted for the permeability coefficient and soil-water characteristic curve of unsaturated VSs, as follows:
(8)kx=kz=ksexpβψ,(9)Se=θ−θrθs−θr=S−SrSs−Sr=expβψ,where β is the parameter related to the intake pressure.Equations (7), (8), and (9) are incorporated into Equation (1) to obtain two-dimensional partial differential equations as follows:
(10)Pexpβψ∂ψ∂t=∂∂xKexpβψ∂ψ∂x+∂∂zKexpβψ∂ψ∂z+βKexpβψ∂ψ∂z−SRx,zΓx,z−h2,where P=βθsn+ηγw1/Hs+1+2vx/Hx.Setting Kirchhoff variables [40]
(11)Φ=expβψ,Equation (10) can be rearranged as
(12)P∂Φβ∂t=kβ∂2Φ∂x2+kβ∂2Φ∂z2+k∂ψ∂z−SRx,zΓx,z−h2.To solve Equation (12), the dimensionless variables Z=βz, X=βx, and T=βkt/P and dependent variables as VX,Z,T are obtained as
(13)V=ψexpZ2+T4−βTPZ22h1k,Uniformly,ψexpZ2+T4−βTPexpZ−H−h1−0.5βTPZ2kexph1−h1−1,Exponential,ψexpZ2+T4−βTPkh12Z3−6HZ26h1+Z2,Triangular,ψexpZ2+T4−TPβkh132Z4+2H−h1Z3+3h1H−H2Z2,Oval.From [12, 13], Equation (12) can be rearranged as
(14)∂V∂T=∂2V∂X2+∂2V∂Z2.


### 2.2. Initial Condition and Boundary Conditions

For the solution of Equation (14), the solution and boundary conditions must be available. If the shape of the solution is too complicated, the two-dimensional seepage analysis cannot obtain an analytical solution. For the two-dimensional seepage analysis area, the rectangular area is used to solve many practical engineering problems. Therefore, the solution area in this paper is a rectangle with a length of L′=βL and a height of H′=βH. The initial condition and four boundary conditions applicable to two-dimensional transient seepage can be expressed as
(15)FX,Z=VX,Z,TT=0=ψiexpZ2,0≤Z<h2,ψiexpZ2−βTPZ22h1k,Uniformly,ψiexpZ2−βTPexpZ−H−h1−0.5βTPZ2kexph1−h1−1,Exponential,ψiexpZ2−βTPkh12Z3−6HZ26h1+Z2,Triangular,ψiexpZ2−TPβkh132Z4+2H−h1Z3+3h1H−H2Z2,Oval,h2≤Z≤h1+h2,∂V∂XX=0=0,∂V∂XX=L=0,VZ=0=expT4,∂V∂Z+0.5VZ=H=qgkexpH2+T4,where qg is the rainfall intensity.


### 2.3. Analytical Solution of Root Zone

#### 2.3.1. Uniformly Roots

Using the Fourier integral transformation method [41], the general solution expressions of the pressure head ψX,Z,T in the uniformly rooted region and the rootless region can be obtained, respectively.
(16)ψX,Z,T=1βlnexp−Z2×βTPZ22h1k+∑m=0∞∑n=0∞K1χm,XK2λn,Z+exp−χm2+λn2+0.25T×F¯χm,λn+∫t=0Texpχm2+λn2tAχm,λn,tdt,where χm=mπ/L′,m=1,2,3,⋯.The kernel function can be expressed as(17)K1χm,X=2L′cosχmX,m=1,2,3,⋯,(18)K1λn,Z=2λn2+0.25H′λn2+0.25+0.51/2sinλnZ,n=1,2,3,⋯,(19)Aχm,λn,t=dK2λn,ZdZZ=0∫x′=0L′K1χm,xexpT4dx′+K2λn,Z0.5Z=H′∫x′=0L′K1χm,x′i−ickexpB2+T4+βTPZ22h1kdx′,(20)F¯χm,λn=∫x′=0L′∫z′=0H′K1χm,x′K2λm,zψiexpz′2dx′dz′,where χm=mπ/L′,m=1,2,3,⋯ and λm=nπ/H′,n=1,2,3,⋯.Substituting formulas (17) and (18) into formulas (19) and (20), it can be changed to
(21)Aχm,λn,t=2λn2+0.25L′H′λn2+0.25+0.5L′1/2expT/4sinnχmL′χmL′×λn+qgkexpH′2+βTPZ22h1ksinλnH′,F¯χm,λn=2λn2+0.25H′λn2+0.25+0.51/2expT/4sinχmL′χmL′×0.5−βsinλnH′expH′0.5−βλn2+0.5−β−λnexpH′0.5−βcosλnH′−1λn2+0.5−β2.Thus, the general solution expression of pore water pressure in the root zone can be obtained as(22)ψx,z,t=1βlnexpTPh1+h22kx−βx+βz2×∑m=0∞∑n=0∞C1cosχm,βxsinλn,βz×exp−χm2+λn2+0.25⋅βktP+C2cosχm,βxsinλn,βz,where
(23)C1=2λn2+0.25β2Lh1+h2λn2+0.25+0.5βL,C2=2λn2+0.25βh1+h2λn2+0.25+0.51/2×expT/4sinχmβLχmβL×0.5-βsinλnβh1+h2expβh1+h20.5-βλn2+0.5-β-λnexpβh1+h20.5-βcosλnβh1+h2-1λn2+0.5-β2+2λn2+0.25β2Lh1+h2λn2+0.25+0.5βL1/2×sinnχmβL2χm2+2λn2χmβL×exp14χm2+4λn2+2Tχm2+λn2-exp14χm2+4λn2×λn+qgkexpβh12+β2TPh1+h22kxsinλnβh1+h2.


#### 2.3.2. Rootless Zone

The pore pressure head of the rootless zone according to the above method is(24)ψx,z,t=1βlnexp−βx−βz2×∑m=0∞∑n=0∞C3cosχm,βxsinλn,βz×exp−χm2+λn2+0.25⋅βkt/P+C4cosχm,βxsinλn,βz,where(25)C3=2λn2+0.25β2h2Lλn2+0.25+0.5βL,C4=2λn2+0.25βh2λn2+0.25+0.51/2expT/4sinχmβLχmβL×0.5−βsinλnβh2expβh20.5−βλn2+0.5−β−λnexpβh20.5−βcosλnβh2−1λn2+0.5−β2+2λn2+0.25β2h2Lλn2+0.25+0.5βL1/2sinnχmβL2χm2+2λn2χmβLexp14χm2+4λn2+2Tχm2+λn2−exp14χm2+4λn2×λn+qgkexpβh22sinλnβh2.According to Robida’s law, the values ofC1,C2, C3,and C4 tend to zero when the values of m and n tend to infinity, and the water pressure tends to converge.


## 2.1. Infiltration Equation

Rainfall infiltration is the major inducement of shallow instability of VS, which is a typical unsaturated fluid-structure coupling phenomenon. As shown in Figure1, the transpiration of VS can reduce the shallow soil water content, increase matrix suction, and improve slope stability; rainfall first infiltrates the shallow root soil and then enters the unsaturated zone under the action of water pressure and unsaturated infiltration. Changes in the seepage field of the soil cause the pore water pressure in the unsaturated zone to increase, forming a temporarily saturated zone and resulting in the instability of VS.Figure 1
Schematic diagram of the hydraulic action of VS.The water movement in the unsaturated zone of the VS can be expressed by Richards’ equation, which satisfies the conservation of mass and Darcy’s law. Considering the suction effect of the vegetation roots, the sink termSRx,z,t is incorporated into the two-dimensional seepage Equation (26):
(1)n∂Sr∂ψ+ηSrγw1+2νxEx+1Ez∂ψ∂t=∂∂xkx∂ψ∂x+∂∂zkz∂ψ∂z+cosα∂kzψ∂z−SRx,zΓx,z−h2,where n is the porosity of the soil; Sr is the saturation of the soil; η is the void factor, η=1−n2/1−no, no is the initial porosity; γw is the water severity; Ex and Ez are the elastic moduli of the soil in the x and z directions, respectively; νx is Poisson’s ratio of the soil; ψ is the pore pressure head; kx and kz are the permeability coefficients of the soil in the x and z directions, respectively, assuming that the soil is isotropic, kx=kz; SRx,z is the water absorption rate of root. Γx,z−h2 is defined as
(2)Γx,z−h2=00≤x≤L,0≤z<h2,10≤x≤L,h2≤x<h1+h2,Rootlesszone,Rootzone,where h1 is the vertical distance between the lower edge of the root zone and the slope surface, h2 is the vertical distance between the lower edge of the root zone and the groundwater level, and L is length of the zone.According to Fatahi et al.’s [25] suggestion, the water absorption of the roots should consider the three aspects of soil water absorption, root distribution, and vegetation transpiration. The water absorption function SRx,z of the roots can be expressed as
(3)SRx,z=α⋅PT⋅Gz,where α is the function of soil transpiration reduction. According to Feddes [37] and Chiu et al. [32], the value of α is 1, Gz is the density function of the vegetation roots, and PT is the transpiration of vegetation.The vegetated transpirationPT is related to the leaf index, and the expression can be as follows:
(4)PT=PET⋅C,where PET is the potential evaporation of the vegetation; C is the canopy coverage factor, C=1−exp−b⋅LAI; b is the radiation parameter of vegetation leaves; and LAI represents the leaf index of vegetation.The roots of vegetation are mainly composed of the primary root and secondary root. The primary root plays the role of skeleton and the fibrous root absorbs water. According to Nyambayo and Potts [23] and Figure 1, the expression of four different root architectures is as follows:
(5)Gz=1h1,Uniformly,expz−H−h1−1exph1−h1−1,Exponential,2h1z−Hh1+1,Triangular,6h13z2+2H−h1z+h1H−H2,Oval.For the solution of the vegetation seepage equation, it is necessary to obtain the water absorption characteristics of the vegetation roots and the hydraulic characteristics of the root soil. Void ratio caused by vegetation roots is proposed by Ni et al. [24], as follows:
(6)e=e0−Rv1+e01−Rv1+e0,where e0 is the porosity of the rootless zone and Rv is the volume ratio of the roots.The relationship between the saturated permeability coefficient and porosity of the VS is proposed by Yin [38], as follows:
(7)ks=λdecexpρ⋅eexpρ⋅e0ks0,where λdec is the increasing porosity ratio due to the increase of the roots, ranging from 1.3 to 6.5; ρ is the fitting parameter; and ks0 is the saturation permeability coefficient.Gardner’s model [39] is adopted for the permeability coefficient and soil-water characteristic curve of unsaturated VSs, as follows:
(8)kx=kz=ksexpβψ,(9)Se=θ−θrθs−θr=S−SrSs−Sr=expβψ,where β is the parameter related to the intake pressure.Equations (7), (8), and (9) are incorporated into Equation (1) to obtain two-dimensional partial differential equations as follows:
(10)Pexpβψ∂ψ∂t=∂∂xKexpβψ∂ψ∂x+∂∂zKexpβψ∂ψ∂z+βKexpβψ∂ψ∂z−SRx,zΓx,z−h2,where P=βθsn+ηγw1/Hs+1+2vx/Hx.Setting Kirchhoff variables [40]
(11)Φ=expβψ,Equation (10) can be rearranged as
(12)P∂Φβ∂t=kβ∂2Φ∂x2+kβ∂2Φ∂z2+k∂ψ∂z−SRx,zΓx,z−h2.To solve Equation (12), the dimensionless variables Z=βz, X=βx, and T=βkt/P and dependent variables as VX,Z,T are obtained as
(13)V=ψexpZ2+T4−βTPZ22h1k,Uniformly,ψexpZ2+T4−βTPexpZ−H−h1−0.5βTPZ2kexph1−h1−1,Exponential,ψexpZ2+T4−βTPkh12Z3−6HZ26h1+Z2,Triangular,ψexpZ2+T4−TPβkh132Z4+2H−h1Z3+3h1H−H2Z2,Oval.From [12, 13], Equation (12) can be rearranged as
(14)∂V∂T=∂2V∂X2+∂2V∂Z2.


## 2.2. Initial Condition and Boundary Conditions

For the solution of Equation (14), the solution and boundary conditions must be available. If the shape of the solution is too complicated, the two-dimensional seepage analysis cannot obtain an analytical solution. For the two-dimensional seepage analysis area, the rectangular area is used to solve many practical engineering problems. Therefore, the solution area in this paper is a rectangle with a length of L′=βL and a height of H′=βH. The initial condition and four boundary conditions applicable to two-dimensional transient seepage can be expressed as
(15)FX,Z=VX,Z,TT=0=ψiexpZ2,0≤Z<h2,ψiexpZ2−βTPZ22h1k,Uniformly,ψiexpZ2−βTPexpZ−H−h1−0.5βTPZ2kexph1−h1−1,Exponential,ψiexpZ2−βTPkh12Z3−6HZ26h1+Z2,Triangular,ψiexpZ2−TPβkh132Z4+2H−h1Z3+3h1H−H2Z2,Oval,h2≤Z≤h1+h2,∂V∂XX=0=0,∂V∂XX=L=0,VZ=0=expT4,∂V∂Z+0.5VZ=H=qgkexpH2+T4,where qg is the rainfall intensity.


## 2.3. Analytical Solution of Root Zone

### 2.3.1. Uniformly Roots

Using the Fourier integral transformation method [41], the general solution expressions of the pressure head ψX,Z,T in the uniformly rooted region and the rootless region can be obtained, respectively.
(16)ψX,Z,T=1βlnexp−Z2×βTPZ22h1k+∑m=0∞∑n=0∞K1χm,XK2λn,Z+exp−χm2+λn2+0.25T×F¯χm,λn+∫t=0Texpχm2+λn2tAχm,λn,tdt,where χm=mπ/L′,m=1,2,3,⋯.The kernel function can be expressed as(17)K1χm,X=2L′cosχmX,m=1,2,3,⋯,(18)K1λn,Z=2λn2+0.25H′λn2+0.25+0.51/2sinλnZ,n=1,2,3,⋯,(19)Aχm,λn,t=dK2λn,ZdZZ=0∫x′=0L′K1χm,xexpT4dx′+K2λn,Z0.5Z=H′∫x′=0L′K1χm,x′i−ickexpB2+T4+βTPZ22h1kdx′,(20)F¯χm,λn=∫x′=0L′∫z′=0H′K1χm,x′K2λm,zψiexpz′2dx′dz′,where χm=mπ/L′,m=1,2,3,⋯ and λm=nπ/H′,n=1,2,3,⋯.Substituting formulas (17) and (18) into formulas (19) and (20), it can be changed to
(21)Aχm,λn,t=2λn2+0.25L′H′λn2+0.25+0.5L′1/2expT/4sinnχmL′χmL′×λn+qgkexpH′2+βTPZ22h1ksinλnH′,F¯χm,λn=2λn2+0.25H′λn2+0.25+0.51/2expT/4sinχmL′χmL′×0.5−βsinλnH′expH′0.5−βλn2+0.5−β−λnexpH′0.5−βcosλnH′−1λn2+0.5−β2.Thus, the general solution expression of pore water pressure in the root zone can be obtained as(22)ψx,z,t=1βlnexpTPh1+h22kx−βx+βz2×∑m=0∞∑n=0∞C1cosχm,βxsinλn,βz×exp−χm2+λn2+0.25⋅βktP+C2cosχm,βxsinλn,βz,where
(23)C1=2λn2+0.25β2Lh1+h2λn2+0.25+0.5βL,C2=2λn2+0.25βh1+h2λn2+0.25+0.51/2×expT/4sinχmβLχmβL×0.5-βsinλnβh1+h2expβh1+h20.5-βλn2+0.5-β-λnexpβh1+h20.5-βcosλnβh1+h2-1λn2+0.5-β2+2λn2+0.25β2Lh1+h2λn2+0.25+0.5βL1/2×sinnχmβL2χm2+2λn2χmβL×exp14χm2+4λn2+2Tχm2+λn2-exp14χm2+4λn2×λn+qgkexpβh12+β2TPh1+h22kxsinλnβh1+h2.


### 2.3.2. Rootless Zone

The pore pressure head of the rootless zone according to the above method is(24)ψx,z,t=1βlnexp−βx−βz2×∑m=0∞∑n=0∞C3cosχm,βxsinλn,βz×exp−χm2+λn2+0.25⋅βkt/P+C4cosχm,βxsinλn,βz,where(25)C3=2λn2+0.25β2h2Lλn2+0.25+0.5βL,C4=2λn2+0.25βh2λn2+0.25+0.51/2expT/4sinχmβLχmβL×0.5−βsinλnβh2expβh20.5−βλn2+0.5−β−λnexpβh20.5−βcosλnβh2−1λn2+0.5−β2+2λn2+0.25β2h2Lλn2+0.25+0.5βL1/2sinnχmβL2χm2+2λn2χmβLexp14χm2+4λn2+2Tχm2+λn2−exp14χm2+4λn2×λn+qgkexpβh22sinλnβh2.According to Robida’s law, the values ofC1,C2, C3,and C4 tend to zero when the values of m and n tend to infinity, and the water pressure tends to converge.


## 2.3.1. Uniformly Roots

Using the Fourier integral transformation method [41], the general solution expressions of the pressure head ψX,Z,T in the uniformly rooted region and the rootless region can be obtained, respectively.
(16)ψX,Z,T=1βlnexp−Z2×βTPZ22h1k+∑m=0∞∑n=0∞K1χm,XK2λn,Z+exp−χm2+λn2+0.25T×F¯χm,λn+∫t=0Texpχm2+λn2tAχm,λn,tdt,where χm=mπ/L′,m=1,2,3,⋯.The kernel function can be expressed as(17)K1χm,X=2L′cosχmX,m=1,2,3,⋯,(18)K1λn,Z=2λn2+0.25H′λn2+0.25+0.51/2sinλnZ,n=1,2,3,⋯,(19)Aχm,λn,t=dK2λn,ZdZZ=0∫x′=0L′K1χm,xexpT4dx′+K2λn,Z0.5Z=H′∫x′=0L′K1χm,x′i−ickexpB2+T4+βTPZ22h1kdx′,(20)F¯χm,λn=∫x′=0L′∫z′=0H′K1χm,x′K2λm,zψiexpz′2dx′dz′,where χm=mπ/L′,m=1,2,3,⋯ and λm=nπ/H′,n=1,2,3,⋯.Substituting formulas (17) and (18) into formulas (19) and (20), it can be changed to
(21)Aχm,λn,t=2λn2+0.25L′H′λn2+0.25+0.5L′1/2expT/4sinnχmL′χmL′×λn+qgkexpH′2+βTPZ22h1ksinλnH′,F¯χm,λn=2λn2+0.25H′λn2+0.25+0.51/2expT/4sinχmL′χmL′×0.5−βsinλnH′expH′0.5−βλn2+0.5−β−λnexpH′0.5−βcosλnH′−1λn2+0.5−β2.Thus, the general solution expression of pore water pressure in the root zone can be obtained as(22)ψx,z,t=1βlnexpTPh1+h22kx−βx+βz2×∑m=0∞∑n=0∞C1cosχm,βxsinλn,βz×exp−χm2+λn2+0.25⋅βktP+C2cosχm,βxsinλn,βz,where
(23)C1=2λn2+0.25β2Lh1+h2λn2+0.25+0.5βL,C2=2λn2+0.25βh1+h2λn2+0.25+0.51/2×expT/4sinχmβLχmβL×0.5-βsinλnβh1+h2expβh1+h20.5-βλn2+0.5-β-λnexpβh1+h20.5-βcosλnβh1+h2-1λn2+0.5-β2+2λn2+0.25β2Lh1+h2λn2+0.25+0.5βL1/2×sinnχmβL2χm2+2λn2χmβL×exp14χm2+4λn2+2Tχm2+λn2-exp14χm2+4λn2×λn+qgkexpβh12+β2TPh1+h22kxsinλnβh1+h2.


## 2.3.2. Rootless Zone

The pore pressure head of the rootless zone according to the above method is(24)ψx,z,t=1βlnexp−βx−βz2×∑m=0∞∑n=0∞C3cosχm,βxsinλn,βz×exp−χm2+λn2+0.25⋅βkt/P+C4cosχm,βxsinλn,βz,where(25)C3=2λn2+0.25β2h2Lλn2+0.25+0.5βL,C4=2λn2+0.25βh2λn2+0.25+0.51/2expT/4sinχmβLχmβL×0.5−βsinλnβh2expβh20.5−βλn2+0.5−β−λnexpβh20.5−βcosλnβh2−1λn2+0.5−β2+2λn2+0.25β2h2Lλn2+0.25+0.5βL1/2sinnχmβL2χm2+2λn2χmβLexp14χm2+4λn2+2Tχm2+λn2−exp14χm2+4λn2×λn+qgkexpβh22sinλnβh2.According to Robida’s law, the values ofC1,C2, C3,and C4 tend to zero when the values of m and n tend to infinity, and the water pressure tends to converge.


## 3. Shallow Slope Stability Analysis

Slope instability of unsaturated soil vegetation caused by rainfall mainly occurs in the shallow parts of the slope. The sliding surface is approximately parallel to the slope surface [42]. The calculation model of the infinite root slope is shown in Figure 2.Figure 2
Schematic diagram of the slope limit equilibrium analysis.According to the shear strength theory of unsaturated soils, the shear strength of unsaturated reinforced soil is obtained to modify the Mohr-Coulomb instability criterion of the infinite slope:(26)τf=cr+c′+σN−σstanϕ′,where cr and c′ are the increasing shear strength caused by the roots and the effective cohesion of the rootless soil; ϕ′ is the effective internal friction angle; σN is the net normal stress; and σs is the suction stress of the VS.The shear strength of the vegetated soil should take into account the role of the roots to strengthen the soil. According to the empirical formula,cr is obtained as
(27)cr=0.48⋅Tr⋅RAR,where Tris the average tensile strength of the roots and RAR is the root area ratio.Normal stress can be expressed as(28)σN=γH−z+γw∫zHθdzcos2α,where γs is the dry weight of the soil and θ is the volumetric capacity of water in the soil.According to the theory of suction stress proposed by Lu et al. [43], the suction stress of soil with VS is
(29)σs=−ua−uw1+δua−uwmm−1/m=−γwψa−γwψw1+δγwψa−ψwmm−1/m,where δ is the reciprocal of the intake pressure and ψa and ψw are the pressure head and pore pressure head, respectively. When the soil is not saturated, the groundwater pressure head is negative: ψa is ψ, ψw is zero. When the soil is saturated, the groundwater pressure head is positive: ψa is zero, ψw is H. The pressure head ψ can be obtained from Equation (1).The FOS of the VS is closely related to the shear strength. According to the limit equilibrium analysis method, FOS is as follows:(30)FS=0.48⋅Tr⋅RAR+c′γH−z+γw∫zHθdzsinαcosα+γwψa−ψwtanϕ′1+δγwψa−ψwmm−1/mγH−z+γw∫zHθdzsinαcosα+tanϕ′tanα.


## 4. Establishment of Numerical Model

### 4.1. Calculation Program

In recent years, many researchers have developed the unsaturated soil seepage theory proposed by Richards. These theories mainly perform coupled flow deformation analysis and capture the interaction between the soil, water, and air. The boundary and initial conditions are different. The analytical method is very difficult to solve the Richards equation. In this paper, the finite element method and the large-scale finite element software ABAQUS are used to analyze the rainfall seepage problem on the VS. For unsaturated soil, the stress field acts on the soil and changes the seepage characteristics of the soil, which causes changes in the head pressure. This paper uses the effective stress equation proposed by Bishop to analyze the seepage behaviour of unsaturated soils.(31)σij′=σij−ua+χua−uwδij,where σij′ is the effective stress between grains, σij is the total stress, ua−uw is the matrix suction, χ is the material property. Khalili et al. [44] obtained the relationship between χ and ψe/ψ through a large number of experimental fittings. The expression can be as follows:
(32)χ=ψψe−0.55,ψ>ψe,1,ψ≤ψe,where ψe is the intake pressure.The modification of the effective stress equation proposed by Cabarkapa et al. [45] and the effective stress equation of unsaturated soil under water absorption hardening and rheology are considered; the effective stress equation is as follows:
(33)G∂2uie∂xi∂xj+λ+G∂2uje∂xi∂xj+χ∂uw∂xi−∂ua∂xi+∂ua∂xi+Fi=0,cm1−χ−ϕa∂uw∂t−∂ua∂t−χ∂∂t∂ui∂x+∂ui∂z=1γw∂∂xKxψ∂ψ∂x+∂∂zKzψ∂ψ∂z−∂Kzψ∂z−SRx,z,ϕaP∂ua∂t+cm1−χ−ϕa∂ua∂t−∂uw∂t−1−χ∂∂t∂ui∂x+∂ui∂z=1γa∂∂xKxψ∂ψ∂x+∂∂zKzψ∂ψ∂z−∂Ka∂z,where G and λ are the Lame function, Fi is the per unit volume force, cm is the compression coefficient, Kw and Ka are the water permeability coefficient and the air permeability coefficient, and P is the absolute pressure.In the finite element method, the stiffness matrix connects the stress and strain at the element nodes. The element stiffness matrix of the unsaturated soil in the roots is as follows:(34)K=Eα1Cα2Cα1CT−θtH1−a11Ma12Mα2CTa21M−θtH2−a22M,where E is the traditional stiffness matrix, C is the coupling between flow and deformation as a function matrix of root water absorption, M is the traditional quality matrix, and a11 and a22 are the water compression coefficient and the air compression coefficient.According to the rainfall boundary conditions, we first assume that the effective rainfall intensity is equal to the permeability coefficient of the root zone. All infiltration at the beginning of the rainfall will not cause stagnant water. The above pore water pressure is substituted into Equation (30) to solve FOS. The pressure head on the slope surface will be greater than zero as the rainfall duration increases. The slope surface is set as the pressure head boundary at this time. The pore water pressure above the groundwater level of the slope was calculated from the new calculation. Because it is difficult to determine the infiltration boundary conditions of the slope in advance, the following formula can be used to judge.
(35)Kxψ∂ψ∂xsinα+Kzψ∂ψ∂z−1cos−TP≤ES.


### 4.2. Model Building

As shown in Figure3, the two-dimensional geometric plane model is established. The roots of the model consider four root architectures distributed in 1.5 deep of the slope. It is assumed that the forces and deformation of the slope are considered two-dimensional plane strain problems. The left and right boundaries of the model are subject to horizontal displacement and pore water pressure constraints. The bottom boundary is rigid. The grid is divided into CPE8P units. In the calculation, the dead weight stress and rainfall surface pore flow are considered. To simulate the exchange of water with the surroundings, the lateral runoff of the slope water during rainfall is considered. A small amount of water (0.06 mm/h) is allowed to be discharged at the foot of the slope.Figure 3
Geometric model of the slope.The yield criterion of slope soil adopts the modified Cambridge model [46]. The shape of the yield surface of the modified Cambridge model is defined by the following formula:
(36)f=q2−M2pp−pc,where M is the slope of the critical state line, pc′ is the preconsolidation pressure, p is the average stress, and q is the generalized shear stress.The multiple hardening reactions of soil-root composites are in the formula for increasing the preconsolidation pressure. This increment consists of three parts: plastic strain increment control, saturation control, and root plastic deformation control. Each part is responsible for a different hardening mechanism.(37)p′¯c=pc′εvp1+eλ−κ−bpc′Sr+Rpp′eεr,where κ is the degree of elastic compression, λ is the degree of plastic compression, b is a partial saturation parameter, and Rp is a constitutive parameter.This numerical analysis is divided into three analysis steps. The first analysis step is to balance the slope ground stress under gravity. The second analysis step is to increase the matrix suction of the slope by the root water absorbed in a steady dry climate. The last step is rainfall.The soil of the model is clay. The physical and mechanical parameters are as follows: elastic modulusE is 50 MPa, dry density γ is 18.6 kg.m-3, saturated volume water content θs is 0.50, residual volume water content θr is 0.15, soil desaturation coefficient β is 0.85, reciprocal of the intake pressure δ is 0.3 kPa-1, e is 1, cohesion c′ is 15 kPa, internal friction angle φ′ is 25°, saturated permeability coefficient ks is 0.030 m·h-1, material parameter m is 5, root tensile strength Tr is 25 kPa, root area ratio RAR is 0.06, PET is 0.24 mm·h-1, and depth of root zone h1 is 1.5 m.To check the grid independence, four different grid systems with total elements of 1,785, 2,035, 3,357, and 3,415 were verified under the same conditions, as shown in Figure4. The root architecture is uniform. The rainfall amount is 400 mm. The time of rainfall is 60 h. From a comparison of FOS profiles, it was found that there was less than 1% deviation among these four cases, indicating that the numerical results used in the model were independent of the grid. Besides, the numerical results showed that further increasing the grid number will not cause obvious variation in FOS when the grid number was 3,415. Thereby, the grid number of 3,415 was selected as the basic grid system to reduce the computational cost.Figure 4
Comparison of FOS with different grids.


### 4.3. Validation

To validate the reliability of the numerical procedures, the numerical results were compared with the experimental results with uniformly distributed roots after an hour of rainfall [27]. Vetiver grass and clay slope was used in the experiment. The slope degree is 30°, c′ is 20 kPa, internal friction angle φ′ is 24°, saturated permeability coefficient ks is 0.031 m·h-1, saturated volume water content θs is 0.45, residual volume water content θr is 0.05, and soil desaturation coefficient β is 0.8. The numerical and experimental results were found to be in close agreement, as shown in Figure 5. The differences between them can be attributed to several aspects, such as root distribution, environmental factors, and the interception of early vegetation canopy, which are ignored in the numerical models. In addition, the differences between the fitted values and measured values of the soil and root-related material properties will also cause certain deviations. It is worth noting that, for VSs, the maximum FOS for rootless areas generally does not exceed 20%. Therefore, the overall difference is acceptable. It can be concluded that this value.Figure 5
Comparison between numerical results and experimental results.


## 4.1. Calculation Program

In recent years, many researchers have developed the unsaturated soil seepage theory proposed by Richards. These theories mainly perform coupled flow deformation analysis and capture the interaction between the soil, water, and air. The boundary and initial conditions are different. The analytical method is very difficult to solve the Richards equation. In this paper, the finite element method and the large-scale finite element software ABAQUS are used to analyze the rainfall seepage problem on the VS. For unsaturated soil, the stress field acts on the soil and changes the seepage characteristics of the soil, which causes changes in the head pressure. This paper uses the effective stress equation proposed by Bishop to analyze the seepage behaviour of unsaturated soils.(31)σij′=σij−ua+χua−uwδij,where σij′ is the effective stress between grains, σij is the total stress, ua−uw is the matrix suction, χ is the material property. Khalili et al. [44] obtained the relationship between χ and ψe/ψ through a large number of experimental fittings. The expression can be as follows:
(32)χ=ψψe−0.55,ψ>ψe,1,ψ≤ψe,where ψe is the intake pressure.The modification of the effective stress equation proposed by Cabarkapa et al. [45] and the effective stress equation of unsaturated soil under water absorption hardening and rheology are considered; the effective stress equation is as follows:
(33)G∂2uie∂xi∂xj+λ+G∂2uje∂xi∂xj+χ∂uw∂xi−∂ua∂xi+∂ua∂xi+Fi=0,cm1−χ−ϕa∂uw∂t−∂ua∂t−χ∂∂t∂ui∂x+∂ui∂z=1γw∂∂xKxψ∂ψ∂x+∂∂zKzψ∂ψ∂z−∂Kzψ∂z−SRx,z,ϕaP∂ua∂t+cm1−χ−ϕa∂ua∂t−∂uw∂t−1−χ∂∂t∂ui∂x+∂ui∂z=1γa∂∂xKxψ∂ψ∂x+∂∂zKzψ∂ψ∂z−∂Ka∂z,where G and λ are the Lame function, Fi is the per unit volume force, cm is the compression coefficient, Kw and Ka are the water permeability coefficient and the air permeability coefficient, and P is the absolute pressure.In the finite element method, the stiffness matrix connects the stress and strain at the element nodes. The element stiffness matrix of the unsaturated soil in the roots is as follows:(34)K=Eα1Cα2Cα1CT−θtH1−a11Ma12Mα2CTa21M−θtH2−a22M,where E is the traditional stiffness matrix, C is the coupling between flow and deformation as a function matrix of root water absorption, M is the traditional quality matrix, and a11 and a22 are the water compression coefficient and the air compression coefficient.According to the rainfall boundary conditions, we first assume that the effective rainfall intensity is equal to the permeability coefficient of the root zone. All infiltration at the beginning of the rainfall will not cause stagnant water. The above pore water pressure is substituted into Equation (30) to solve FOS. The pressure head on the slope surface will be greater than zero as the rainfall duration increases. The slope surface is set as the pressure head boundary at this time. The pore water pressure above the groundwater level of the slope was calculated from the new calculation. Because it is difficult to determine the infiltration boundary conditions of the slope in advance, the following formula can be used to judge.
(35)Kxψ∂ψ∂xsinα+Kzψ∂ψ∂z−1cos−TP≤ES.


## 4.2. Model Building

As shown in Figure3, the two-dimensional geometric plane model is established. The roots of the model consider four root architectures distributed in 1.5 deep of the slope. It is assumed that the forces and deformation of the slope are considered two-dimensional plane strain problems. The left and right boundaries of the model are subject to horizontal displacement and pore water pressure constraints. The bottom boundary is rigid. The grid is divided into CPE8P units. In the calculation, the dead weight stress and rainfall surface pore flow are considered. To simulate the exchange of water with the surroundings, the lateral runoff of the slope water during rainfall is considered. A small amount of water (0.06 mm/h) is allowed to be discharged at the foot of the slope.Figure 3
Geometric model of the slope.The yield criterion of slope soil adopts the modified Cambridge model [46]. The shape of the yield surface of the modified Cambridge model is defined by the following formula:
(36)f=q2−M2pp−pc,where M is the slope of the critical state line, pc′ is the preconsolidation pressure, p is the average stress, and q is the generalized shear stress.The multiple hardening reactions of soil-root composites are in the formula for increasing the preconsolidation pressure. This increment consists of three parts: plastic strain increment control, saturation control, and root plastic deformation control. Each part is responsible for a different hardening mechanism.(37)p′¯c=pc′εvp1+eλ−κ−bpc′Sr+Rpp′eεr,where κ is the degree of elastic compression, λ is the degree of plastic compression, b is a partial saturation parameter, and Rp is a constitutive parameter.This numerical analysis is divided into three analysis steps. The first analysis step is to balance the slope ground stress under gravity. The second analysis step is to increase the matrix suction of the slope by the root water absorbed in a steady dry climate. The last step is rainfall.The soil of the model is clay. The physical and mechanical parameters are as follows: elastic modulusE is 50 MPa, dry density γ is 18.6 kg.m-3, saturated volume water content θs is 0.50, residual volume water content θr is 0.15, soil desaturation coefficient β is 0.85, reciprocal of the intake pressure δ is 0.3 kPa-1, e is 1, cohesion c′ is 15 kPa, internal friction angle φ′ is 25°, saturated permeability coefficient ks is 0.030 m·h-1, material parameter m is 5, root tensile strength Tr is 25 kPa, root area ratio RAR is 0.06, PET is 0.24 mm·h-1, and depth of root zone h1 is 1.5 m.To check the grid independence, four different grid systems with total elements of 1,785, 2,035, 3,357, and 3,415 were verified under the same conditions, as shown in Figure4. The root architecture is uniform. The rainfall amount is 400 mm. The time of rainfall is 60 h. From a comparison of FOS profiles, it was found that there was less than 1% deviation among these four cases, indicating that the numerical results used in the model were independent of the grid. Besides, the numerical results showed that further increasing the grid number will not cause obvious variation in FOS when the grid number was 3,415. Thereby, the grid number of 3,415 was selected as the basic grid system to reduce the computational cost.Figure 4
Comparison of FOS with different grids.


## 4.3. Validation

To validate the reliability of the numerical procedures, the numerical results were compared with the experimental results with uniformly distributed roots after an hour of rainfall [27]. Vetiver grass and clay slope was used in the experiment. The slope degree is 30°, c′ is 20 kPa, internal friction angle φ′ is 24°, saturated permeability coefficient ks is 0.031 m·h-1, saturated volume water content θs is 0.45, residual volume water content θr is 0.05, and soil desaturation coefficient β is 0.8. The numerical and experimental results were found to be in close agreement, as shown in Figure 5. The differences between them can be attributed to several aspects, such as root distribution, environmental factors, and the interception of early vegetation canopy, which are ignored in the numerical models. In addition, the differences between the fitted values and measured values of the soil and root-related material properties will also cause certain deviations. It is worth noting that, for VSs, the maximum FOS for rootless areas generally does not exceed 20%. Therefore, the overall difference is acceptable. It can be concluded that this value.Figure 5
Comparison between numerical results and experimental results.


## 5. Results and Analysis

### 5.1. Effect of Rainfall Duration

From Figure6, it can be seen that the FOS decreases with increasing depth in the root zone at the initial state, and FOS does not change with increasing depth in the rootless zone. FOS is abruptly changed for uniform roots and oval roots at the boundary between the root zone and the rootless area. This phenomenon is because the roots are distributed differently at the root zone and rootless zone boundary. After drying for 12 h, the FOS of the root zone increased significantly and the FOS of the rootless area increased near the interface, indicating that roots transpiration can affect the rootless zone. From the FOS curves, it can be seen that there is no linear correlation between the FOS and the rainfall duration. The FOS of the four root architectures with a rainfall duration of 12 h are greater than 1.0, which is safe; the FOS of the slope depth with a rainfall duration of 24 h is less than 1.0 within a depth of 2.5 m, which is unsafe; when the rainfall duration is 36 h, the FOS of the uniform slope are less than 1.0 at a depth of 0.25-1.1 m and 1.5-2 m, the FOS of the exponential root slope are less than 1.0 at a depth of 0.27-1.9 m, the FOS of the triangular root slope are less than 1.0 at a depth of 0.2-2.1 m, and the FOS of oval root slope less than 1.0 at a depth of 0.2-1.7 m, which is usually shallow landslides. With the increase of depth, the FOS of the four root architectures showed a trend of decreasing first and then increasing, and FOS was the smallest at the root zone and rootless zone boundary (1.5 m), indicating that this boundary was most prone to instability. The FOS increases with increasing rainfall duration when the rainfall duration is more than 24 h.Figure 6
The FOS of different root architectures for rainfall amount of 500 mm with different depths: (a) uniformly, (b) exponential, (c) triangular, and (d) oval.
(a)(b)(c)(d)


### 5.2. Effect of Rainfall Amount

Figure7 shows the FOS of the uniform roots for a depth of 1.5 m with a total rainfall of 400 mm, 500 mm, and 600 mm. It can be seen that FOS increases with the time of drying for 12 h in a stable, indicating that the water absorption of roots leads to an increase in the suction of the substrate, which increases the shear strength of the slope. For rainfall, the FOS decay rate increases with decreasing rainfall duration under the same rainfall conditions. The FOS increases first and then decreases with the increasing of rainfall duration. The FOS in each period is greater than 1.0, and the shallow slope is stable when the rainfall is 400 mm; the rainfall duration with FOS less than 1.0 are 24 h, 36 h, and 48 h when the rainfall is 500 mm. The instability time is 46 h, 42 h, 22 h, and 12 h occurred for rainfall duration of 24 h, 36 h, 48 h, and 60 h curves for the rainfall amount of 600 mm. The results show that with the same rainfall duration, the greater the rainfall intensity, the more likely the shallow instability of the VS.Figure 7
The FOS of uniform roots slope with the depth of 1.5 m with different rainfall durations: (a) 400 mm, (b) 500 mm, and (c) 600 mm.
(a)(b)(c)Figure8 shows the pressure head of the uniform roots for the rainfall duration of 60 h with a rainfall amount of 500 mm. It can be seen that the pressure head of the shallow layer decreased and the matric suction increased after 12 h drying, which indicated that the matric suction was significantly increased in the presence of roots under drying conditions. The pressure head of the shallow layer increased and the matric suction decreased relative to the rainfall duration. The underground water table increased after 12 h of rainfall, indicating that the FOS of the slope decreases. When there is saturated zone runoff, the groundwater table decreased slowly after 60 h of rainfall, indicating that the FOS of the slope increases.Figure 8
Pressure head distributed on the uniform roots for rainfall duration of 60 h with rainfall amount of 500 mm: (a) drying 12 h, (b) rainfall 12, (c) rainfall 36, and (d) rainfall 60 h.
(a)(b)(c)(d)


### 5.3. Effect of Runoff and Transpiration

To reflect the relationship between the boundary conditions and the root hydraulic effect, the cumulative runoff and cumulative evaporation of 500 mm rainfall at different times on the uniform VS were obtained, as shown in Figure9. It can be seen from Figure 9(a) that the runoff is getting larger and larger with increasing rainfall time and the accumulated water is increasing. The runoff is getting smaller and smaller when the rainfall is over. The runoff growth and stagnant water increase with increasing rainfall intensity in the early stage. As shown in Figure 9(b), the cumulative evaporation slowly increases during the drying period. The shorter the rainfall duration, the faster the cumulative evaporation increases, the greater the substrate suction, and the faster the FOS increases. The possibility of VS failure is slight under the condition of large rainfall intensity. At the same rainfall, runoff increases with the increase of rainfall duration. The increase of dry evaporation accelerates the infiltration of the accumulated water at the end of the rainfall period. During prolonged rainfall periods, water flowing from the slope boundary will not trigger slope damage.Figure 9
Cumulative runoff and cumulative evaporation for the rainfall amount of 500 mm at different times: (a) cumulative runoff and (b) cumulative evaporation.
(a)(b)


### 5.4. Effect of Vegetation Root Architecture

Figure10 shows the FOS of four root architectures at 500 mm rainfall. It can be seen that the rainfall intensity is higher, and the decline rate of FOS is faster, but FOS changes with different root architectures. The time that the oval roots induced shallow slope instability is 29 h during rainfall duration 24 h and 33 h during the rainfall duration 36 h. The time that the exponential roots induced shallow slope instability is 23 h during rainfall duration 24 h and 31 h during rainfall duration 36 h, and 51 h during the rainfall duration 48 h. The time that the triangle roots induced shallow slope instability is 31 h during the rainfall duration 24 h and 35 h during rainfall duration 36 h. Only the FOS of oval roots slope is greater than the one when the rainfall duration is 48 h. The results in this paragraph indicate that the root architecture also affects the stability of the shallow slope.Figure 10
FOS of the VS for the rainfall amount of 500 mm with different root architectures: (a) oval, (b) exponential, and (c) triangle.
(a)(b)(c)


### 5.5. Landslide Threshold

According to the previous results, it can be known that the root architecture, rainfall amount, and rainfall duration may affect FOS. Therefore, we have considered the above three factors and obtained four landslide threshold curves. Based on the four root architectures, landslide thresholds are established based on different landslide amount and rainfall duration as shown in Figure11. As can be seen from Figure 11, the rainfall threshold of the four root architectures divides the calculation area into a safe area and an unsafe area. If rainfall and duration are unsafe areas to the right of the threshold curve, the VS is unstable. Conversely, if the rainfall is a safe area to the left of the threshold curve, VS is stable. Taking 500 mm and 48 h rainfall as an example, landslides are induced in the triangular and exponential roots slopes, while the same rainfall and duration do not cause landslides in oval and uniform roots.Figure 11
Landslide threshold curve of rainfall-induced slope instability of four root architectures.It can be known from Figure11 that the landslide threshold of the four root architectures includes the upper and lower boundaries. For example, the rainfall duration of the lower bound is 12.5 h and the rainfall duration of the upper bound is 43 h when the shallow instability of the uniform roots, indicating that the slope is unsafe and prone to landslides when the rainfall lasts from 12.5 h to 43 h. For a certain rainfall amount, the heavy rainfall amount with the rainfall duration below the lower bound and the weak rainfall amount with the rainfall duration above the upper limit will not cause the shallow instability. It is due to the effect of boundary runoff and transpiration suction of the roots. The matric suction increases with decreasing the pressure head. It is shown that the pore water pressure dissipation caused by the effect of boundary runoff on the boundary and the transpiration of roots can increase the FOS. Unlike one-dimensional simulation results provided by Wu et al. [40], VS always slips when the duration is greater than the lower bound.


## 5.1. Effect of Rainfall Duration

From Figure6, it can be seen that the FOS decreases with increasing depth in the root zone at the initial state, and FOS does not change with increasing depth in the rootless zone. FOS is abruptly changed for uniform roots and oval roots at the boundary between the root zone and the rootless area. This phenomenon is because the roots are distributed differently at the root zone and rootless zone boundary. After drying for 12 h, the FOS of the root zone increased significantly and the FOS of the rootless area increased near the interface, indicating that roots transpiration can affect the rootless zone. From the FOS curves, it can be seen that there is no linear correlation between the FOS and the rainfall duration. The FOS of the four root architectures with a rainfall duration of 12 h are greater than 1.0, which is safe; the FOS of the slope depth with a rainfall duration of 24 h is less than 1.0 within a depth of 2.5 m, which is unsafe; when the rainfall duration is 36 h, the FOS of the uniform slope are less than 1.0 at a depth of 0.25-1.1 m and 1.5-2 m, the FOS of the exponential root slope are less than 1.0 at a depth of 0.27-1.9 m, the FOS of the triangular root slope are less than 1.0 at a depth of 0.2-2.1 m, and the FOS of oval root slope less than 1.0 at a depth of 0.2-1.7 m, which is usually shallow landslides. With the increase of depth, the FOS of the four root architectures showed a trend of decreasing first and then increasing, and FOS was the smallest at the root zone and rootless zone boundary (1.5 m), indicating that this boundary was most prone to instability. The FOS increases with increasing rainfall duration when the rainfall duration is more than 24 h.Figure 6
The FOS of different root architectures for rainfall amount of 500 mm with different depths: (a) uniformly, (b) exponential, (c) triangular, and (d) oval.
(a)(b)(c)(d)


## 5.2. Effect of Rainfall Amount

Figure7 shows the FOS of the uniform roots for a depth of 1.5 m with a total rainfall of 400 mm, 500 mm, and 600 mm. It can be seen that FOS increases with the time of drying for 12 h in a stable, indicating that the water absorption of roots leads to an increase in the suction of the substrate, which increases the shear strength of the slope. For rainfall, the FOS decay rate increases with decreasing rainfall duration under the same rainfall conditions. The FOS increases first and then decreases with the increasing of rainfall duration. The FOS in each period is greater than 1.0, and the shallow slope is stable when the rainfall is 400 mm; the rainfall duration with FOS less than 1.0 are 24 h, 36 h, and 48 h when the rainfall is 500 mm. The instability time is 46 h, 42 h, 22 h, and 12 h occurred for rainfall duration of 24 h, 36 h, 48 h, and 60 h curves for the rainfall amount of 600 mm. The results show that with the same rainfall duration, the greater the rainfall intensity, the more likely the shallow instability of the VS.Figure 7
The FOS of uniform roots slope with the depth of 1.5 m with different rainfall durations: (a) 400 mm, (b) 500 mm, and (c) 600 mm.
(a)(b)(c)Figure8 shows the pressure head of the uniform roots for the rainfall duration of 60 h with a rainfall amount of 500 mm. It can be seen that the pressure head of the shallow layer decreased and the matric suction increased after 12 h drying, which indicated that the matric suction was significantly increased in the presence of roots under drying conditions. The pressure head of the shallow layer increased and the matric suction decreased relative to the rainfall duration. The underground water table increased after 12 h of rainfall, indicating that the FOS of the slope decreases. When there is saturated zone runoff, the groundwater table decreased slowly after 60 h of rainfall, indicating that the FOS of the slope increases.Figure 8
Pressure head distributed on the uniform roots for rainfall duration of 60 h with rainfall amount of 500 mm: (a) drying 12 h, (b) rainfall 12, (c) rainfall 36, and (d) rainfall 60 h.
(a)(b)(c)(d)


## 5.3. Effect of Runoff and Transpiration

To reflect the relationship between the boundary conditions and the root hydraulic effect, the cumulative runoff and cumulative evaporation of 500 mm rainfall at different times on the uniform VS were obtained, as shown in Figure9. It can be seen from Figure 9(a) that the runoff is getting larger and larger with increasing rainfall time and the accumulated water is increasing. The runoff is getting smaller and smaller when the rainfall is over. The runoff growth and stagnant water increase with increasing rainfall intensity in the early stage. As shown in Figure 9(b), the cumulative evaporation slowly increases during the drying period. The shorter the rainfall duration, the faster the cumulative evaporation increases, the greater the substrate suction, and the faster the FOS increases. The possibility of VS failure is slight under the condition of large rainfall intensity. At the same rainfall, runoff increases with the increase of rainfall duration. The increase of dry evaporation accelerates the infiltration of the accumulated water at the end of the rainfall period. During prolonged rainfall periods, water flowing from the slope boundary will not trigger slope damage.Figure 9
Cumulative runoff and cumulative evaporation for the rainfall amount of 500 mm at different times: (a) cumulative runoff and (b) cumulative evaporation.
(a)(b)


## 5.4. Effect of Vegetation Root Architecture

Figure10 shows the FOS of four root architectures at 500 mm rainfall. It can be seen that the rainfall intensity is higher, and the decline rate of FOS is faster, but FOS changes with different root architectures. The time that the oval roots induced shallow slope instability is 29 h during rainfall duration 24 h and 33 h during the rainfall duration 36 h. The time that the exponential roots induced shallow slope instability is 23 h during rainfall duration 24 h and 31 h during rainfall duration 36 h, and 51 h during the rainfall duration 48 h. The time that the triangle roots induced shallow slope instability is 31 h during the rainfall duration 24 h and 35 h during rainfall duration 36 h. Only the FOS of oval roots slope is greater than the one when the rainfall duration is 48 h. The results in this paragraph indicate that the root architecture also affects the stability of the shallow slope.Figure 10
FOS of the VS for the rainfall amount of 500 mm with different root architectures: (a) oval, (b) exponential, and (c) triangle.
(a)(b)(c)


## 5.5. Landslide Threshold

According to the previous results, it can be known that the root architecture, rainfall amount, and rainfall duration may affect FOS. Therefore, we have considered the above three factors and obtained four landslide threshold curves. Based on the four root architectures, landslide thresholds are established based on different landslide amount and rainfall duration as shown in Figure11. As can be seen from Figure 11, the rainfall threshold of the four root architectures divides the calculation area into a safe area and an unsafe area. If rainfall and duration are unsafe areas to the right of the threshold curve, the VS is unstable. Conversely, if the rainfall is a safe area to the left of the threshold curve, VS is stable. Taking 500 mm and 48 h rainfall as an example, landslides are induced in the triangular and exponential roots slopes, while the same rainfall and duration do not cause landslides in oval and uniform roots.Figure 11
Landslide threshold curve of rainfall-induced slope instability of four root architectures.It can be known from Figure11 that the landslide threshold of the four root architectures includes the upper and lower boundaries. For example, the rainfall duration of the lower bound is 12.5 h and the rainfall duration of the upper bound is 43 h when the shallow instability of the uniform roots, indicating that the slope is unsafe and prone to landslides when the rainfall lasts from 12.5 h to 43 h. For a certain rainfall amount, the heavy rainfall amount with the rainfall duration below the lower bound and the weak rainfall amount with the rainfall duration above the upper limit will not cause the shallow instability. It is due to the effect of boundary runoff and transpiration suction of the roots. The matric suction increases with decreasing the pressure head. It is shown that the pore water pressure dissipation caused by the effect of boundary runoff on the boundary and the transpiration of roots can increase the FOS. Unlike one-dimensional simulation results provided by Wu et al. [40], VS always slips when the duration is greater than the lower bound.


## 6. Discussion

It is a challenging task to predict or control landslides by using suitable measures. The vegetation roots can increase the shallow stability of the VS, mainly due to the hydraulic characteristics of vegetation [47]. Root water uptake, vegetation interception, and root increment have effects on soil strength at time scale and seasons [48, 49]. The results indicate that vegetation type affects the landslide threshold of slope failure, which is consistent with the root water uptake of different vegetation. The shade-tolerant tree species had higher root specific hydraulic conductivity, and the difference of hydraulic conductivity among different root rings was small [50, 51]. Different root architectures have different water absorption, water holding capacity, and the ability to improve matric suction, so the rainfall threshold is different. The rainfall interception by vegetation canopies plays an important role in slope stabilization, which is related to the vegetation coverage [52]. The interception is not involved in this paper, and further research in this area should be strengthened in the future.Previous studies showed that vegetation can control soil erosion [53]. The results indicate that vegetation type affects the runoff, which is consistent with the soil erosion of different vegetation [54, 55]. The Vicia sativa L. reduced runoff discharge from 32.87 to 13.68% [56]. Vegetation can be a sustainable solution to reduce soil runoff and erosion [57, 58]. The positive effects of vegetation on rainfall in the shallow slope are shown in two aspects: the reduction of slope runoff and sediment transport [59, 60]; the second is that vegetation coverage gradually decreases with the rate of soil erosion within a certain number of years [61]. The runoff depends on the topographic factor. The larger plot has a smaller runoff per unit area, and the runoff coefficient decreases with the decrease of slope degree [62]. Therefore, the amount of runoff from soil erosion, which depends on a variety of factors, must be considered together. In addition to different root architectures, different conditions (vegetation coverage, rainfall intensity, slope degree, and soil texture) have effects on slope runoff and sediment yield. The amount of runoff and sediment reduction increased with the increase of vegetation coverage, and the effect of vegetation on sediment reduction was greater than that of runoff control under heavy rainfall. Sediment reduction effect of grassland is generally better than that of shrub and tree vegetation types [63].


## 7. Conclusion

In this paper, the two-dimensional Richards equation and the method of limit equilibrium slope stabilization are used to simulate the shallow instability of VS. Based on this, a numerical model is developed to analyze the shallow stability of VS. The method is used to analyze the effects of rainfall and root architecture on the shallow slope stability. The main results are as follows:(1)
FOS increases first and then decreases with the increasing of rainfall duration and depth. The interface between the root zone and the rootless zone has the smallest FOS and is most prone to landslides(2)
The rainfall runoff increases with the increasing of rainfall duration. The increase of dry evaporation accelerates the infiltration of the accumulated water, the probability of shallow landside increases as the rainfall amount increases, or the rainfall duration decreases(3)
Slope boundary water outflow and roots water absorption may have a positive effect on stability of shallow slopes. Threshold curves for four root architectures were established, revealing the upper and lower bounds of rainfall duration(4)
The arrangement of root architecture in safe areas from large to small is uniform roots, oval roots, triangular roots, and exponential roots. The time of the safe area is 35 h for rainfall amount 500 m in the uniformly roots slope



---

*Source: 6682113-2021-06-09.xml*