Effect of Rainfall, Runoff and Infiltration Processes on the Stability of Footslopes

Abstract

To analyze the effect of runoff on shallow landslides, a model coupling one-dimensional rainfall–runoff and two-dimensional infiltration was established to simulate rainfall, infiltration, and runoff processes. Based on Bishop’s limit equilibrium method, the slope failure of a hypothetical footslope was studied. First, conditions with and without inflow were compared. The results reveal a remarkable difference in factors of safety (FS) between the two conditions, suggesting that considering the effect of runoff is crucial for landslide modeling. In terms of a series of tests of the various magnitudes, durations, lag-time, and peak position of the hydrograph, analyses show that larger inflow leads to more accumulated infiltration and triggers landslides earlier. A long-term duration inflow decreases the stability more than short intensive inflow does. With subsequent surface inflow, slope failure may occur after rainfalls stop, owing to the inflow, and the shape of inflow hydrographs could slightly affect the variance in FS. Results also indicate the necessity of considering the surface runoff when using a numerical model to analyze landslide, particularly on a footslope.

1. Introduction

Landslides occur in numerous locations worldwide and present serious threats to lives and the economy. There is a general consensus that rainfall is the most common external factor that stimulates shallow landslides. Therefore, the mechanism underlying rainfall-induced landslides has been widely studied and discussed in recent decades [1,2].

A relationship between water and landslides, known as “landslide hydrology,” has been proposed [2]. In the hydrological cycle, water can originate from many sources, such as precipitation, snowmelt, surface water and groundwater. Water infiltrating soil surfaces leads to low matric suction and a reduction in the shear strength of soil, thereby increasing the instability. The groundwater seepage has been broadly and deeply analyzed using numerical models, wherein different soil conditions and facilities have been considered [3,4,5,6,7]. However, it is difficult to completely reproduce the mechanism underlying rainfall-triggered shallow landslides using a numerical model. With the mechanism greatly simplified in the earliest studies, landslide models have gradually progressed in recent decades and currently include more mechanisms. For example, Iverson [8], Baum et al. [9], and Tsai and Yang [10] considered pore-pressure diffusion as a near-saturation response. Tarantino and Bosco [11], Collins and Znidarcic [12], and Tsai et al. [13] used the Richards equation and the extended Mohr–Coulomb failure criterion [14] to develop shallow landslide models for both saturated and unsaturated soils. Based on previous research, there has been progress in the clarification of landslide behavior, and the improved modeling capabilities allow researchers to explore details.

The influence of rainfall infiltration on landslide has been studied using numerical models. Tsai [15] used a modified version of Iverson’s one-dimensional (1D) model to assess the influence of different rainfall patterns on the groundwater table and slope stability in near-saturated soils. Tsai and Wang [16] demonstrated the effect of different patterns of rainfall on shallow landslides in unsaturated soils. Chen et al. [17] further investigated the effects of the rainfall duration, rainfall amount, and lateral flow-induced slope failures using a vertical two-dimensional (2D) numerical landslide model. The relationship between rainfall infiltration and slope stability has been discussed extensively.

Recently, hydrological and landslide susceptibility models were combined [18], and the results revealed that the runoff depth may be an appropriate analysis factor, instead of the rainfall depth or maximum rainfall intensity. Chiu et al. [19] combined a landslide model with a surface flow simulation of a small catchment and observed that the accumulated surface water increases the possibility of landslide, especially in the downstream slope. The latter study emphasized the importance of considering the effect of runoff on the slope stability during the analysis. However, the simulation was based on a regional watershed, and it was difficult to control the runoff conditions at an affected slope. Therefore, different surface flow conditions have not been studied yet.

The objective of this study is to analyze the effect of runoff on the slope failure under different conditions including magnitudes, durations, lag time, and peak positions of hydrographs. This study aimed to analyze the variety in footslope stability with respect to time. A footslope refers to a slope near the toe of a hill (Figure 1). The geographical position allows the converged water from the upper areas of the catchment to flow through the footslope; therefore, compared to other slopes, a footslope is more likely to be affected by runoff. A physical-based model was developed to simulate the rainfall–runoff and slope instability based on the limit equilibrium method. The model is governed by the coupling of the 1D kinematic wave equation and 2D Richards equation. A case study demonstrating the different results of conditions with and without surface inflow on the slope stability of shallow landslides is presented, and various surface inflow hydrographs have been adopted to test the slope stability in different conditions.

2. Method and Materials

2.1. Mathematical Basis

