Influence of Internal Erosion on Rainfall-Induced Instability of Layered Deposited-Soil Slopes

Abstract

Layered deposited-soil slopes are widely distributed in mountainous terrain. The rainfall-induced instability of layered deposited-soil slopes is not only controlled by the unsaturated infiltration process but also by the seepage-induced internal-erosion process within the deposited soils. In this paper, the main physical processes within two-layer deposited-soil slopes under rainfall infiltration are summarized, and a coupled seepage–erosion finite element model is established to analyze the interactions between the rainfall infiltration process and the internal-erosion process within layered deposited-soil slopes. This finite element model was validated by simulating the coupled seepage–erosion process in a one-dimensional layered soil column. Then, serials of two-dimensional coupled seepage–erosion simulations were conducted to investigate the rainfall-induced seepage–erosion patterns, as well as their impact on the stability evolution of the layered deposited slopes. It was shown that the rainfall-induced seepage–erosion accelerate the water infiltration into the slope and facilitate the generation of subsurface stormflow near the layer interface, which will weaken the soils around the layer interface and accelerate the slope failure process inevitably. Special attention should be paid to the rainfall-induced seepage–erosion effect when evaluating the stability of layered deposited-soil slopes.

1. Introduction

Geological disasters occur frequently in mountainous regions with complex terrains and geological conditions, such as the Alps and Himalaya areas [1,2]. In addition to the direct damage caused by the geological disasters themselves, these disasters also lead to the accumulation of large amounts of loose rock and soil debris deposited on the slope surfaces, which provide favorable loose-material source conditions for the development of secondary geological disasters, such as rainfall-induced landslides and debris flows [3,4,5]. Compared with the well-consolidated soils within the original slope, the soils deposited on the slope surface generally belong to a kind of wide-graded, weakly consolidated soil, usually characterized by a loose structure, high porosity, and low vegetation coverage [6,7]. Under rainfall conditions, these deposited soils are highly susceptible to seepage-induced internal erosion [8,9], which further triggers instability in the layered deposited-soil slopes. The disasters caused by the instability of these layered deposited-soil slopes have a wide impact range and could last for very long durations [10,11], making it one of the main forms of secondary geological disasters in mountainous areas.

In recent years, a significant amount of research has been undertaken to investigate the hydro-mechanical response of unsaturated layered soils during rainfall infiltration. For example, Yang et al. (2006) [12] performed serials of infiltration tests in a soil column apparatus to evaluate the effect of the rainfall intensity and duration on the infiltration in a layered soil column; Wu et al. (2012 [13], 2016 [14]) conducted both analytical and numerical analyses of the 1D coupled water infiltration and deformation in two-layer unsaturated soils; Tsai et al. (2013) [15] investigated the importance of the soil layer distribution for rainfall-induced shallow landslides by extending a physical-based model to consider the layered infinite slope; Gan and Zhang (2020) [16] performed physical-model tests as well as hydro-mechanically coupled finite element modeling to investigate the rainfall-induced failure behavior of layered slopes; Girardi et al. (2023) [17] also investigated the failure and post-failure behavior of layered slopes under rainfall infiltration with the material-point method. According to these studies, due to the distinct hydro-mechanical properties of soils within different layers, different soil layer distributions could lead to distinct seepage saturation paths and instability failure modes; for the classic two-layer slopes with higher permeable upper layers, significant positive pore-water pressure can accumulate above the soil layer interface due to the delayed drain-out of infiltrated water into deeper layers, which could form a potential weak layer for slope failure along this layer interface [18,19].

Although the seepage process and failure modes of layered soil slopes upon rainfall infiltration have been well studied in the literature, there are no reports in the research on the associated seepage-induced internal-erosion process within layered slopes, which is prone to take place within loose deposited-soil layers. Accordingly, the impact of internal erosion on the rainfall infiltration process, as well as on the stability of layered deposited slopes, remains unclear. Therefore, it was the object of this research to investigate the rainfall-induced internal-erosion process within layered deposited slopes quantitatively with a proper physical-based mathematical model, so as to provide a better understanding of the intrinsic failure mechanism and more comprehensive safety evaluation approach for layered deposited-soil slopes under rainfall infiltration. Toward this goal, the main physical processes within the classic two-layer deposited-soil slopes under rainfall infiltration are summarized firstly. Then, a mathematical framework based on the mixture theory is presented for capturing the seepage-induced internal-erosion process within unsaturated deposited soils; within this framework, the rainfall-induced seepage–erosion process in a one-dimensional layered soil column, including the erosion, transport, and deposition of fine soil particles, was investigated. Further, serials of two-dimensional coupled seepage–erosion finite element simulations were conducted to investigate the rainfall-induced seepage–erosion patterns, as well as their impact on the stability evolution of layered deposited slopes.