The proposed model can be divided into two parts: (1) hydrological module and (2) soil failure module. The former is used to simulate the rainfall, infiltration, and runoff, whereas the latter is used to calculate the slope stability. These two modules are physical-process-based to determine the physical mechanisms triggering landslides. The theory underlying 1D rainfall–runoff, 2D infiltration equation, and slope stability analysis is briefly described in this section.

2.1.1. Hydrological Module

The hydrological module mainly simulates the rainfall infiltration and overland flow processes. Surface flow is obtained from kinematic wave approximation, ignoring inertia and pressure terms. The infiltration simulation considers both saturated and unsaturated conditions to compute the pore-water pressures by solving the 2D Richards equation. The kinematic wave and infiltration equations are briefly described below.

Kinematic Wave Equation

The actual physical flow processes may be quite complex, but a simplified hydrodynamic model is sufficient for practical conditions (Henderson [22], Gunaratnam and Perkins [23], and Miller [24]). Kinematic wave equations [25] are often used to simulate the rainfall–runoff process in small- and average-sized basins with steep slopes. The details and the derivation of these equations are provided by Chow et al. [26] and Henderson [22].

The equation governing flow over a plane under the assumption of a kinematic wave is

$$\frac{\partial q}{\partial X}+\frac{\partial h_s}{\partial t}=i-f,$$

(1)

where $h_s$ denotes the depth of the slope surface flow; $i$ is the rainfall intensity; $f$ is the infiltration rate; $t$ is the time; $X$ is the distance downslope; and $q$ is the flow rate per unit width and can be obtained as follows

$$q=\frac{\sqrt{S_0}}{n}\times h_s^{5/3},$$

(2)

where n is Manning’s roughness coefficient, and $S_0$ is the bed slope.

Infiltration Equation

Infiltration flow is generally considered to be an essential triggering mechanism that must be addressed when studying the slope failure following a period of heavy rainfall. In this study, both saturated and unsaturated infiltration flows were considered when evaluating the pore-water pressure by solving the 2D Richards equation. The equation that can be used to calculate the groundwater flow in response to rainfall infiltration on a hillslope with a local rectangular Cartesian coordinate system is [27]

$$\frac{\partial \psi }{\partial t}\frac{\partial \theta }{\partial \psi }=\frac{\partial }{\partial x}\left[K_x\left(\psi \right)\left(\frac{\partial \psi }{\partial x}\right)\right]+\frac{\partial }{\partial z}\left[K_z\left(\psi \right)\left(\frac{\partial \psi }{\partial z}+1\right)\right],$$

(3)

where $\theta$ is the moisture content;$\psi$ is the pore-water pressure head; and $K_x$ and $K_z$ are the hydraulic conductivities in the $x$ and $z$ directions, respectively, a function of the soil properties presented in Equation (8). Equation (3) can be solved by using the appropriate initial and boundary conditions.

Regarding an initial steady state with a water table of $d_Z$ in the vertical direction, the initial condition in terms of the pressure head can be expressed as

$$\psi =Z-d_Z,$$

(4)

where $Z$ is the elevation. In addition, the surface of the hillslope subjected to the rainfall yields

$$f=-K_x\left(\psi \right)\left(\frac{\partial \psi }{\partial x}\right)\sin \alpha -K_z\left(\psi \right)\left(\frac{\partial \psi }{\partial z}-1\right)\cos \alpha ,$$

(5)

where $\alpha$ is the slope angle. When the ground surface is saturated, the pressure head $\psi$ of the ground surface is equal to the depth of flow; thus, the boundary condition can be expressed as

$$\psi =h_s$$

(6)

Solving Equations (3)–(6) requires a relationship among pressure, moisture content, and hydraulic conductivity. The water retention curve proposed by van Genuchten [28] was used in this study for the hydraulic conductivity $K=K_x=K_z$ for isotropic soils

$$S=\frac{\theta -{\theta }_r}{{\theta }_s-{\theta }_r}={\left[\frac{1}{1+{\left(\xi \psi \right)}^N}\right]}^M$$

(7)

$$\frac{K\left(\theta \right)}{K_s}={\left(\frac{\theta -{\theta }_r}{{\theta }_s-{\theta }_r}\right)}^{1/2}{\left\{1-{\left[1-{\left(\frac{\theta -{\theta }_r}{{\theta }_s-{\theta }_r}\right)}^{1/M}\right]}^M\right\}}^2,$$

(8)

where $S$ denotes the degree of saturation; $K_s$ is the saturated hydraulic conductivity; ${\theta }_s$ is the saturated moisture content; ${\theta }_r$ is the residual moisture content; and $\xi$, $N$, and $M$ are the fitting parameters ($M=1-1/N$). The calculated pore-water pressures were used to evaluate the factor of safety (FS) of a hillslope.

2.1.2. Soil Failure Module (Slope Stability Analysis)

The slope stability is commonly analyzed using the limit equilibrium method of slices [29] and based on the consideration of the moment equilibrium of the soil mass. The limit equilibrium methods assume that the potential failure surface is governed by Mohr–Coulomb relationships between shear strength and the normal stress on the failure surface. The soil mass can be subdivided into a number of vertical slices that exhibit a width b above an assumed circular slip surface (Figure 2). The maximum shear stress $\tau$ of the unsaturated soil acting at the lower boundary of a slice is related to the total normal stress $\sigma$, according to the Mohr–Coulomb failure criterion [30], and can be expressed as follows

$$\tau =c^{\prime }+\left(\sigma -u_{\alpha }\right)\tan {\varphi }^{\prime }+\left(u_{\alpha }-u_w\right)\tan {\varphi }^b,$$

(9)

where $c^{\prime }$ denotes the effective cohesion; ${\varphi }^{\prime }$ and ${\varphi }^b$ are the shearing resistance and friction angles with respect to the matric suction, respectively, and $u_{\alpha }$ and $u_w$ are the pore-air and pore-water pressures, respectively. In Bishop’s method [29], the force transmitted between adjacent slices is assumed to be strictly horizontal. The equilibrium of moments relative to the center of the slip circle can be expressed by equating the sum of the moments of the weight of each slice relative to the center of the circle to the sum of the moments of the shearing forces at the bottom of the slices. The depth-averaged unit weight $\gamma$ in a slice can be expressed as

$$\gamma =\frac{1}{h}{{\displaystyle \int }}_0^h\left[\left(1-\theta \right)·{\gamma }_w·G_s+\theta ·{\gamma }_w\right]dh,$$

(10)

where $G_s$ denotes the specific gravity of a solid; $h$ is the height of each slice; and ${\gamma }_w$ is the unit weight of water.

When using Equation (10) together with the shear strength of unsaturated soil given by Equation (9) and assuming that the pore-air pressure is atmospheric, the Factor of safety ($\mathrm{FS})$ for the equilibrium of moments based on Bishop’s method [29], which was derived by Chen [17], can be expressed as

$$\mathrm{FS}=\frac{{\sum }^{}\left[\frac{\left(c^{\prime }+\gamma ·htan{\varphi }^{\prime }-{\gamma }_w·{\psi }_{c}tan{\varphi }^{\prime }-{\gamma }_w{\psi }_p\tan {\varphi }^b\right)}{cos\alpha \left(1+tan\alpha tan{\varphi }^{\prime }/FS\right)}\right]}{{\sum }^{}\gamma ·hsin\alpha },$$

(11)

where the soil is unsaturated when the groundwater pressure head becomes negative and ${\psi }_c$ equals $\psi$, which can be obtained using Equation (3), whereas ${\psi }_p$ is zero. Conversely, when the groundwater pressure head becomes positive, ${\psi }_p$ is identical to $\psi$, ${\psi }_c$ is zero, and the soil is saturated. When FS is less than one, it means the slope is unstable. Equation (11) indicates that slope failure occurs not only in saturated soil due to an increase in the positive groundwater pressure head, but also in unsaturated soil due to a decrease in the negative groundwater pressure head, indicating the dissipation of the matric suction. As shown in Figure 2, Bishop slip circle is to be determined by a center and radius which has the overall minimum factor of safety calculated by Equation (11).

2.1.3. Computational Procedure

The numerical model is solved using the finite difference method; the discretization of the kinematic wave equation is computed using the Preissmann four-point finite difference scheme [30], whereas the Richards equation is computed with a fully implicit scheme [19,31]. An iterative procedure (Figure 3) is adopted to solve Equations (1) to (11). The iterative approach alternates the solutions of the infiltration and overland flow equations at each time step until convergence is reached. The groundwater pressure head of a hillslope is obtained by assuming that the infiltration rate equals the rainfall intensity, as presented in Equation (5), and the surface runoff is solved in succession using Equations (1) and (2). If the pressure head on the slope surface is less than or equal to zero, surface runoff does not occur, and the calculation progress to the next time step. However, if the calculated pressure head or water depth on the slope surface is larger than zero, that is, ${\psi }_0>0$ or $h_s>0$, and the error between the pressure head and the water depth is assessed. If the error is greater than the minimum permissible error, $\left({\psi }_0+h_s\right)/2$ transfer into a boundary condition of infiltration for the recalculation within the same time step. Finally, the calculated pore-water pressures at the slopes are substituted into Equation (11), and FS is obtained. The physical-based method for the landslide simulation has a high accuracy under field and laboratory conditions. However, it ignores other factors, such as the vegetation.