2. Physical Processes and Mathematical Modeling Framework

2.1. Overview of the Main Physical Processes

The loose deposited soils formed by earth surface processes, such as landslides and debris flows, are generally wide-grading and poorly consolidated [6,7]. Due to their relatively short deposition time and loose soil structure, the fine particles (fines) within these deposited soils are highly prone to migrate within the pores formed by the coarse particles when subjected to interstitial water flow [20,21]. In contrast, the beneath-soil layers covered by the deposited soils are usually well consolidated and less permeable than the loose deposited soils [22]. The main physical processes within classic deposited-soil slopes with two-layer soil structures are schematically illustrated in Figure 1. Upon rainfall, water infiltrates into the unsaturated soils gradually and flows from the upper and surface parts of the slope toward the deeper and toe parts of the slope under gravity [23,24] (Figure 1a). Because the deposited soils are prone to internal erosion, fine soil particles within the upper layer can be detached by the seepage forces and migrate along with the seepage water as its suspension [25] (Figure 1b). Near the slope surface, the rainfall-induced fines migration leads to the so-called “surface coarsening” effect [7], in which the loss of fine particles increases the soil porosity and, hence, its permeability, leading to the acceleration of rainfall infiltration. As the underlying layer is usually less permeable than the loose deposited soils, a perched water table is formed above the soil layer interface gradually when the infiltrated water reaches the interface [14,24] (Figure 1c). This causes a rapid rise in the pore-water pressure head (rapid decrease in suction), and a reduction in the shear strength of the partially saturated soils. This combined effect could lead to the failure of the upper layer along the interface. In such a case, a stable seepage path could be formed along the soil layer interface (i.e., subsurface stormflow [26]). What is worse, when the fine particles eroded from the upper layer migrate to the layer interface, they are prone to deposition due to the distinct permeability and pore size characteristics between the two layers [27,28]. This fine deposition process blocks the downward seepage path and could further accelerate the accumulation of perched water. Meanwhile, the stable water seepage along the interface (i.e., subsurface stormflow) will inevitably lead to serious erosion in the loose deposited soils near the interface with time, which could further weaken the stability of the upper soil layer.

2.2. Mathematical Framework

Mathematically, the unsaturated deposited soil can be conceptualized as a three-phase, multi-species, porous medium [29,30], as shown in Figure 2. The porous medium consists of the solid phase ($S$), the liquid phase ($L$), and the gas phase ($G$), with the coarse particle species ($c^{\prime }$), solid fine particle species ($e^{\prime }$), liquid fine particles ($e$), water species ($w$), and gas species ($g$) as their corresponding components. The superscript ‘‘$’$’’ indicates the species that are attached to the solid phase, whereas no superscript is used to refer to the species within the fluid phases. Due to the erosion and deposition processes, the solid fine species ($e^{\prime }$) and fluidized fine species ($e$) can exchange mass with each other.

Assuming all the species are intrinsically incompressible, the total volume of the representative volume element (RVE) is calculated as the sum of the individual volumes occupied by each species, which can also be expressed as the sum of the individual volumes occupied by each phase. Based on the mixture theory, each phase is assumed to occupy the whole RVE volume with its corresponding apparent density (${\rho }^{\alpha }$), which is linked to its corresponding intrinsic density (${\rho }_{\alpha }$) through its respective volume fraction ($n_{\alpha }$) (i.e., ${\rho }^{\alpha }=n_{\alpha }{\rho }_{\alpha }$) (top row of Figure 2b). At the same time, it also can be treated as the superimposition of all the species particles within the RVE with their apparent densities (${\rho }^j$) (bottom row of Figure 2b). The erosion and deposition phenomena of fine particles are treated as the phase-transfer processes between the solid fine species ($e^{\prime }$) and fluidized fine species ($e$) by assuming that their intrinsic densities in both phases are identical. To be consistent with the classic terminology used in soil mechanics, the volume fraction of the whole pore space (liquid fraction plus gas fraction) is referred to as the porosity ($n$), whereas the relative volume fraction occupied by the liquid within these pores is referred to as the liquid saturation ($S_L$). Further, the volumetric concentration ($c_e$) is used to refer to the relative volume fraction of the fluidized fine species ($e$) in the liquid phase, and the mass content ($m_{e^{\prime }}$) is adopted to represent the mass of solid fine species ($e^{\prime }$) per unit RVE volume prior to deformation.