2.2. Description of the Study Case

A hypothetical slope is used for demonstration. Figure 4 presents the geometry of the slope with an angle of 31° and height of 6 m. The red dash line in Figure 4 is a schematic line of potential circular slip that has the overall minimum factor of safety. The slope is set as a deep slope, where the gradient is larger than the friction angle. The boundary conditions are as follows: $\overline{\mathrm{ab}}$ is an impervious bottom (rock mass); $\overline{\mathrm{ae}}$ and $\overline{\mathrm{bd}}$ are the Dirichlet boundaries with total water heads; boundary $\overline{\mathrm{df}}$ has no flux; and the ground surface $\overline{\mathrm{ef}}$ is subjected to rainfall. The initial groundwater table $\overline{\mathrm{cd}}$ was 4 m below the ground surface of the hillslope. The soil is assumed to be sandy loam [28], and the parameters are as follows: ${\theta }_s=0.47$; ${\theta }_r=0.17$; $K_s=0.031\mathrm{m}/\mathrm{h}$; $N=2$; $M=0.5$; $\xi =0.01$; $c^{\prime }=500\mathrm{N}/{\mathrm{m}}^2$; ${\varphi }^{\prime }=26°$; ${\varphi }^b=13°$; and $G_s=2.65$. All symbols refer have been defined in Appendix A. This hypothetical slope is adopted to all the tests presented in the study with various surface inflow conditions.

3. Results and Discussions

Firstly, two cases were selected for comparison: (1) a uniform rainfall event with an intensity of 25 mm/h and a surface inflow of 0.0001 m²/s over 12 h and (2) the same rainfall conditions but without surface inflow. The results of FS, accumulated infiltration and pore pressures are presented in Figure 5 and Figure 6. In Figure 5a, when the inflow is considered, the FS decreases faster and a landslide is triggered at the 9th hour; however, in the condition wherein the inflow is not considered, the FS is always larger than 1.0. The lowest FS of the two cases are 1.03 and 0.73, respectively. The variation in accumulated infiltration is presented in Figure 5b, indicating that, in the case considering inflow, the infiltration is higher and the difference of the two cases increases with time. In the case without considering the inflow, the final accumulated infiltration is 34,985 cm², however, in the case that considers the inflow, it is 48,769 cm², i.e., 1.4 times the other case. Figure 6 presents the pore pressures of the two cases. The pore pressures in Figure 6c,d increase faster than those in Figure 6a,b, that is, they are more likely to cause landslide. The abovementioned results indicate that, firstly, the hillslope inflow may increase the slope instability, and the phenomenon is more likely to occur on footslopes than on upper slopes due to the possibility of runoff. In addition, from the numerical modeling perspective, the mechanism considering the runoff is crucial for these footslopes.

According to the previous results, the effect of the inflow on the slope stability is notable. In situ, the magnitude, duration, lag time, and peak of the hydrograph [32] may vary according to the storm conditions and catchment characteristics. Therefore, different conditions of inflow and associated effects on the slope stability are discussed in the following sections.

3.1. Magnitude of Surface Inflow

The size of the upper catchment area and the magnitude of rainfall cause different magnitudes of surface inflow. Four different magnitudes of surface inflow (0.00001 m²/s, 0.00002 m²/s, 0.00004 m²/s, and 0.0001 m²/s) were selected in this study to analyze the change in the FS with time. These values are in a reasonable range as a result of estimating the length of upper slope are about tens to hundreds of meters. In addition, the slope is subjected to a uniform rainfall event (200 mm, 12 h). As presented in Figure 7a, a landslide is triggered in all cases, except for the one with the smallest inflow. The times at which landslides are triggered is 11.7, 9.5, and 9 h for the cases with 0.00002, 0.00004, and 0.0001 m²/s inflow, respectively. In addition, although these cases all exhibit that the lowest FS appears at the same time (12th h), the time at which the landslides are triggered differ, depending on the magnitude of the inflow. Figure 7b presents the variance in accumulated infiltration with time, and reveal that larger inflow leads to more accumulated infiltration. In the tests, the accumulated infiltrations are 30,638, 36,304, 45,383, and 48,390 cm² to the inflow, from small to large. The pore pressure also reveals the effects of different magnitudes of inflow, as presented in Figure 8, indicating that the magnitude of inflow accelerates the wetting front propagating downward. Therefore, a larger inflow increases the probability of triggering a slope failure and may bring forward the trigger time.

3.2. Duration

Three scenarios with different durations but the same inflow volume were considered: (i) Scenario 1: 0.00016 m²/s for 3 h, (ii) Scenario 2: 0.00008 m²/s for 6 h, and (iii) Scenario 3: 0.00004 m²/s for 12 h. The designs were used to observe whether an intense, short inflow affect the FS more, or a mild, longterm inflow does. We firstly analyzed the variance in FS with time, as presented in Figure 9a. In the first 3 h, the differences in FS among three scenarios are small; the values of FS are 1.13, 1.14, 1.15, respectively, for Scenario 1, 2, and 3, and Scenario 1 has the lowest FS. From 3 to 6 h, Scenario 2 yields the lowest FS. Similarly, because the surface flow increases the infiltration rate, Scenario 1 has a higher accumulated infiltration in the first 3 h (Figure 9b). However, after the runoff stops, for example, at the 4th h, the infiltration of Scenario 1 is still higher than that of Scenario 3, but the FS rebounds from 1.13 to 1.192 in a short time and exhibits a different trend than the infiltration, due to surface pressure dissipation. Figure 10 presents the results for the pore pressure that explains the rebound of FS. Between the 3rd h (Figure 10a) and the 4th h (Figure 10b), the surface pressure of Scenario 1 changes from 0 to −30 cm, whereas those of Scenarios 3 and 2 remain at 0 (Figure 10c,d). The results emphasize the effect of continuous inflow on the stability, and long-term duration is more likely to trigger a landslide.

3.3. Delayed Surface Inflow

In this section, the start and end times of the inflow of all scenarios differ by 2 h (Figure 11). This arrangement was used to analyze whether or not landslides may occur after the rain stops, attributed to the subsequent runoff. The slope is subjected to a uniform rainfall event (200 mm, 12 h), the same as the previous sections.

In scenarios with delayed surface inflow (delay: 0, 2, and 4 h), landslides are triggered during the rainfall event (within 12 h), as shown in Figure 12. However, in the other three scenarios, landslides are not triggered during the rainfall event but rather when the slope receives the subsequent surface flow. The pore-pressure variation in Figure 13 indicates that the pressure under the surface continues to increase after the rainfall stops, which is due to surface water, and the wetting front propagates downward. The results indicate that landslides can occur after the rain stops and delayed runoff is more likely to cause landslide.

3.4. Shape of Inflow

Different basin shapes may lead to different types of hydrograph [32]. In this study, the inflow hydrograph is simplified into three types based on the peak position of the hydrograph, as shown in Figure 14: (1) peak appears at the beginning, (2) peak appears at middle, and (3) peak appears at the end. Figure 15a shows that the differences among the three types gradually expand until about the 2.6 h, the FS of Type 1, Type 2, and Type 3 are 1.20, 1.23 and 1.26, respectively. Moreover, the times at which landslides occurred are 9.3, 9.7, and 10.0 h. The FS of the Type 1 decreases faster than that of the others and triggers the landslide slightly earlier than the others. However, if the peak value of the inflow is a lower value, such as 0.00008 m²/s, Figure 15b indicates the possibility that only Type 3 triggers a landslide. However, compared with the tests in the previous sections, only slight differences can be found in the FS trends among the three types.

4. Conclusions

This study focused on the effect of surface water on the slope stability. In this study, a 2D landslide model coupling the 1D rainfall–runoff and 2D infiltration equations was developed. Based on the physical-based model, the difference between considering inflow and not considering inflow on the FS was elucidated. With inflow consideration, the FS decreased much faster than with no inflow consideration. This revealed that the consideration of the runoff is crucial for modeling, particularly when footslopes are analyzed. In addition, based on a series of simulations of the hypothetical slope, the magnitude, duration, and lag time affect the FS and time of landslide occurrence, although the influence of the peak position of hydrograph isminor. The magnitude tests showed that a larger inflow has a lower FS and triggers landslide earlier. Long-term duration inflow decreases the stability more than short intensive inflow does, in the duration tests. In addition, the delay surface inflow tests indicate that landslides may occur after the rainfall stops, attributing to the lingering inflow.

In practice, the result of the study supports the usages of drainages and agrees that hillslope management should contrive to reduce runoff, because uncontrolled runoff can cause slope failure. Besides, when using numerical model to analyze the shallow landslide on footslopes, we strongly suggest modelling the surface flow and considering the inflow from upper slopes as well, and the effect of runoff on the boundary condition should be calculated in the model, owing to the remarkable influence we exhibited in the study.