Following the work of Lei (2022) [30], the coupled seepage–erosion process within unsaturated soils can be described by the following set of governing equations, by assuming zero air pressure. The governing equations include the mass conservation equation of the liquid phase (Equation (1)), the mass conservation equation of fluidized fines (Equation (2)), the erosion rate equation (Equation (3)), and the deposition rate equation (Equation (4)). The corresponding four main unknowns are selected as the liquid pore pressure ($p_L$), the volumetric concentration of fluidized fines ($c_e$), the mass content of fines at the pore surface ($m_{er}$), and the mass content of fines at the pore throats ($m_{de}$).

$$\frac{\left(1-S_L\right)}{{\rho }_{e^{\prime }}}\frac{\partial }{\partial t}\left(m_{er}+m_{de}\right)+n\frac{\partial S_L}{\partial p_L}\frac{\partial p_L}{\partial t}+\nabla \cdot J^L=0$$

(1)

$$S_{L}n\frac{\partial }{\partial t}{c}_e+\frac{1-c_e}{{\rho }_{e^{\prime }}}\frac{\partial }{\partial t}\left(m_{er}+m_{de}\right)+J^L\cdot \nabla c_e=0$$

(2)

$$\frac{\partial }{\partial t}\frac{m_{er}}{{\rho }_{e^{\prime }}}=-k_e\cdot \left(x_{e^{\prime }}-x_{e^{\prime }\infty }\right)\cdot \Vert v^{wS}\Vert$$

(3)

$$\frac{\partial }{\partial t}\frac{m_{de}}{{\rho }_{e^{\prime }}}=k_d\cdot c_e\cdot \Vert v^{wS}\Vert$$

(4)

In the above equations, the total solid-fines mass content ($m_{e^{\prime }}$) is divided into two parts with $m_{e^{\prime }}=m_{er}+m_{de}$ (i.e., (1) the solid fines at the pore surface ($m_{er}$), which are available for erosion; (2) the solid fines at the pore throats ($m_{de}$), which refer to the fines deposited at the pore throats). The porosity ($n$) is updated with $n=n_0+\left({\dot{m}}_{e^{\prime }}/{\rho }_{e^{\prime }}\right)\Delta t$ during each step; $J^L$ is the seepage volume flux of the liquid; $\Vert v^{wS}\Vert$ represents the norm of the velocity vector of the pore water relative to the soil skeleton, in which the relative velocity ($v^{wS}$) can be approximated as $v^{wS}=J^L/n$; $k_e$ and $k_d$ are the material constants to reflect a specific soil’s sensitivity for fines erosion and deposition; $x_{e’}$ is the current mass fraction of fines available for erosion (erodible fines), $x_{e’}=m_{er}/\left(m_{e^{\prime }}+m_{c^{\prime }}\right)$; $x_{e’\infty }$ is the final mass fraction of the erodible fine particles under the seepage velocity ($v^{wS}$), which can be expressed by Cividini and Gioda (2004) [31]’s fitting expression:

$$x_{e’\infty }\left(v^{wS}\right)=\{\begin{array}{l}{x}_{e’0}-\left(x_{e^{\prime }0}-x_{e^{\prime }\infty }^{*}\right)\frac{v^{wS}}{v^{*}}\mathrm{if}0<v^{wS}<v^{*}\\ x_{e^{\prime }\infty }^{*}-{\alpha }_{er}\log \left(\frac{v^{wS}}{v^{*}}\right)\mathrm{if}{v}^{wS}>v^{*}\end{array}$$

(5)

where $v^{*}$ is the reference steady-flow velocity of the pore water, under which the erosion process is very slow; $x_{e^{\prime }\infty }^{*}$ is the final mass fraction of fine particles at this reference velocity after an infinite long time; $x_{e’0}$ is the initial mass fraction of fine particles before erosion; ${\alpha }_{er}$ is a parameter controlling the decreasing rate of $x_{e’\infty }$ with the increasing flow velocity ($v^{wS}$).

The term $\partial S_L/\partial p_L$ in Equation (1) accounts for the soil water retention property, which is usually represented by a specific soil water characteristic curve (SWCC). Any expression that can best fit the relationship between the liquid saturation ($S_L$) and matric suction ($p_c$) of the soils under investigation would be suitable. In this paper, the widely used Van Genuchten model [32] is adopted and presented below, as there are ready-to-use fitting parameters in the literature for our subsequent numerical examples:

$$S_L=S_{\min }+\left(S_{\max }-S_{\min }\right){\left[1+{\left({\alpha }_v{p}_c\right)}^{n_v}\right]}^{-\left(1-1/n_v\right)}$$

(6)

where $S_{\min }$ and $S_{\max }$ are the soil’s residual and maximum liquid saturation; $p_c$ is the matric suction, and $p_c=-p_L$ by assuming zero air pressure; ${\alpha }_v$ and $n_v$ are the fitting parameters.

The liquid-seepage flux ($J^L$) is usually calculated using the Darcy law in cases of laminar flow:

$$J^L=-\frac{K_h}{{\gamma }_L}\left(\nabla p_L-{\gamma }_L\right)$$

(7)

where ${\gamma }_L$ is the specific weight of the liquid phase; $K_h$ is the soil hydraulic conductivity, which can be expressed as follows:

$$K_h=K_{h0}{K}_r=K_{h0}\exp \left(-{\alpha }_v{p}_c\right)\frac{{\left[(1-k_r{m}_{de})n/n_0\right]}^{n_r}}{1+2.5c_e}$$

(8)

In the above equation, $K_{h0}$ is the initial hydraulic conductivity at full saturation; $n_r$ and $k_r$ are the material parameters related to the intrinsic permeability changes. The permeability ratio ($K_r$) is introduced to account for the permeability variation. The unsaturated effect is taken into account by Gardner (1958) [33]’s equation, in which a higher suction ($p_c$) will lead to a lower permeability ratio. The combined effects of internal erosion are accounted for via a modified Carman–Kozeny equation [34], in which the erosion effect is taken into account via the porosity change caused by erosion with $n=n_0+\left({\dot{m}}_e/{\rho }_{e^{\prime }}\right)\Delta t$; the pore-fines clogging effect is related to the fines deposited at the pore throats ($m_{de}$); the viscosity of the seepage liquid is also set to depend on the volumetric concentration of liquidized fines ($c_e$).

3. 1D Simulation of Seepage–Erosion Process in Layered Soil Column

The coupled seepage–erosion governing equations proposed in Section 2.2 represent the fundamental mathematical description for the coupled seepage–erosion process within each single soil layer. They were implemented in the open-source finite element code Bil v2.0 (http://perso.lcpc.fr/dangla.patrick/bil/, accessed on 16 October 2023) and validated and verified in our previous works [28]. For the two-layered deposited slopes considered in this paper, the coupled seepage–erosion processes within these two soil layers are described separately via these mathematical equations with two different sets of material parameters. To validate the code’s capability in dealing with unsaturated seepage–erosion among different soil layers, the rainfall infiltration process of a one-dimensional layered soil column, originally analyzed by Wu et al. (2016) [14], was simulated first; further, the rainfall-induced seepage–erosion process in this one-dimensional layered soil column, including the erosion, transport, and deposition of fine particles, was simulated and analyzed in detail.

A simplified model of a layered unsaturated soil column with a height of 2 m is presented in Figure 3, where the upper soil layer has a thickness of 0.5 m. According to field experience, the soils in the lower layer are usually well consolidated, and they possess smaller pore sizes and lower permeability than the upper soil layer [21]. The two sets of material parameters adopted for these two soil layers are provided in

Table 1. Material parameters for the seepage–erosion modeling of the 1D layered soil column.

Parameter	Upper Layer	Lower Layer	Parameter	Upper Layer	Lower Layer
Initial saturated conductivity $K_0$ (m/s)	3.21 × 10−6	8.42 × 10−7	Initial fine fraction $x_{e’0}=m_{e’0}/{\rho }_S$	0.61	0.71
Initial porosity $n_0$	0.5	0.5	Erosion coefficient $k_e$	0.7	0.0
Maximum saturation $S_{sat}$	1	1	Soil erosion sensitivity ${\alpha }_{er}$	1.5	-
Residual saturation $S_{res}$	0.157	0.28	Reference velocity $v^{*}$ (m/s)	6 × 10−8	-
SWCC fitting parameter ${\alpha }_v$ (${\mathit{kPa}}^{-1}$)	0.105	0.086	Reference ultimate fine fraction $x_{e’\infty }^{*}$	0.58	-
SWCC fitting parameter $n_v$	1.547	1.611	Deposition coefficient $k_d$	0.0	0.4
Soil intrinsic density ${\rho }_S$ (${\mathit{kg}/m}^3$)	1370	1500	Permeability damage parameter $k_r$	0.05	0.05
Initial erodible-fines mass fraction $m_{e’0}$ (${\mathit{kg}/m}^3$)	831	1072	Permeability damage parameter $n_r$	3	3

. For validation reasons, the physical and hydraulic parameters are taken directly from Wu et al. (2016) [14] and Ma et al. (2010) [35], in which the initial erodible-fines mass fractions are calculated as the mass fractions of particles with diameters smaller than 0.075 mm according to the soil layers’ grain size distribution curves. The parameters related to fines migration are mainly adjusted from Lei et al. (2017) [28], in which the same erosion law is adopted. The reference velocity is calculated by $v^{*}=K_0{i}^{*}$, with $i^{*}$ being the critical hydraulic gradient. The critical hydraulic gradient is set to 0.1 in Lei et al. (2017); in this section, the critical hydraulic gradients for the upper and lower layers are taken as 0.02 and 0.07, respectively. To simulate the multiple rainfall-induced seepage–erosion effects over a long period in a short-time, single rainfall event, the soil erosion sensitivity parameters were increased to accelerate the seepage–erosion process.

In the example, the following three simulation scenarios were designed: (1) the reference (R) case, in which the internal erosion is neglected (all the erosion parameters ($k_e^{}$) and deposition parameters ($k_d^{}$) in Table 1 are set to zero ($k_e^{up}=k_e^{low}=k_d^{up}=k_d^{low}=0$)); (2) the erosion (E) case, in which the fine particles are eroded, but without deposition (all the deposition coefficients in Table 1 are set to zero ($k_e^{up}\ne 0$, $k_e^{low}=k_d^{up}=k_d^{low}=0$)); (3) the erosion–deposition (ED) case, in which both the fines erosion and deposition processes are taken into account ($k_e^{up}\ne 0$,$k_d^{low}\ne 0$, $k_e^{low}=k_d^{up}=0$). Initially, the soil column is unsaturated, and it is equilibrated with a uniform suction of 20 kPa. At time t = 0, a liquid flux of 8.3 × 10⁻⁷ m/s (3 mm/h, being the same value used in [14]) is applied on the top surface of the column, which triggers the rainfall infiltration, as well as the associated coupled seepage–erosion process within the unsaturated soil column.

The evolution profiles of the calculated field variables corresponding to different cases are presented in Figure 4. Along with the rainfall infiltration, the soils are saturated with water gradually, which increases the pore liquid pressure at the same time (Figure 4a,b). As the permeability of the upper layer is higher than that of the lower layer, water will transport in the upper layer at a relatively faster rate. Accordingly, water accumulates on the interface of the soil layers with time, which is reflected by the faster increase in the saturation above the layer interface (with the same porosity). For the reference (R) case, the numerical results of Wu et al. (2016) [14] are available for comparison. As can be seen in Figure 4a, our liquid pressure profile in the R case is very close to that of Wu et al. (2016) [14], which demonstrates the code’s capability in dealing with unsaturated seepage in layered soils.

In the E case, the fines erosion process within the upper layer is taken into account. As shown in Figure 4c, the mass content of the erodible fines within the upper layer decreases along with the rainfall infiltration, which increases the corresponding soil porosity simultaneously based on the mass conservation law of the solid phase (Figure 4d) (corresponding to the surface coarsening phenomenon reported experimentally in Cui et al. (2014) [7]). As the porosity increases, the permeability of these surface soils increases at the same time, according to Equation (8), which is reflected explicitly in the sharp increase in the permeability ratio near the column top (Figure 4f). The increase in the permeability near the surface of the upper layer leads to accelerated water infiltration into the deeper soil column, which is then partially trapped above the layer interface, as the lower layer is less permeable and is not affected by the fines erosion. Accordingly, the liquid pressure and saturation around the layer interface in the E case accumulate much faster than those in the R case, especially after t = 20 h, when the erosion is becoming serious.

In the ED case, the fines deposition process at the layer interface is further accounted for. Due to the distinct pore structure and size difference between the two layers of soils, the liquidized fines transported with water are prone to be captured around the layer interface, which is achieved by setting a non-zero deposition coefficient for the lower soil layer in the numerical model. In such a way, the fine particles eroded from the upper layer will be partially captured by the lower soil layer. With time, more and more fines will deposit around the layer interface (Figure 4e), which block the pores and reduce the soil permeability significantly according to Equation (8) (Figure 4f); accordingly, it becomes more and more difficult for water to infiltrate into the lower soil layer. As shown in Figure 4a,b, the liquid pressure and saturation on the interface at t = 40 h for the ED case are apparently higher than those for the E case, and much higher than those for the R case, which implies that it would become the most serious scenario for a layered deposited slope under rainfall.

4. 2D Simulation of Rainfall-Induced Seepage–Erosion Process within Layered Slopes

To further investigate the rainfall-induced seepage–erosion patterns, as well as their impact on the stability state of the layered deposited slopes, serials of two-dimensional finite element simulations of the coupled seepage–erosion processes within a layered slope under rainfall conditions were conducted and are described in this section.

In Figure 5, the finite element mesh of the simulated two-layer deposited-soil slope with a slope angle of 35° is presented. The slope section is 20 m long and was meshed with 25 segments in the horizontal direction (Lh = 0.8 m); the thickness of the upper layer was set as 1 m and was meshed with 8 segments in the vertical direction (Lv_up = 0.125 m); the thickness of the lower layer was chosen as 2 m, and it was meshed with 10 segments with the progression ratio 1.1 from the layer interface to the slope bottom in the vertical direction (Lv_low varied from 0.126 m to 0.296 m). A brief mesh sensitivity analysis is presented in Appendix A.

The slope is initially unsaturated with an initial pore-water pressure of −45 kPa at the model bottom and an initial pore-water pressure of −75 kPa at the slope surface. At time t = 0, an inward liquid flux of 1.8 × 10⁻⁶ m/s (6.5 mm/h) is imposed on the slope surface, which triggers rainfall infiltration, as well as the subsequent coupled seepage–erosion process.

Infinite-slope stability analysis is a commonly used tool for evaluating shallow landslides because of its simplicity and practicability [15]. This approach is valid in landside cases with small depths compared to lengths and widths, and it was adopted by Tsai and Chiang (2013) [15] and Wu et al. (2016) [14] for evaluating the rainfall stability of two-layered slopes. A classic infinite-slope model is illustrated in Figure 6. Lu and Godt (2008) [36]’s infinite-slope stability analysis model developed for unsaturated slopes is presented in Equation (9). Based on the Mohr–Coulomb failure criterion, an infinite slope is assumed to fail at an certain depth ($H_s(z)$) when the dimensionless safety factor ($F_{\mathrm{s}}$), representing the ratio of the resisting stress (due to friction, cohesion, and suction) to the acting stress (due to gravity), is smaller than 1:

$$F_{\mathrm{s}}(z,t)=\frac{\tan \phi }{\tan \beta }+\frac{2c}{\gamma H_s(z)\sin 2\beta }-\frac{S_L(z,t)p_L(z,t)}{\gamma H_s(z)}(\tan \beta +\cot \beta )\tan \phi$$

(9)

In the above equation, $\beta$ is the slope angle; $\phi$ is the soil internal-friction angle; $c$ is the cohesion of the soil; $\gamma$ is the average unit weight of the soil; $H_s(z)$ is the vertical depth of the current evaluation position ($x$, $z$) from the slope surface. The first term on the right-hand side of the above equation is due to the internal frictional resistance of the soil, the second term is due to the cohesion, and the third term is due to the suction stress ($-p_L$). Because suction stress is tensile in nature, the overall effect is to hold the soil together, thereby increasing the factor of safety; in contrast, the effect of positive pore-water pressure on the factor of safety is generally the opposite [36]. To evaluate the impact of the soil properties and rainfall conditions on the slope stability, the pore-water pressure ($p_L$) and saturation ($S_L$) fields obtained via numerical modeling are integrated into this infinite-slope model to account for the rainfall-induced suction changes; further, to account for the erosion-induced damage, the cohesion ($c$) and internal-friction angle ($\phi$) in the above equation are updated along with the newly calculated porosity ($n$) according to $c=c_0(1-n)/(1-n_0)$ and $\phi ={\phi }_0(1-n)/(1-n_0)$, with $c_0$ and ${\phi }_0$ being the initial cohesion and internal-friction angle corresponding to the initial porosity ($n_0$).

As described in Section 3, three types of scenarios are taken into account based on their internal-erosion characteristics. Namely, the R case serves as the reference case, in which both the erosion and deposition processes are neglected ($k_e^{up}=k_e^{low}=k_d^{up}=k_d^{low}=0$); the E case intends to investigate the pure erosion effect on the seepage process ($k_e^{up}\ne 0$, $k_e^{low}=k_d^{up}=k_d^{low}=0$); the ED case takes both the fines erosion and deposition processes into account ($k_e^{up}\ne 0$, $k_d^{low}\ne 0$, $k_e^{low}=k_d^{up}=0$). The soil material parameters necessary for these simulations are presented in

Table 2. Material parameters for the seepage–erosion modeling of the 2D layered slope.

Parameter	Upper Layer	Lower Layer	Parameter	Upper Layer	Lower Layer
Initial saturated conductivity $K_0$ (m/s)	3.66 × 10−5	1.42 × 10−6	Erosion coefficient $k_e$	0.022	0
Initial porosity $n_0$	0.35	0.35	Soil erosion sensitivity ${\alpha }_{er}$	4.76	-
Maximum saturation $S_{sat}$	1	1	Reference velocity $v^{*}$ (m/s)	3.6 × 10−7	-
Residual saturation $S_{res}$	0.17	0.16	Reference ultimate fine fraction $x_{e’\infty }^{*}$	0.193	-
SWCC fitting parameter ${\alpha }_v$ (${\mathit{kPa}}^{-1}$)	0.04	0.105	Deposition coefficient $k_d$	0	0.4
SWCC fitting parameter $n_v$	1.27	1.547	Permeability damage parameter $k_r$	0.05	0.05
Soil intrinsic density ${\rho }_S$ (${\mathit{kg}/m}^3$)	2240	2240	Permeability damage parameter $n_r$	3	3
Initial erodible-fines mass fraction $m_{e’0}$ (${\mathit{kg}/m}^3$)	448	448	Initial cohesion $c_0$ ($\mathit{kPa}$)	2.5	2.5
Initial fine fraction $x_{e’0}=m_{e’0}/{\rho }_S$	0.2	0.2	Initial internal-friction angle ${\phi }_0$ ($\circ$)	25	25

, in which the seepage–erosion-related parameters are referenced from Lei et al. (2017) [28], and the physical and hydro-mechanical parameters are taken from Wu et al. (2016) [14] and Wang et al. (2017) [37].

The evolution profiles of some key variables at section A in different simulation scenarios are presented in Figure 7; the contours of the liquid pressure, liquid flux, and porosity at the end of the simulation (t = 20 h) are illustrated in Figure 8. As shown in Figure 7a, the slope was initially equilibrated with a pore-water pressure of −45 kPa at the slope bottom and a pore-water pressure of −75 kPa at the slope surface at t = 0. When the simulation started, rainfall infiltrated into the slope via the slope surface with the imposed liquid flux and moved downwardly under gravity according to Equation (7), which increased the pore-water pressure gradually with time (Figure 7a). The liquid saturation is linked to the pore-water pressure bijectively via the assumed soil water characteristic curve (SWCC) (Equation (6)). As the soils in the two layers have distinct water retention properties (SWCC parameters), there are significant discontinuities in the saturation profiles at the soil layer interface shown in Figure 7b.

As the permeability of the lower soil layer was set to be much smaller than the one of the upper layer (which is usually the case for deposited slopes), perched water started to accumulate above the layer interface at around t = 15 h (indicated by the positive water pressure and full saturation above the interface in Figure 7a,b). As shown in Figure 8a, the perched-water region featured with positive liquid pressure can be visibly identified on the 2D inclined slope, in which the perched-water depth gradually increased along the direction of the downward slope due to the accumulation of the upstream flux. Moreover, obvious subsurface stormflow [26] was generated right above the layer interface (Figure 8b), with the flow velocity becoming larger along the direction of the downward slope due to gravity.

As the soils in the upper layer are prone to internal erosion, this continuous subsurface stormflow can lead to significant erosion around the interface in the E case and ED case. As shown in Figure 7c, in the early stages (t = 5 h, 10 h), the erosion was mainly caused by the rainfall infiltration, which was concentrated near the slope surface with a relatively large liquid flux and long seepage time, leading to the so-called “surface coarsening phenomenon” [7] near the slope surface (as already illustrated in Figure 4c of Section 3). After t = 15 h, the soils above the layer interface were nearly saturated, and a continuous stream flow was generated along the interface. As the liquid flux of this subsurface stormflow was much larger than the rainfall infiltration flux (see Figure 8b), the soils in these regions were eroded to a more serious degree than the soils near the slope surface with time (Figure 7c). What is more, the soils around the downstream region were eroded more seriously than those around the upstream region in the 2D inclined slope (shown by the larger porosities in Figure 8c), as they were exposed to larger and longer stormflow fluxes (Figure 8b).

Due to the pore size and grading difference, the fines eroded from the upper layer will inevitably be blocked by the lower layer. In the ED case, the fines deposition process was taken into account. As shown in Figure 7d, the fines deposited around the soil interface increased with time gradually, and this became significant after t = 15 h. The fines deposited around the interface blocked the pores and reduced the soil permeability of the lower layer significantly (as already demonstrated in Figure 4f), which restricted the further infiltration of water into the lower soil layer. Accordingly, the perched water in the ED case accumulated faster than that in the E case and R case, as shown in Figure 7b and Figure 8a. Along with the gradual permeability damage of the lower layer, more and more infiltrated water flowed downwardly along the layer interface (Figure 7b), which, in turn, led to more serious soil erosion around the layer interface than in the E case (Figure 7c and Figure 8c). As a direct result of the erosion, the porosity of the soils increased with the erosion simultaneously (Figure 7e and Figure 8c), which weakened the soils along the layer interface gradually.

In Figure 7f, the evolution profiles of the safety factors along section A in different scenarios are illustrated. As can be seen, prior to the generation of the continuous subsurface stormflow (t < 15 h), the safety factor profiles in all the cases were almost identical to each other. Later, the pore-liquid pressure around the layer interface increased rapidly due to the accumulation of perched water; meanwhile, the generated stream flows weakened the soils along the layer interface continuously. Both these processes are negative for slope safety according to Equation (9). Accordingly, the ED case with both the fines erosion and deposition processes taken into account would be the most dangerous scenario for layered deposited-soil slopes. In the E case, the soil erosion was weaker, and the perched water was less accumulated around the layer interface; therefore, the E case with the pure erosion process accounted for is the second most dangerous scenario. As shown in Figure 7f, the safety factor profile in the ED case crossed the critical safety line ($F_{\mathrm{s}}=1$) at the region around the layer interface (instable region) at t = 20 h, while the safety factor profile in the ED case just reached the safety line. As for the R case, though its safety factor profile did not reach the safety line during the simulation time, the slope will probably fail after a long enough rainfall infiltration with the low-permeability lower layer adopted in this example.

It should be noted that, as the internal-erosion-induced weakening of soils can be accumulated during the seepage–erosion process, layered deposited-soil slopes prone to internal erosion may fail gradually after multiple short-duration rainfall events. Further, the fines deposition process may turn a relatively highly permeable lower layer into a less permeable soil layer, which could also lead to the gradual failure within the layered slope with an initial highly permeable lower layer after a long duration. As the internal fines erosion and deposition processes within these layered slopes are difficult to identify and evaluate, the failure of these layered deposited-soil slopes is difficult to predict, and therefore special attention has to be paid to their stability analysis.

5. Conclusions

In this paper, the main physical processes within the classic two-layer deposited-soil slopes under rainfall infiltration are summarized, and a coupled seepage–erosion finite element model is established to analyze the interactions between the rainfall infiltration process and the seepage-induced internal-erosion process within layered deposited-soil slopes. Serials of coupled seepage–erosion simulations, including a scenario without internal erosion, a scenario with pure fines erosion, and a scenario with both fines erosion and deposition, were undertaken; the evolution profiles of the key field variables within these scenarios, as well as the slope safety factors, were compared and analyzed. The research indicates the following:

(a)

The seepage-induced internal-erosion process within layered deposited slopes will impose a significant influence on their failure patterns and mechanisms under rainfall. Except for the surface coarsening effect with accelerated infiltration reported previously (e.g., in Cui et al. (2014) [7] and Lei (2022) [30]) for single-layer deposited slopes, the fines migration around the soil layer interface can lead to the generation of a continuous subsurface stormflow above the layer interface, which would increase the pore-water pressure and weaken the soil strength and, hence, decrease the slope stability along the layer interface continuously;

(b)

The layered deposited-soil slope prone to internal erosion may fail gradually after sufficiently long rainfall infiltration due to the accumulated weakening of the soils. Further, the fines deposition process may turn a relatively highly permeable lower layer into a less permeable soil layer, which could also lead to the gradual failure within the layered slope with an initial highly permeable lower layer;

(c)

Due to the concealed and impactful nature of seepage-induced internal erosion, special attention should be paid when analyzing the stability of these layered deposited-soil slopes. The coupled seepage–erosion FEM presented in this paper could serve as a useful mathematical tool to analyze the intrinsic failure mechanism and stability evolution of such slopes under rainfall infiltration.