Investigations of Unsaturated Slopes Subjected to Rainfall Infiltration Using Numerical Approaches—A Parametric Study and Comparative Review

Abstract

In the past 30 years, research on rainfall-induced landslides has grown remarkably. The contribution of matric suction to soil strength and the physics of water flow in unsaturated soils are widely accepted phenomena among researchers. However, the adoption of unsaturated soil mechanics in geotechnical engineering practice has been relatively slow, in part due to the practicality of design solutions available to the engineer. This paper conducts a parametric study on unsaturated silty slopes under a vertical steady flow rate to identify the suitable slope and hydrologic conditions to incorporate unsaturated conditions for preliminary stability analysis. Notably, the contribution of suction is most significant for silt/clay slopes with a water table located below the mid-height of the slope. For slopes with slope height ≥20 m and a fairly high water table, the slope height is a primary controlling factor of slope stability. Two case studies based on distinct failure mechanisms are presented to review the application of common geotechnical software in rainfall seepage and stability analyses of unsaturated slopes. Focus is placed on the pre-failure and failure stages of each case study. The slip surface search method, failure mode, and coupling approach integrated into each computer program caused notable differences in output results.

1. Introduction

Recent decades have witnessed an increased frequency of rainfall-induced landslides, particularly in tropical and sub-tropical climatic regions [1]. Along with other factors such as geomorphological conditions, seismic activity, and land use policies, the effects of climate change are expected to further exacerbate the problem with a projected rise in high-intensity short-duration rainfall events [2]. Rainfall-induced slope failures are mostly shallow in nature and mainly attributed to the progression of a wetting front and reduction of shear strength due to a decrease in matric suction in the unsaturated soil layer, which warrants the incorporation of matric suction in the stability analysis [1].

Despite advances in unsaturated soil mechanics, several hurdles limit its application in standard geotechnical engineering practice. Firstly, obtaining unsaturated soil parameters requires advanced laboratory tests which can be time-consuming and expensive. Although this has been improved by the development of pedotransfer functions (PTFs) that can indirectly quantify unsaturated soil properties, caution should be taken due to concerns regarding the reliability of collected data and estimation errors [3]. Secondly, the intricate nature of infiltration in unsaturated soils requires a solution of partial differential equations that describe the physical response of unsaturated soil to rainfall seepage. Many analytical infiltration models may involve complicated formulations that include infinite series or integrals or may only be applicable to simplified systems [4]. In this case, engineers should strike a balance between the practicality and accuracy of the infiltration model. Thirdly, more sophisticated modeling of rainfall infiltration warrants the use of advanced numerical models to capture the highly variable nature of soil properties and boundary flux conditions. This calls for a certain level of expertise and several input parameters for rainfall seepage analyses which have limited the widespread use of numerical simulation in common slope stability practices [5].

At the design stage, geotechnical engineers focus on the quantification of soil and slope geometric variables through a series of analyses from which the engineers must select a design solution. This normally takes the form of a parametric study [6]. For the unsaturated slope problem, various studies have discussed the relationship between suction forces and slope stability [7,8,9,10]; however, not much attention has been shown to quantifying the contribution of suction to slope stability over a wide range of slope and hydrological parameters [11]. This information could be beneficial to engineers during the design process.

The present paper is divided into two parts: Part I and Part II. Part I conducts a comprehensive parametric study of unsaturated slopes under rainfall infiltration using a practical and widely used one-dimensional (1-D) steady flow rate formulation by Lu and Griffiths (2004) [12]. The calculated pore water pressure profiles are exported to two-dimensional (2D) finite element limit analysis (FELA) software; Optum G2 [13] for slope stability analyses to obtain the factor of safety $FoS$. Several guidelines on the suitable application of unsaturated conditions in slope stability analysis are proposed which can be used by practicing engineers in an initial study.

Part II reviews the use of widely used and accepted commercial software programs RS2 [14], Slide 2D [15], SEEP/W, SLOPE/W [16,17], PLAXIS (LE) [18], PLAXIS 2D [19], and Optum G2 [13], on seepage and stability analysis of partially saturated slopes subjected to rainfall infiltration. Two case studies which range in complexity from an assumed slope case by Rahardjo et al. (2007) [9] to a full-scale landslide flume test conducted by Lee et al. (2019) [20] are presented and the output results are validated on the basis of recorded results and observations. To date, quite a number of studies have used computer software programs to simulate slope-related problems mainly due to rainfall seepage [10,21,22,23,24,25,26] yet, only a handful have presented a comparative view of geotechnical software programs to general slope problems [27,28], let alone rainfall-induced landslide cases. It is envisaged that a discussion on the results and limitations of each software program can serve as a design aid for practicing engineers and improve judgment of slope stability problems.

This paper focuses on the pre-failure and failure stages of unsaturated slopes which are relevant to the engineering design phase. Recently developed numerical methods such as the material point method (MPM), smooth particle hydrodynamics (SPH), and other meshless methods have been shown to effectively simulate multi-phase interactions involving large deformations such as the post-failure stages of landslides [29,30,31]. However, these novel methods are yet to gain large-scale use in the industry [32,33,34].

2. Theoretical Background

2.1. Infiltration Analysis in Soil Slopes

To simulate rainfall infiltration, traditional slope methods ignored suction forces in the unsaturated layer and assumed that the water table rises to the slope surface as a worst-case scenario. Although this is a simplified approach, it is unrealistic and overly conservative [1]. Over the past century, special efforts have been dedicated to developing infiltration models that can capture water flow in unsaturated and saturated flow systems by combining the principles of Darcian flow law [35], Equation (1), and the conservation of fluid mass [36].

$$Q=-k_w\frac{\partial h}{\partial y}$$

(1)

where Q = water flow rate, $k_w$ = coefficient of permeability with respect to the water phase, $\frac{\partial h}{\partial y}$ = hydraulic head gradient in the vertical direction.

2.1.1. Infiltration Conceptual Models

Green and Ampt (1911) [37] proposed an analytical solution for 1-D water flow into an initially dry uniform column of soil. The proposed model was based on a sharp wetting front hypothesis that assumed that the soil suction beyond the wetting front is constant with space and time and the water content and corresponding hydraulic conductivity before the wetting front is also constant with space and time. However, the sharp wetting front observation is not always valid, especially for fine-grained soils or soil with initial moist conditions [38]. The Green and Ampt model, initially developed for level ground with constant-head or ponded infiltration boundary conditions, was expanded to constant [39] and unsteady rainfall flux conditions [40], and also to sloping ground conditions [41].

Notably, Lumb (1962) [42] applied the wetting front concept to slope failures in Hong Kong. Under prolonged and heavy rainfall, the depth of the wetting front, $z_w$ is provided in Equation (2). Despite its use to locate the wetting front at the end of a wet season [10], Lumb’s equation does not account for slope angle inclination, rainfall intensity, or the nonlinear relationship between soil permeability coefficient and suction.

$$z_w=\frac{k_{s}t}{n\left(S_f-S_0\right)}$$

(2)

where $k_s$ = saturated hydraulic conductivity (m/s),$S_f$ = final degree of saturation, $S_0$ = initial degree of saturation, n = soil porosity, $t$ = time (s).

2.1.2. Analytical Solution of Richard’s Equation

For most slope infiltration problems, the governing equation for two-dimensional transient flow in unsaturated soil shown in Equation (3) is used [43]. This formulation known as Richard’s equation can also be extended to one-dimensional or three-dimensional flow problems.

$$\frac{\partial }{\partial x}\left(k_x\frac{\partial h}{\partial x}\right)+\frac{\partial }{\partial y}\left(k_y\frac{\partial h}{\partial y}\right)+q=\frac{\partial {\theta }_w}{\partial t}=m_w{\rho }_{w}g\frac{\partial h}{\partial t}$$

(3)

where $h$ = total hydraulic head (m), $k_{x,y}$ = unsaturated hydraulic conductivity in the x- and y-directions respectively (m/s), $q$ = applied unit flux (m/s), ${\theta }_w$ = volumetric water content, $m_w$ = coefficient of volume change (slope of SWCC), ${\rho }_w$ = density of water ($\raisebox{1ex}{$\mathrm{kg}$}\!\left/ \!\raisebox{-1ex}{${\mathrm{m}}^3$}\right.$), $g$ = acceleration due to gravity (m/s) and $t$ = time (s).

Equation (3) states that the difference between the flux in and out of an elemental volume per unit of time equals the change in storage of the soil system, ensuring the conservation of fluid mass during motion. For steady-state flow, the water storage portion (right side) of Equation (3) reduces to zero.

The analytical solutions to Richard’s equation can only be applied to a limited number of cases involving homogeneous soils, simplified initial and boundary conditions, and simple constitutive relationships that detail unsaturated hydraulic properties, mainly water content and permeability change with soil suction [4]. The non-linear differential equation can be linearized using various transformation methods such as Fourier’s or Laplace’s transformation methods. Gardner’s exponential model [44] which relates hydraulic conductivity with the pressure head, as shown in Equation (4), is also used for linearization.

$$k_u=k_s{e}^{-\alpha \left(u_a-u_w\right)}$$

(4)

where $k_u$ is the unsaturated hydraulic conductivity (m/s), $k_s$ = saturated hydraulic conductivity (m/s), α = inverse of air-entry pressure (kPa⁻¹), $u_a$ = pore air pressure (kPa), $u_w$ = pore water pressure (kPa), $\left(u_a-u_w\right)$ = matric suction (kPa).

Many downward infiltration analytical solutions have been proposed for Richard’s equation [45,46,47,48]. Simunek (2005) [4] observed that while there exist many analytical and semi-analytical solutions to unsaturated flow problems, some are widely used and many others appear to be only of academic interest. For slope stability problems, Lu and Griffiths’s (2004) [12] vertical steady infiltration solution has found wide use [49,50,51,52]. Additionally, only two pore size parameters, α and n, are required to account for unsaturated hydraulic behavior. This study adopts Lu and Griffiths’s 1-D analytical solution and combines it with 2-D slope stability analysis to quantify the contribution of suction to slope stability. Owing to the simplification of the infiltration process, the results can be used for an initial study.

2.1.3. Numerical Solution of Richard’s Equation

Numerical methods have the advantage of incorporating advanced models of unsaturated hydraulic behavior, complex geometries, initial and boundary conditions, and hydrological conditions during infiltration analysis [4]. An increase in computational technology has led to the development of computer programs based on finite element and finite difference codes. The FE seepage analysis procedure is illustrated in Figure 1. The soil mass is discretized into elements consisting of nodes located at the corners or midpoints of the element sides. For two-dimensional problems, triangular or quadrilateral elements are normally adopted. The main objective is to solve for the primary unknown (hydraulic head) which constitutes the primary driving potential for water flow in saturated and unsaturated soils. The FEM formulation for transient seepage is developed through spatial discretization of an element and applying the Galerkin principle of weighted residuals [53] while satisfying nodal compatibility which ensures that each node shared by adjacent elements has the same hydraulic head. A finite difference technique is used to approximate the time derivative relating to the nodal heads of an element at two successive time steps.

The element equations are assembled into a set of global flow equations that are solved on an iterative basis considering that the coefficient of permeability varies with matric suction (related to nodal hydraulic heads) in unsaturated soil systems. During the first iteration, the coefficient of permeability is approximated based on initial boundary conditions (hydraulic heads or water table) to evaluate the first set of nodal hydraulic heads. The calculated hydraulic heads are used to compute the average matric suction in the element. In the next iteration, the coefficient of permeability is adjusted based on the new set of matric suction values. This process is repeated until the hydraulic head and permeability differences are smaller than a specified tolerance between successive iterations or the maximum number of iterations is reached. Other quantities such as pore water pressures and hydraulic gradients can also be calculated [8].

The computer programs reviewed in the second part of this study, mainly Slide 2D/RS2, PLAXIS LE, PLAXIS 2D, and SEEP/W have a similar seepage/infiltration procedure as described above. Detailed explanations of the FE seepage mechanism can be found in the respective manuals [14,15,16,18,19].

2.2. Slope Stability Analysis in Soil Slopes

2.2.1. Coupled and Uncoupled Analysis

Two main methods are used to assess the physical response of unsaturated soils to rainfall infiltration, the coupled and uncoupled approachs. In the uncoupled approach, seepage analysis is first conducted to estimate pore water pressures which are used as input in deformation or stability analysis [54,55]. The accuracy of seepage solutions is dependent on the chosen time steps, hence small time increments are recommended for transient seepage analysis [8]. The coupled approach considers the interdependence between seepage and deformation by simultaneously solving pore water pressures and deformation during rainfall infiltration. This way, the coupled analysis method demonstrates a more realistic physical relation between water flow and stress-strain variation within unsaturated soils [56].

Table 1. Analysis methods in seepage and slope stability analysis.

Software	Analysis Method
RS2	Partially coupled
Slide 2D	Uncoupled
SEEP/W
PLAXIS LE
PLAXIS 2D	Fully coupled

summarizes the approaches implemented in the software programs used in this study. The partially coupled analysis in RS2 means that changes in pore pressures and effective stress can affect deformation whereas changes in deformation cannot affect pore pressures [14].

2.2.2. $FoS$ and Failure Surface Estimation

Slope analysis methods are mainly distinguished by the method used to locate the critical failure surface and the $FoS$ definition. The limit equilibrium method (LEM) based on static equilibrium assumes a number of slip surfaces along which sliding may occur until the critical slip surface with the lowest $FoS$ is determined. In unsaturated soil analysis, the additional shear strength due to suction is included in the shear strength mobilized at the base of the slices. Commercial software provides different slip surface search criteria such as the failure surface shape (circular or non-circular), the center of the failure circle, and depth of the failure surface, which can be manipulated by the user to obtain the desired result. For an accurate solution of the $FoS$, an iterative procedure with certain convergence criteria is adopted. The convergence criteria may be in form of maximum allowable changes in calculated $FoS$ for successive iterations or limitation on the magnitude of force and moment imbalance or a combination of both [57]. Commercial software programs PLAXIS LE, SLOPE/W, and Slide 2D are based on the LEM.

The finite element analysis (FEA) method locates the critical failure surface based on a calculation of failure shear strain zones which have been shown to coincide with rupture surfaces in experimental tests [58]. The $FoS$ is calculated following two methods: the strength reduction finite element analysis method (SRFEA) [59,60] and the finite element limit analysis method (FELA), based on upper and lower bound techniques [61,62]. SRFEA is used in PLAXIS 2D and RS2 [14,19] while FELA is employed in Optum G2 [13]. FELA has recently been applied to a wide range of geotechnical problems [63,64].

3. Stability Analysis of Unsaturated Slopes Subjected to 1-D Steady-State Infiltration—Part I

This section proposes guidelines on the contribution of apparent cohesion in unsaturated slopes subjected to vertical steady-state rainfall infiltration. As an extension of a study by Li et al. 2022 [11] which considered hydrostatic conditions, this parametric study considers the effect of environmental conditions, in the form of steady-state rainfall, on the long-term stability of unsaturated slopes. Lu and Griffiths’s (2004) [12] formulation based on the suction stress concept is used to calculate the vertical apparent cohesion profiles for different flow rates.

3.1. Suction Stress Concept

Lu et al. (2010) [65] proposed a closed-form expression for suction stress for the full range of matric suction, as seen in Equations (5) and (6).

$${\sigma }^s=-\left(u_a-u_w\right) u_a-u_w\le 0$$

(5)

$${\sigma }^s=-\frac{\left(u_a-u_w\right)}{{\left(1+{\left[\alpha \left(u_a-u_w\right)\right]}^n\right)}^{\left(n-1\right)/n}}{u}_a-u_w>0$$

(6)

where $\alpha$ (approximately the inverse of the air-entry value) and $n$ (related to the breadth of soil pore size) are fitting parameters used in Van Genuchten’s (1980) [66] soil water characteristic curve model.

Based on Darcy’s law [35] and Gardner’s model (1958) [44], Lu and Griffiths (2004) [12] developed a 1-D steady-state matric suction profile expression shown in Equation (7), which accounts for the factors such as soil type, environmental conditions, water table depth, and slope configuration.

$$\left(u_a-u_w\right)=\frac{-1}{\alpha } \ln \left(1+\raisebox{1ex}{$q$}\!\left/ \!\raisebox{-1ex}{$k_s$}\right.\right)e^{-\alpha {\gamma }_{w}z}-\raisebox{1ex}{$q$}\!\left/ \!\raisebox{-1ex}{$k_s$}\right.$$

(7)

where q is the steady infiltration rate (m/s) denoted as negative, ${\gamma }_w$ is the unit weight of water ($\raisebox{1ex}{$\mathrm{kN}$}\!\left/ \!\raisebox{-1ex}{${\mathrm{m}}^3$}\right.$), $\alpha$ is the inverse of the air-entry value (kPa⁻¹), $k_s$ is the saturated hydraulic conductivity (m/s), and z is the vertical distance from the water table (m). Positive upwards z = 0 at the water table.

The apparent cohesion of unsaturated soil, $C_{app}$, due to matric suction can be calculated by combining Equations (6) and (7), as presented in Equations (8) and (9).

$$C_{app}=-{\sigma }^s$$

(8)

$$C_{app}= \frac{\left(u_a-u_w\right)}{{\left(1+{\left[\alpha \left(u_a-u_w\right)\right]}^n\right)}^{\left(n-1\right)/n}}$$

(9)

where ${\varphi }^{\prime }$ (˚) is the soil friction angle.

3.2. Development of Numerical Model and Slope Parameters

The numerical model shown in Figure 2 is adopted in this study. Details on the boundary conditions, mesh and discretization, and the validation of the numerical model can be found in a previous study by Li et al. (2022) [11]. The Optum G2 software program was used for numerical analysis.

The soil and slope parameters are presented in

Table 2. Soil and slope parameters for unsaturated stability analysis.

Category	Parameter	Range
Mechanical properties	Effective cohesion, c′ (kPa)	1, 5
	$\mathrm{Internal} \mathrm{friction} \mathrm{angle}, {\varphi }^{\prime }$	25–35
	$\mathrm{Unit} \mathrm{weight}, \mathsf{\gamma } (\raisebox{1ex}{$\mathrm{kN}$}\!\left/ \!\raisebox{-1ex}{${\mathrm{m}}^3$}\right.$)	18
Hydraulic properties
Van Genuchten fitting parameters	α (kPa−1)	0.01, 0.1
	$n$	2–4
Hydrological properties
	$\raisebox{1ex}{$q$}\!\left/ \!\raisebox{-1ex}{$k_s$}\right.$(Infiltration)	0–(−0.9) *
Slope Geometry
	Slope height, H (m)	10–30
	Slope angle, β (˚)	25–40
	$\mathrm{Water} \mathrm{table} \mathrm{depth}, D_{sw}$(m)	5–30

* Infiltration is denoted as negative (−).

based on typical ranges for silty soil shown by Lu and Griffiths (2004) [12]. Li et al. 2022 [11] showed that most slopes with H ≤ 30 m and β = 25–40°, water table depths of 5–30 m are of engineering importance having $FoS$ ranging from 1.0–1.5. For such unsaturated slopes under hydrostatic conditions, the contribution of apparent cohesion to $FoS$ is considered significant for slope height, H = 10–30 m and water table depth ≥5 m [11]. Figure 3 categorizes different soil types based on pores size parameters α and $n$. Silty soil pore size parameters denoted α and $n$ range from (0.01–0.1) and (2–4), respectively. A sample with α = 0.01 and $n$ = 2 borders on silt/clay soil while with α = 0.1 and $n$ = 4 borders on silt/sand soil.

3.3. Examination of Key Parameters That Affect the Stability of Unsaturated Slopes

The criterion used in this study to quantify the contribution of apparent cohesion to slope stability is based on a comparison between the $FoS$ of cases that consider suction to cases without suction. If the $FoS$ difference (suction/without suction) is ≥5%, the contribution of apparent cohesion is deemed to be significant. For an $FoS$ difference (suction/without suction) of <5%, the apparent cohesion may be ignored [8]. Apparent cohesion profiles are sampled from the slope crest to the water table location. Preliminary findings revealed that the influence of apparent cohesion is mainly controlled by pore size parameters α and $n$, the hydrological loading, $\raisebox{1ex}{$q$}\!\left/ \!\raisebox{-1ex}{$k_s$}\right.$, the water table depth, $D_{sw}$, and the slope height, H, which are discussed in the next section. Other parameters, such as c′, ${\varphi }^{\prime }$ and β, primarily affect the $FoS$. Only one set of these parameters are presented (c′ = 5 kPa, ${\varphi }^{\prime }$= 25°, β = 30°) to ensure the influence of the main controlling parameters on apparent cohesion is captured.

3.3.1. Influence of Pore Size Parameters α and $n$ on Apparent Cohesion Profile

The α parameter is indicative of the pore size of a soil sample. Smaller α values signify that the soil has smaller pore sizes and a greater zone of water retention above the water table. Figure 4a shows apparent cohesion profiles for a light infiltration rate $\raisebox{1ex}{$q$}\!\left/ \!\raisebox{-1ex}{$k_s$}\right.$ = −0.3 with varying α. For a larger α = 0.1 value, bordering sandy soil, the maximum apparent cohesion value (3 kPa) occurs about 2 m from the water table. For a smaller α = 0.01 value, the maximum apparent cohesion is significantly larger (30 kPa) at about 17.5 m from the water table. The pore sizes determine the zone of influence of apparent cohesion above the water table. The n parameter, related to pore size distribution, has a less pronounced effect on the shape of the apparent cohesion profile, as seen in Figure 4b. The maximum apparent cohesion values occur at about the same depth, 3–5 m from the water table, and range from 6–7 kPa.

3.3.2. Influence of Hydrological Loading $\raisebox{1ex}{$q$}\!\left/ \!\raisebox{-1ex}{$k_s$}\right.$ on the Apparent Cohesion Profile

This study categorizes the infiltration rate into light $\raisebox{1ex}{$q$}\!\left/ \!\raisebox{-1ex}{$k_s$}\right.$ = 0–(−0.3), moderate $\raisebox{1ex}{$q$}\!\left/ \!\raisebox{-1ex}{$k_s$}\right.$ = −0.3–(−0.6), and heavy infiltration $\raisebox{1ex}{$q$}\!\left/ \!\raisebox{-1ex}{$k_s$}\right.$ = −0.6–(−0.9). For borderline silt/sand soil (α = 0.1, n = 4), a moderate infiltration rate can increase the apparent cohesion due to an increase in the water contact area between soil particles. A heavy infiltration rate ($\raisebox{1ex}{$q$}\!\left/ \!\raisebox{-1ex}{$k_s$}\right.=$ −0.7, −0.9) reduces the apparent cohesion significantly, as can be seen in Figure 5a. For borderline silt/clay soil (α = 0.01, n = 2), Figure 5b, an increase in infiltration rate $\raisebox{1ex}{$q$}\!\left/ \!\raisebox{-1ex}{$k_s$}\right.$ corresponds to a fall in apparent cohesion. Due to the small pore size and larger zone of water retention, the borderline silt/clay soil still possesses considerable apparent cohesion (5–15 kPa) at a heavy infiltration rate.

Figure 6a,b present failure surfaces based on shear strain for a silt/clay slope with H = 10 m, $D_{sw}$ = 10 m, subjected to light and heavy infiltration rates, $\raisebox{1ex}{$q$}\!\left/ \!\raisebox{-1ex}{$k_s$}\right.$, respectively. The increase in infiltration rate results in a shallower failure surface due to the loss of matric suction and wetting of the failure zone. As seen in Figure 6c, at $\raisebox{1ex}{$q$}\!\left/ \!\raisebox{-1ex}{$k_s$}\right.$ = −0.9, the effective degree of saturation, $S_e$ is close to one which translates to near full saturation of the unsaturated layer.

3.3.3. Influence of Water Table Location on Apparent Cohesion and $FoS$

The location of the water table determines the depth of the unsaturated zone and the failure surface propagation. Three locations of the water table are presented: above the slope toe, at the toe, and below the toe. The slope parameters include H = 10 m, α = 0.1, 0.01, n = 2, 4, and infiltration $\raisebox{1ex}{$q$}\!\left/ \!\raisebox{-1ex}{$k_s$}\right.$ = 0–(−0.9). Figure 7a,b show that the contribution of apparent cohesion for the studied parameters is relevant for water table depth ≥10 m, since $FoS$% difference ≥5%. For the borderline silt/sand slope, Figure 7a, the apparent cohesion is relevant under moderate infiltration conditions, $\raisebox{1ex}{$q$}\!\left/ \!\raisebox{-1ex}{$k_s$}\right.$ = −0.3–(−0.6). The $FoS$% difference between suction and no suction condition for $D_{sw}$ = 10 m ranges between 5–8%. For borderline silt/clay slopes, Figure 7b, the contribution of apparent cohesion is more significant for a wider range of $\raisebox{1ex}{$q$}\!\left/ \!\raisebox{-1ex}{$k_s$}\right.$ = 0–(−0.7), with the $FoS$% difference between suction and no suction ($D_{sw}$ = 10 m) ranging from 9–23%.

The failure surfaces based on shear strain shown in Figure 8a–c show the influence of the water table location on failure surface propagation. The slope properties include H = 10 m, α = 0.1, n = 4, $D_{sw}$ = 5 m, 10 m, 15 m, and infiltration $\raisebox{1ex}{$q$}\!\left/ \!\raisebox{-1ex}{$k_s$}\right.$ = −0.3. The failure mechanism for all cases is predominantly a toe failure. For a water table depth of 5 m, Figure 8a, the failure surface cuts slightly beneath the toe. At a water table depth of 10 m, Figure 8b, the path of the failure surface is similar but slightly shallower. The failure surface of a water table depth of 15 m in Figure 8c is shallower and does not extend beyond the toe due to the increased apparent cohesion around the toe area (See Figure 8d). The strengthened soil pushes the failure surface upwards.

3.3.4. Influence of the Slope Height on the Contribution of Apparent Cohesion and $FoS$

As a general rule for all slope heights sampled, the contribution of apparent cohesion is more significant for water table location below the mid-height of the slope due to the strengthening of the soil by suction forces. Thus, for larger slope heights, H > 10 m, the water table should be deeper for the contribution of apparent cohesion to be significant. Notably, for H = 20 m and 30 m slopes, only borderline silt/lay soils (α = 0.01, $n$ = 2–4) with a significant amount of apparent cohesion yield a $FoS$% difference >5% between suction and no suction cases. In Figure 9a, a silty slope with parameters H = 20 m, $\raisebox{1ex}{$q$}\!\left/ \!\raisebox{-1ex}{$k_s$}\right.$ = 0–(−0.9) is presented. The minimum water table depth to consider apparent cohesion is 14 m for $\raisebox{1ex}{$q$}\!\left/ \!\raisebox{-1ex}{$k_s$}\right.$ = 0–0.5. Similarly, for a slope with H = 30 m, in Figure 9b, the minimum water table depth to consider apparent cohesion is about 20 m.

Figure 10 shows the sensitivity of $FoS$ to input parameters in this study. The baseline parameters for the sensitivity analysis include H = 20 m, β = 25°, $D_{sw}$ = 10 m, c′ = 5 kPa, ${\varphi }^{\prime }$ = 25°, α = 0.1, $n$ = 4, $\raisebox{1ex}{$q$}\!\left/ \!\raisebox{-1ex}{$k_s$}\right.$ = −0.1. Each input parameter was changed by a certain percentage while holding others constant, then the slope $FoS$ was calculated and compared to the baseline $FoS$ in percentage form. The steeper the curve, the more influence the parameter has on the $FoS$. As demonstrated in Figure 10, for slopes with H ≥ 20 m and a fairly high water table (above the mid-slope height, $\raisebox{1ex}{$H_w$}\!\left/ \!\raisebox{-1ex}{$H$}\right.\ge$ 0.5), the slope height, H is a primary controlling factor of the $FoS$. Other parameters, with exception of β and ${\varphi }^{\prime }$, have a limited effect on the $FoS$.

3.4. Identification of Suitable Conditions for Consideration of Apparent Cohesion

This section identifies several slope and environmental parameters whereby the contribution of apparent cohesion to slope stability is significantly based on the criterion that there is an $FoS$% difference of ≥5% between suction and no suction cases. These parameters are summarized in

Table 3. Guidelines on consideration of apparent cohesion for H = 10 m.

Slope Height (m)	10
$\raisebox{1ex}{$q$}\!\left/ \!\raisebox{-1ex}{$k_s$}\right.$	0.3–0.6	0–0.7
$\mathsf{\alpha }\left(kP{a}^{-1}\right)$	0.1	0.01
n	2–4	2–4
Soil Type	Sandy silt	Silt/Silty Clay
$\mathrm{Minimum} \mathrm{Depth} \mathrm{of} \mathrm{Water} \mathrm{Table} D_{sw}$(m)	10	10
$\left(\raisebox{1ex}{$H_w$}\!\left/ \!\raisebox{-1ex}{$H$}\right.\right)$	(0.5)	(0.5)

,

Table 4. Guidelines on consideration of apparent cohesion for H = 20 m.

Slope Height (m)	20
$\raisebox{1ex}{$q$}\!\left/ \!\raisebox{-1ex}{$k_s$}\right.$	0–0.5	0–0.7
$\mathsf{\alpha }\left(kP{a}^{-1}\right)$	0.01	0.01
n	2–4	2–4
Soil Type	Silt/Silty Clay	Silt/Silty Clay
$\mathrm{Minimum} \mathrm{Depth} \mathrm{of} \mathrm{Water} \mathrm{Table} D_{sw}$(m)	14	20
$\left(\raisebox{1ex}{$H_w$}\!\left/ \!\raisebox{-1ex}{$H$}\right.\right)$	(0.3)	(0)

and

Table 5. Guidelines on consideration of apparent cohesion for H = 30 m.

Slope Height (m)	30
$\raisebox{1ex}{$q$}\!\left/ \!\raisebox{-1ex}{$k_s$}\right.$	0–0.5	0–0.7
$\mathsf{\alpha }\left(kP{a}^{-1}\right)$	0.01	0.01
n	2–4	2–4
Soil Type	Silt/Silty Clay	Silt/Silty Clay
$\mathrm{Minimum} \mathrm{Depth} \mathrm{of} \mathrm{Water} \mathrm{Table} D_{sw}$(m)	20	30
$\left(\raisebox{1ex}{$H_w$}\!\left/ \!\raisebox{-1ex}{$H$}\right.\right)$	(0.3)	(0)

. For the practicing engineer, the first step to use the guidelines is to identify the slope parameters, H and $D_{sw}$. Secondly, the soil sample is classified using α and $n$ pore size parameters. Figure 3 can be used for reference. The infiltration rate, $\raisebox{1ex}{$q$}\!\left/ \!\raisebox{-1ex}{$k_s$}\right.$, is then estimated to check if the slope case falls within the specified conditions. Since these guidelines are exclusive to the parameters chosen for this study, proper engineering judgment should be applied on what conditions to consider the contribution of apparent cohesion to slope strength as part of an initial study.

In Table 3, several conditions are identified where the contribution of apparent cohesion may be considered for an H = 10 m slope. For a borderline sand/silt slope (α = 0.1, $n$ = 2–4), apparent cohesion can be considered in stability analysis for moderate infiltration rate conditions and water table depth, $D_{sw}$ ≥ 10 m. Outside these parameter boundaries, the contribution of apparent cohesion is insignificant. For a borderline silt/clay slope (α = 0.01, $n$ = 2–4), apparent cohesion may be considered for a wider range of infiltration rates due to larger suction forces and an expanded zone of water retention.

Table 4 and Table 5 present parameter boundaries for H = 20 m and 30 m slopes. The predominant soil types are silt/clay (α = 0.01, $n$ = 2–4), which possess appreciable apparent cohesion in the unsaturated layer. It is noted that the minimum depth of the water table to consider apparent cohesion varies from the lower 1/3 H of the slope to the slope toe ($\raisebox{1ex}{$H_w$}\!\left/ \!\raisebox{-1ex}{$H$}\right.$ = 0.3–0) depending on the applied infiltration rate. For H = 20 m and 30 m slopes with the water table at the slope toe ($\raisebox{1ex}{$H_w$}\!\left/ \!\raisebox{-1ex}{$H$}\right.$ = 0), the contribution of apparent cohesion remains considerable over a wide range of infiltration rates, $\raisebox{1ex}{$q$}\!\left/ \!\raisebox{-1ex}{$k_s$}\right.$ 0–0.7. At a heavy infiltration rate, $\raisebox{1ex}{$q$}\!\left/ \!\raisebox{-1ex}{$k_s$}\right.>$ 0.7, most of the suction is dissipated and apparent cohesion is insignificant.

4. Description and Numerical Analysis of Case Studies—Part II

4.1. Case 1—Assumed Case

4.1.1. Description of Case Study

This slope case is derived from a parametric study carried out by Rahardjo et al. (2007) [9] to highlight the significance of different controlling factors on the instability of homogeneous slopes subjected to rainfall infiltration. Transient seepage and slope stability analyses are conducted using different FEM programs to demonstrate their ability to predict the hydrological response of the unsaturated slope subjected to rainfall. The pore water pressure profiles, water table locations, and $FoS$ variations with time are then checked against the reported results. The slope model and boundary conditions are illustrated in Figure 11.

Table 6. Soil mechanical and hydraulic properties.

Category	Parameter	Silty Soil
Mechanical properties (Mohr coulomb model)	Effective cohesion, c′ (kPa)	10
	$\mathrm{Unit} \mathrm{weight}, \mathsf{\gamma } (\raisebox{1ex}{$\mathrm{kN}$}\!\left/ \!\raisebox{-1ex}{${\mathrm{m}}^3$}\right.$)	20
	$\mathrm{Effective} \mathrm{angle} \mathrm{of} \mathrm{internal} \mathrm{friction}, {\varphi }^{\prime }$	26
	$\mathrm{Matric} \mathrm{suction} \mathrm{angle}, {\varphi }^b$	26
	Young’s modulus, E′ (kPa) a	20,000
	Poisson’s ratio, 𝜐 a	0.3
Hydraulic properties
Permeability function	$\mathrm{Saturated} \mathrm{coefficient} \mathrm{of} \mathrm{permeability}, k_s$(m/s)	10−5
	$p$	4
SWCC fitting parameters	$a_f$(kPa)	50
	$n_f$	1
	$m_f$	1

^a Assumed stiffness values used in RS2 and PLAXIS 2D.

shows the soil shear strength and hydraulic parameters adopted based on typical values encountered in a few Singapore sites. A rainfall flux (36 mm/h) was applied on the slope surface for a one-day period.

The SWCC and permeability fitting functions by Fredlund and Xing (1994) [67] and Leong and Rahardjo (1997) [68] in Equations (10) and (12) were used, respectively.

$${\theta }_w={\theta }_{s}C\left(\psi \right) \left\{\frac{1}{{\left\{ln\left[e+{\left(\frac{\left(u_a-u_w\right)}{a_f}\right)}^{n_f}\right]\right\}}^{m_f}}\right\}$$

(10)

The correction factor, $C\left(\psi \right)$, is defined in Equation (11) as

$$C\left(\psi \right)=1-\frac{ln\left(1+\raisebox{1ex}{$\left(u_a-u_w\right)$}\!\left/ \!\raisebox{-1ex}{${\left(u_a-u_w\right)}_r$}\right.\right)}{ln\left[1+\left(\raisebox{1ex}{${10}^6$}\!\left/ \!\raisebox{-1ex}{${\left(u_a-u_w\right)}_r$}\right.\right)\right]}$$

(11)

where ${\theta }_w$ = volumetric water content, ${\theta }_s$ = saturated volumetric water content, $a_f$ = fitting parameter related to the air-entry value, $n_f$ = fitting parameter related to the rate of soil desaturation, $m_f$ = fitting parameter related to the residual water content, $e$ = natural number, $\left(u_a-u_w\right)$ = matric suction (kPa), $u_a$ = pore air pressure, $u_w$ = pore water pressure, and ${\left(u_a-u_w\right)}_r$ = residual matric suction corresponding to residual water content.

The permeability function,

$$k_w=k_s{\Theta }^p$$

(12)

where $k_w$ = coefficient of permeability with respect to water for unsaturated soil, $k_s$ = saturated coefficient of permeability, $p$ = fitting parameter related to the slope of permeability function, and $\Theta$ = $({\theta }_w/{\theta }_s)$ = dimensionless form of the SWCC.

Rainfall infiltration analysis comprised two stages: (1) the seepage analysis and (2) the stability analysis. Two-dimensional transient water flow due to rainfall follows Richard’s governing differential equation [43], which was solved using SEEP/W software (Geo-studio 1998) [69].

For the initial conditions, a hydrostatic condition based on the initial water table location was assumed. A limiting negative pore pressure cutoff of −75 kPa was imposed to prevent the generation of unrealistic pore pressures. The pore water pressures obtained from the seepage analysis were then incorporated into the stability analysis to calculate the factor of safety ($FoS$). Fredlund et al. (1978) [70] used the unsaturated shear strength expression, Equation (13), to account for the contribution of negative pore pressure to the soil shear strength.

$${\tau }_f={\mathrm{c}}^{\prime }+\left({\sigma }_n-u_a\right)tan{\varphi }^{\prime }+\left(u_a-u_w\right)tan{\varphi }^b$$

(13)

where ${\tau }_f$ = unsaturated shear strength of soil, $c^{\prime }$ = effective cohesion, $\left({\sigma }_n-u_a\right)$ = net normal stress, ${\sigma }_n$ = total normal stress, $u_a$ = pore air pressure, $\left(u_a-u_w\right)$ = matric suction, $u_w$ = pore water pressure, and ${\varphi }^b$ = matric suction angle.

SLOPE/W (Geo-studio 1998) [71] was used for slope stability analysis. Bishop’s simplified method based on the LEM was adopted to calculate the $FoS$.

4.1.2. Seepage and Stability Modeling

FE Seepage Modeling

The assumed case is modeled using five different software programs: Slide2D, RS2, PLAXIS LE, Geo-studio (SEEP/W, SLOPE/W), and PLAXIS 2D. The first stage of rainfall infiltration analysis is the seepage analysis. Slide 2D, RS2, PLAXIS LE, and PLAXIS 2D have in-built FE seepage modules to run seepage problems. SLOPE/W has a companion software SEEP/W that conducts FE seepage analysis. In this study, the parallel seepage and stability analyses using SEEP/W and SLOPE/W, respectively, serve as a reference point for suction values and failure surface propagation which are not presented in the original paper [9].

The mesh and discretization details are summarized in

Table 7. Mesh and discretization details for FEM software.

Software	Element Type	No. of Elements/Nodes	Mesh Element Length (m)
			Region 1	Region 2
RS2	6-noded triangle	2725/5618	-	-
Slide 2D	6-noded triangle	1987/4106	-	-
PLAXIS 2D	15-noded triangle	2463/20031	1.0	3.0
PLAXIS LE	3-noded triangle	3970/2135	0.5	1.0
SEEP/W	4-noded quadrilaterals and 3-noded triangles	2451/2557	0.5	1.0

. With reference to Figure 11, the slope model is divided into two regions. The mesh in Region 1 (above the water table) is refined to enhance the accuracy of the infiltration process.

Table 8. Water content and permeability functions in FE seepage analysis.

Software	SWCC Model	Permeability Model
RS2	Fredlund and Xing (1994)	Fredlund and Xing (1994)
Slide 2D
SEEP/W
PLAXIS LE		Leong and Rahardjo (1997)
PLAXIS 2D	Van Genuchten (1980)	Van Genuchten (1980)

shows the water content and permeability functions used by each software for seepage analysis. Notably, Leong and Rahardjo’s (1997) [68] permeability model is only integrated into PLAXIS LE.

Stability Modeling

For LEM-based programs PLAXIS LE, Slide 2D, and SLOPE/W, slope stability analysis is performed using Bishop’s simplified method. The slip surface search criteria and $FoS$ calculation methods are shown in

Table 9. Slip surface search criteria and $FoS$ calculation methods.

Software	$\mathit{F}\mathit{o}\mathit{S}\mathbf{Definition}$	Slip Surface Search Method
RS2	Strength reduction method	Finite Element Analysis (Shear strain)
PLAXIS 2D
Slide 2D	LE—Bishop’s simplified method	Auto-refine search [15]
SLOPE/W		Entry and exit points [17]
PLAXIS LE		Slope search [72]

. FE-based programs, RS2 and PLAXIS 2D, calculate $FoS$ using the strength reduction method.

Table 10. Unsaturated shear strength equations used in Case 1.

Reference	Unsaturated Shear Strength Equation	Commercial Software
Bishop (1959) [73]	${\tau }_f={\mathrm{c}}^{\prime }+\left(\sigma -u_a\right)tan{\varphi }^{\prime }+\chi \left(u_a-u_w\right)tan{\varphi }^{\prime }$	PLAXIS 2D
Fredlund et al. (1978) [70]	${\tau }_f={\mathrm{c}}^{\prime }+\left({\sigma }_n-u_a\right)tan{\varphi }^{\prime }+\left(u_a-u_w\right)tan{\varphi }^b$	PLAXIS LE, RS2, Slide 2D, SLOPE/W

Note: $\chi =$ effective stress coefficient which varies from 0–1 as a function of the degree of pore-water saturation, $S_e$ = ($S-S_{res}/1-S_{res}$), where $S$ = degree of saturation, and $S_{res}$ = residual degree of saturation.

provides the unsaturated shear strength equations employed by different software for this case study.

4.1.3. Comparison of Results from Geotechnical Software

The suction profile in Figure 12a represents hydrostatic conditions drawn from an initial water table location. The suction cut-off, applied at the initial stage of seepage analysis to prevent the development of excess negative pore pressures, is visible in the SEEP/W and PLAXIS LE pore pressure profiles. In PLAXIS 2D, it is not possible to specify a suction cut-off value. In RS2 and Slide 2D, even though the suction cut-off is not displayed due to software limitations, it is still implemented with regard to the calculation of effective stress and unsaturated shear strength. Additionally, a steady-state groundwater analysis is applied at the initial stage in RS2 and Slide 2D which slightly alters the initial water table location compared to the other software programs as seen in Figure 12a.

The suction profiles presented in Figure 12a–f demonstrate progressive seepage with rainfall infiltration. After 6 h of rainfall, the PLAXIS 2D pore pressure profile in Figure 12b shows a well-defined and rapidly progressing wetting front compared to other computer programs, with the upper 5 m of the slope having pore pressure close to 0 kPa, resulting in critical conditions being reached earlier [54]. The PLAXIS 2D model collapses after 8 h of rainfall and analysis abruptly stops. Further results are not available.

At the 12 h and 24 h mark, pore pressure values fall below −40 and −20 kPa, respectively. Generally, RS2 and Slide 2D use the same finite element engine for seepage analysis which explains the close similarity of the pore water pressure profiles [14]. In Figure 12c,d, it can be deduced that SEEP/W predicts faster dissipation of matric suction in the upper 1/3 H of the slope as a result of rainfall infiltration. This phenomenon is most visible around the slope crest after 12 h of rainfall. At 24 h, the SEEP/W pore water pressure at the slope crest reduces to 0 kPa due to the accumulation of a film of rainwater on the slope surface.

In the mid 1/3 H of the slope, PLAXIS LE suction values are slightly larger by 25–50% (around 5–10 kPa) compared to the other three software programs as seen in Figure 12c,d. Pore water pressure values for RS2 and Slide 2D for the bottom 1/3 H of the slope are mostly above 0 kPa, signifying a rise of the water table. Figure 13 presents the water table location predicted by each computer program after 24 h of rainfall. The water table rise in PLAXIS LE is minimal in contrast to RS2, Slide 2D, and SEEP/W, whereby the groundwater surges to the mid-slope height range along the sloping surface. This phenomenon may have arisen due to the element types used in each software as shown in Table 7. Generally, higher-noded elements (triangular, quadrilateral, or both) provide a more robust and accurate solution due to a higher number of integration points during the numerical integration process in FE seepage analysis [74]. However, further investigation is needed. The post-rainfall period is characterized by a rapid recovery rate of the matric suction and almost identical pore water pressure profiles for each software program.

The $FoS$ variation with time is presented in Figure 14. Apart from PLAXIS 2D, other software programs record the minimum $FoS$ after 24 h of rainfall, in agreement with Rahardjo’s et al. (2007) results. PLAXIS 2D based on a fully coupled formulation yields a lower critical $FoS$ (<1) at the 6th–8th h mark. The slope model collapses shortly afterward. The minimum $FoS$ for PLAXIS LE is about 15% larger than Rahardjo’s result. This is attributed to the bigger suction values recorded in PLAXIS LE during seepage analysis as seen in Figure 12b,c. Although RS2 and Slide 2D yield similar pore pressure profiles, there is a notable difference in calculated $FoS$. This is most pronounced at the initial and closing stages where the $FoS$ difference is about 30%.

Figure 15a–f show failure surfaces generated by PLAXIS 2D, RS2 (maximum shear strain contours), PLAXIS LE, SLOPE/W, and Slide 2D. The failure surface by RS2 cuts deeper and extends from the slope toe to the back of the slope crest compared to Slide 2D at the initial and end stages. Apart from the disparity in the method used to locate the critical failure surface between RS2 and Slide 2D, the two computer programs consider tension in the slope quite differently. In RS2, elements can fail in shear or tension. It is assumed that if tension failure occurs in a material before shear failure, the shear strength is reduced to the residual shear strength parameters for that material [14]. However, in this case, perfectly plastic material properties were assumed. The residual shear strength parameters equal peak shear strength parameters and tension failure do not affect shear strength properties. In Slide 2D, the effect of tension is included by the use of tension cracks or discarding surfaces with tension [15]. This study considered the latter option by setting tensile strength to 0 kPa, which adjusts the local $FoS$ on a slice such that the effective normal stress at the base of a slice is zero. The overall effect could be a larger calculated $FoS$ in RS2 than in Slide 2D.

For verification purposes, RS2 pore water pressure contours were manually imported into the Optum G2 software program, and converted to total cohesion, c, following Equation (14).

$$\mathrm{C}={\mathrm{c}}^{\prime }+\left(u_a-u_w\right)tan{\varphi }^b$$

(14)

A slope stability analysis was conducted and the results are summarized in

Table 11. $FoS$ comparison between RS2 and Optum G2.

Time	$\mathbf{Factor} \mathbf{of} \mathbf{Safety} \mathit{F}\mathit{o}\mathit{S}$
	Rahardjo	RS2	Optum G2 *
Initial	2.13	2.53 (19%)	2.05 (−4%)
12 h	1.40	1.67 (19%)	1.40 (0%)
24 h	1.03	1.08 (5%)	1.03 (0%)
48 h	1.23	1.41 (15%)	1.28 (4%)
240 h	1.89	2.31 (22%)	1.89 (0%)

* The $FoS$ is the average of upper and lower bound results.

. Optum G2 yields very close results to Rahardjo’s, with an $FoS$% difference (in brackets) of ≤4%. Shear strain contour diagrams in Figure 16a–d show that the failure surface in Optum G2 is shallower at the initial and end stages compared to RS2 despite having similar pore water pressure distribution. This implies that the main difference may lie in the failure mode obtained in the RS2 software. At critical conditions, Figure 16c, both RS2 and Optum G2 yield similar $FoS$ and failure surface.

The failure surfaces shown in Figure 15a–d become shallower with the ingress of rainfall, a common attribute of rainfall-induced landslides [7]. PLAXIS 2D and RS2 failure surfaces are based on shear strain. From previous results, PLAXIS 2D predicts slope failure after 6 h of rainfall. Correspondingly, the failure surface in Figure 15b is relatively shallower than in other software programs. The LEM circular failure surfaces are mainly influenced by the specified search criteria and the suction values along the failure surface. Initially, the SLOPE/W failure surface is deeper than other LEM-based software programs since it matches the manually specified entry and exit slip surface points. During rainfall, Figure 15b–d, the SLOPE/W failure surface quickly decreases in depth due to the rapid loss of matric suction highlighted earlier. After 24 h, failure surfaces of SLOPE/W and Slide 2D seemingly coincide. PLAXIS LE failure surface cuts slightly deeper due to increased cohesion brought about by slightly larger suction values around the failure surface zone. After 24 h, all the software programs have similar failure surfaces at critical conditions. Post-rainfall, the failure surface shape in the software programs recovers to resemble that of initial conditions due to a similarity in suction values. See Figure 15e,f.

4.2. Case 2—Full-Scale Landslide Flume Test

4.2.1. Description of Case Study

The second case involves a full-scale flume test conducted by Lee et al. (2019) [20] to evaluate the failure mechanisms of a rainfall-induced landslide. The large-scale nature of the experiment replicated the actual size of landslides and enabled the reproducibility of the landslide phenomenon by minimizing scaling effects [75]. The flume consisted of a rectangular channel that was 20 m long, 4 m wide, and 2.5 m deep, with four sections, named 1–4. The lower slope (1) was a 5 m-long section for which the inclination angle was 5°. The mid-upper slope (2) was a 5 m-long section that was inclined to 15°. The main slope (3) was a 10 m-long section that was inclined to 35°. The occurrence angle of 60% of landslides in The Republic of Korea ranges from 30–40° [76]. Globally, landslides occur on slopes with inclination angles of 16–44° [77,78]. The uppermost section (4) was fixed and could not be inclined. Volumetric water content and matric suction sensors were installed at various depths along the slope. Figure 17 is a schematic representation of the slope model and sensor locations.

Weathered granite soil, classified as poorly graded sand (SP) according to the Unified Soil Classification System (USCS), was used in the large-scale experiment.

Table 12. Soil parameters reported in the experiment.

Category	Parameter	Weathered Granite
Physical properties	Effective cohesion, c′ (kPa)	2.5
	$\mathrm{Unit} \mathrm{weight}, \mathsf{\gamma } (\raisebox{1ex}{$\mathrm{kN}$}\!\left/ \!\raisebox{-1ex}{${\mathrm{m}}^3$}\right.$)	18
	$\mathrm{Effective} \mathrm{angle} \mathrm{of} \mathrm{internal} \mathrm{friction}, {\varphi }^{\prime }$	30
	Optimum moisture content (%)	21.7
	USCS	SP
Hydraulic properties	$\mathrm{Saturated} \mathrm{coefficient} \mathrm{of} \mathrm{permeability}, k_s$(m/s)	4.1 × 10−6
	$\mathrm{Saturated} \mathrm{volumetric} \mathrm{water} \mathrm{content}, {\theta }_s$(%)	37.75
Van Genuchten fitting parameters	$\mathrm{Residual} \mathrm{volumetric} \mathrm{water} \mathrm{content}, {\theta }_r$(%)	3.75
	α, n, m	0.058, 2.7, 1.15

shows the basic physical and hydraulic soil properties based on direct shear, field density, permeability, and sieve analysis tests. The laboratory data on water content and permeability variation with matric suction were best fitted using the Van Genuchten (1980) [66] and Mualem (1976) [79] models shown in Equations (15) and (16), respectively. See Figure 18a,b.

$${\theta }_w={\theta }_r+\left({\theta }_s-{\theta }_r\right) {\left[\frac{1}{1+{\left(\alpha \Psi \right)}^n}\right]}^m$$

(15)

$$k_u=k_s\frac{{\left\{1-{\left(\alpha \psi \right)}^{n-1}{\left[1+{\left(\alpha \psi \right)}^n\right]}^{-m}\right\}}^2}{{\left[1+{\left(\alpha \psi \right)}^n\right]}^{\raisebox{1ex}{$m$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}$$

(16)

where α, $n$ and $m$ are fitting parameters, $\psi$ is the matric suction, and $k_u$ is the unsaturated coefficient of permeability.

Artificial rainfall was applied in two stages. The first stage was preceding rainfall with an intensity of 30 mm/h for a period of 3 h in order to increase the water content of the soil. The second stage consisted of extreme rainfall 1 and 2 applied at an intensity of 70 mm/h for a period of 24 h. On the application of extreme rainfall 2 (19 October 1000 h), surface run-off resulted in scouring at the bottom which progressed toward the top of the slope through a shallow fracture. This lowered the stability of the upper part of the slope resulting in failure on 19 October 1935 h, in the form of a circular arc at the top of the slope.

4.2.2. Development of Numerical Model

The full-scale landslide test was analyzed using three commercial software programs: SEEP/W (SLOPE/W), RS2, and PLAXIS 2D, based on uncoupled, partially coupled, and fully coupled formulations, respectively. The numerical model used and applied the boundary conditions illustrated in Figure 19. The Mohr–Coulomb constitutive model under drained conditions was used to capture the behavioral response of weathered granite soil to rainfall conditions since, under rainfall infiltration, excess pore water pressures could be fully dissipated [80].

The steel flume was modeled as an impermeable bedrock of infinite strength, similar to Yang et al. (2021) [24]. The mesh discretization settings are similar to the first case study. The upper layer of weathered granite was refined to increase the accuracy of the transient infiltration process and failure surface propagation. In RS2 and PLAXIS 2D, standard fixities were applied to simulate displacement boundary conditions. Displacement at the bottom was fully restrained in the horizontal and vertical directions (x- and y-direction) using fixed supports. At the sides, roller supports only allowed vertical movement (y-direction) while the slope surface was unrestrained. For flow conditions, both sides of the model were defined as seepage boundaries to model the free flow of water. The bottom was defined as a no-flow boundary ($q=0 \raisebox{1ex}{${\mathrm{m}}^3$}\!\left/ \!\raisebox{-1ex}{$\mathrm{s}$}\right.$) and the input infiltration flux was applied on the surface to model rainfall infiltration.

To account for unsaturated shear strength in this case study, the equations shown in

Table 13. Unsaturated shear strength equations used in Case 2.

Reference	Unsaturated Shear Strength Equation	Commercial Software
Bishop (1959) [73]	${\tau }_f={\mathrm{c}}^{\prime }+\left(\sigma -u_a\right)tan{\varphi }^{\prime }+\chi \left(u_a-u_w\right)tan{\varphi }^{\prime }$	PLAXIS 2D
Fredlund et al. (1978) [70]	${\tau }_f={\mathrm{c}}^{\prime }+\left({\sigma }_n-u_a\right)tan{\varphi }^{\prime }+\left(u_a-u_w\right)tan{\varphi }^b$	RS2
Vanapalli et al. (1996) [81]	$\tau ={\mathrm{c}}^{\prime }+\left(\sigma -u_a\right)tan{\varphi }^{\prime }+\left[\left(\frac{{\theta }_w-{\theta }_r}{{\theta }_s-{\theta }_r}\right)tan{\varphi }^{\prime }\right]$	SLOPE/W

are used.

${\varphi }^b$ was assumed to be 15° (=0.5 ${\varphi }^{\prime }$). The SWCC and permeability function was defined using the Van Genuchten (1980) [66] and Mualem (1976) [79] models integrated into each software program.

The c′ and ks experimental values for weathered granite soil obtained through direct shear and constant head permeability tests were 2.5 kPa and 4.1 × 10⁻⁶ m/s, respectively. The c′ and ks used in the numerical analysis were calibrated through a trial and error approach to match the experimental failure time and failure mechanism and obtained a reasonable FoS vs. time variation. For the partially coupled and fully coupled numerical models in RS2 and PLAXIS 2D, respectively, the best-fit c′ and ks parameters were found to be 2.6 kPa and 3.0 × 10⁻⁶ m/s. It should be noted that it is difficult to measure exact values of ks due to variations in sample size and collection techniques, measuring techniques, flow rate, and system dimensions. Repeated constant and falling head permeability tests, particularly on sandy soils, may yield results that vary significantly (coefficient of variation, CoV = 40–100%). Thus, the chosen numerical ks is within a reasonable range [82,83]. For the partially coupled and fully coupled numerical models in this study, the ks is the main controlling parameter of the hydro-mechanical response and failure time of the slope. Changes in c′ only affected the calculated FoS values. The uncoupled model (SEEP/W, SLOPE/W) is not sensitive to changes in ks for the range of values tested. Spencer’s method which satisfies both moment and force equilibrium was used for stability analysis in SLOPE/W.

A total of 15 stress points were selected to monitor the slope hydrological response to applied rainfall. To obtain similar initial water content conditions to the experiment, a small rainfall flux was applied at the slope surface for a prolonged period [84]. The parameters used in the numerical analysis are shown in

Table 14. Soil parameters used in the numerical analysis.

Category	Parameter	Weathered Granite	Steel Flume
Physical properties	Effective cohesion, c′ (kPa)	2.6	-
	$\mathrm{Internal} \mathrm{friction} \mathrm{angle}, {\varphi }^{\prime }$	30	-
	$\mathrm{Unit} \mathrm{weight}, \mathsf{\gamma } (\raisebox{1ex}{$\mathrm{kN}$}\!\left/ \!\raisebox{-1ex}{${\mathrm{m}}^3$}\right.$)	18 a	25
	Stiffness, E′ (kPa)	15 × 103 b	200 × 106
	Poisson’s ratio, 𝜐	0.3 b	0.2
Hydraulic properties	$\mathrm{Saturated} \mathrm{coefficient} \mathrm{of} \mathrm{permeability}, k_s$(m/s)	3.0 × 10−6	-
	$\mathrm{Saturated} \mathrm{volumetric} \mathrm{water} \mathrm{content}, {\theta }_s$(%)	37.75	-
	$\mathrm{Residual} \mathrm{volumetric} \mathrm{water} \mathrm{content}, {\theta }_r$(%)	3.75	-
Van Genuchten fitting parameters	α, n, m	0.058, 2.7, 1.15	-

^a PLAXIS 2D γ (unsaturated) = 18.57, γ (saturated) = 19.58, ^b assumed values.

.

4.2.3. Comparison of Numerical Results with Experimental Results

Volumetric Water Content vs. Time

The numerical model was validated by checking the volumetric water content variation with time, failure surface propagation, and failure time, against the reported experimental results. Lee et al. 2019 [20] mentioned that the matric suction sensors had limited capacity and could not measure values below 9 kPa. The suction vs. time experiment profiles showed little to no sensitivity to rainfall. For this reason, this study used the volumetric water content (VWC) to validate the numerical model. For each software program used, the initial volumetric water content values fell within the 17–30% range as reported in the experiment.

The variation in ${\theta }_w$ with time for the different monitoring points is shown in Figure 20, Figure 21 and Figure 22. The numerical results exhibit similar trends to those observed in the experiment. In Figure 20a,b, the VWC at Section 1 is comparatively higher to other sections due to transverse water flow resulting in accumulation of rainwater around the slope toe. The monitoring points in Sections 1–3, particularly 0.6 m from the slope surface, are more sensitive to precedent rainfall, similar to observations made by Lee et al. (2019) [20]. At 0.9 m, 1.2 m, 1.5 m, and 1.8 m from the surface, Figure 21 and Figure 22b–e, the volumetric water content response to preceding rainfall predicted by the numerical software programs is less pronounced similar to the experiment. This may be due to the longer time required for the wetting front to reach the lower depths of the sensors. The rise in VWC between and after the preceding rainfall phases could be due to the low $k_s$ (=3.0 × 10⁻⁶ m/s) used in numerical analysis which prolonged the time for rainwater infiltration. This phenomenon is most visible in RS2 and to a lesser extent in PLAXIS 2D.

The VWC values from the numerical software programs and the flume experiment during the extreme rainfall stage are in good agreement. The VWC values show a strong response to rainfall infiltration due to increased intensity and rainfall duration. The $q/k_s$ ratio (=6.48) means a lot of the rainwater is discharged as runoff which is a contributing factor to the slope failure. Generally, the SEEP/W seepage model is most sensitive to extreme rainfall compared to other software. At the peak of extreme rainfall duration, SEEP/W VWC values from the selected monitoring points ranged from 37–37.75% ($S$ ~ 100%). At the time of failure, (19 October 1935 h), VWC readings from most monitoring points in the numerical software programs and the experiment ranged from 35–37.75% ($S$ = 90–100%). Tohari et al. (2007) [85] noted that a level of soil saturation of around 85% is sufficient to initiate slope failure. Immediately after the failure time, some experiment readings, as seen in Figure 21b,c and Figure 22b,c,e, decreased sharply because of soil movement resulting in the sensors being exposed to the air.

Factor of Safety $FoS$ vs. Time

The $FoS$ vs. time chart is shown in Figure 23. For this case, the preceding rainfall, which is meant to represent antecedent conditions, has minimal effect on the stability of the slope. During the extreme rainfall (1 and 2), it is observed that the fully coupled PLAXIS numerical model predicts a lower critical $FoS$ compared to RS2 and SLOPE/W. The PLAXIS $FoS$ drops below 1, ($FoS$ = 0.97) around the same time the slope failed in the experiment (19 October 1940 h). RS2 and SLOPE/W show $FoS$ values of 1.17 and 1.348, respectively, at the observed failure time. SLOPE/W based on the LEM assumes that shear strength is mobilized simultaneously along the critical slip surface which may provide a higher $FoS$ value than FE-based software PLAXIS and RS2.

To confirm the $FoS$ results at the failure time (19 October 1940 h), the suction contours from PLAXIS 2D, RS2, and SEEP/W, are manually imported into Optum G2. The negative pore pressures are converted to cohesion values using the distinct unsaturated strength equations shown in Table 13 for slope stability analysis. The $FoS$ results between Optum G2 and the other 3 software programs show close similarity as seen in

Table 15. $FoS$ comparison at observed failure time.

Software	$\mathbf{Factor} \mathbf{of} \mathbf{Safety} \left(\mathit{F}\mathit{o}\mathit{S}\right)$
		Optum G2
SLOPE/W	1.340	1.278 a (4.6%) b
RS2	1.170	1.149 (1.8%)
PLAXIS 2D	0.970	1.069 (−10.2%)

^a $FoS$ is the average of the upper and lower bound, ^b $FoS$% difference between Optum and other software in brackets.

. The failure surfaces of Optum in Figure 24a–c are quite similar to the other three software programs.

It can be deduced that the observed $FoS$ and failure surface differences between SLOPE/W, RS2, and PLAXIS 2D, are due to the differences in pore water pressure resulting from the coupling approach and the failure surface search criteria implemented in each software program.

Slope Failure Mechanism and Deformation

Lee et al. (2019) [20] noted that the slope failure mechanism was a complex combination of progressive shallow failure from the bottom to the top of the slope and a circular slide at the upper part of the slope. The failure surfaces from PLAXIS 2D, RS2 (based on incremental shear strain), and SLOPE/W (based on the LEM) at the recorded time of failure, are presented in Figure 24a–c. Even though the circular slide at the top of the slope is not captured due to possible limitations of the FE programs [86], this phenomenon can be seen in the total displacement contours in Figure 25a,b which show deformation progressing upwards from the slope toe and at the upper part of the slope.

In SLOPE/W, the location of the critical slip surface is dependent on the slip search criteria specified by the user. The critical slip surface in Figure 24c shows the lowest $FoS$. If the slip search criteria in SLOPE/W are localized at the toe or crest area to match the experiment observations, the resultant $FoS$ is unusually high. The practicing engineer can bypass this limitation by exporting the pore pressure profiles in form of apparent cohesion to FE software, such as Optum G2, for slope stability analysis as demonstrated in case studies 1 and 2. This would yield a more realistic failure surface based on a calculation of stresses.

Figure 25c displays the total displacement vs. time results for the slope toe (S1), mid-slope (S2), and slope crest (S4), in PLAXIS 2D and RS2. For the fully coupled PLAXIS 2D model, displacement commences at the slope toe (S1) at the end of extreme rainfall 1, 18 October 2230 h. Hours later, from 19 October 1430 h onwards, there is a sharp rise in displacement that continues up to the recorded failure time, 19 October 1940 h. For the other two monitoring points S2 (mid-slope) and S4 (slope crest), a significant rise in displacement is observed at the end of extreme rainfall 2, 19 October 1655 h up to the failure time, 19 October 1940 h. The monitored PLAXIS 2D displacements demonstrate slope deformation progressing from the toe area to the upper parts of the slope with time, similar to what was observed in the experiment. Displacements recorded in the partially coupled RS2 model are minimal and barely change with time. LE-based programs such as SLOPE/W provide no information on deformation.

Soil Stress Paths

The 2-D stress paths at 1 m from the slope surface (slope toe and crest) are plotted in Figure 26a–d to show the soil stress state at different stages of the flume experiment. For convenience purposes, only six stages are presented for different wetting and drying cycles. Stages 1–2, 2–3, and 4–5 are wetting cycles, while 3–4 represent 5–6 and are drying cycles, the latter being the last cycle prior to slope failure. Although 1–2 is termed as a ‘wetting cycle’, it represents the preceding rainfall stage marked by wet and dry periods. In PLAXIS 2D, Figure 26a,c, the 1–2 ‘wetting cycle’ has minimal effect on the stress path. In RS2, Figure 26b,d, the stress paths at the slope toe and crest slightly shift to the right for the ‘wetting cycle’ 1–2. This may be due to a decrease in water content during the drying period before the start of extreme rainfall.

In both software, wetting cycles 2–3 and 4–5 shift the stress paths toward the left due to loss of suction as rainfall infiltrates the soil, leading to a decrease in effective stress and shear strength. In PLAXIS 2D, the stress paths 2–3 and 4–5 moved left and touched the failure envelope, $k_f$ line, indicating critical conditions. The stress path 5–6 in Figure 26a slightly shifts upper right since it is a drying cycle and there is an accumulation of soil moving downward towards the toe and crossing the $k_f$ line, representing failure. At the slope crest, Figure 26c, the stress path moves lower left in the q-p′ space. The loss in deviatoric stress, q, may be due to a decrease in overburden as a result of the circular slide observed at the top of the slope.

The corresponding stress path 5–6 predicted by RS2, Figure 26b,d, indicate critical conditions but are quite different to PLAXIS 2D. Despite demonstrating the effect of wetting cycles 2–3 and 4–5 on the soil stress state at the slope toe, the stress paths 5–6 in RS2, Figure 26b,d, seem to not match experimental observations, particularly with regards to slope deformation. Similarly, at the slope crest, Figure 26d, the stress paths 3–4 and 4–5 shift vertically upwards and downwards, respectively. There is no change in effective mean stress, p′.

For the landslide flume case study, the fully coupled PLAXIS numerical model shows a comparatively more realistic representation of the recorded slope behavior in terms of VWC variation with time, failure time, deformation, and soil stress path at different monitor points. Ng et al. (2020) [87] discussed the state-dependency of unsaturated soil behavior linking several variables such as incremental volumetric strain, deviator strain, degree of saturation to incremental mean net stress, deviator stress, and suction. PLAXIS 2D follows Boit’s consolidation theory [88], enabling the simultaneous calculation of deformation and seepage with time-dependent boundary conditions in unsaturated soils. The elastoplastic behavior of the soil skeleton and suction dependency of the degree of saturation and relative permeability is fully considered. To describe the mechanical behavior of unsaturated soil in PLAXIS, Bishop’s stress and suction are used as state variables [89,90]. The fully coupled model also analyses the position of the phreatic line with the time and adjusts the soil weight accordingly [19].

Notably, the effect of hydraulic hysteresis due to the wetting and drying cycles is not considered. Additionally, the Mohr–Coulomb (M–C) model was used, which assumes constant stiffness for the soil layer, hence the stress dependency of soil stiffness is not accounted for. The M–C model is suitable for the first estimation of soil deformation. For more advanced analysis, the hardening soil model which considers the stress dependency of soil stiffness can be used as documented in previous coupled rainfall infiltration studies [23,24].

For this case, the SEEP/W and SLOPE/W computer programs based on the uncoupled approach could be used for preliminary analysis in order to understand the flow regime and stability response of the slope with time. This may be followed by partially coupled (RS2) and fully coupled seepage analysis (PLAXIS 2D), the latter of which explicitly considers the interdependence of pore pressures and deformation.

5. Conclusions

This study conducts a parametric study on unsaturated silty slopes under a 1-D steady flow rate. It is carried out to identify when it is appropriate to consider apparent cohesion in slope stability analysis. A comparative review of widely used computer programs in rainfall seepage unsaturated slope problems is presented. The following findings are made:

• Part I

(a)

For unsaturated silty slopes under 1-D, steady state infiltration, the main controlling parameters determining the contribution of apparent cohesion include the water table location, $D_{sw}$ the infiltration rate, $\raisebox{1ex}{$q$}\!\left/ \!\raisebox{-1ex}{$k_s$}\right.$, pore size parameters, α and n, and slope height, H. The contribution of suction and apparent cohesion may be considered for cases with a water table located below the mid-height of the slope ($\raisebox{1ex}{$H_w$}\!\left/ \!\raisebox{-1ex}{$H$}\right.<$ 0.5). For slopes with a high water table ($\raisebox{1ex}{$H_w$}\!\left/ \!\raisebox{-1ex}{$H$}\right.\ge$ 0.5), the contribution of suction to slope stability may be ignored. The slope height, H, is a primary controlling factor for the $FoS$ of slopes with H ≥ 20 m and a fairly high water table.

(b)

In borderline silt/sand slopes under moderate infiltration rate, $\raisebox{1ex}{$q$}\!\left/ \!\raisebox{-1ex}{$k_s$}\right.$ = −0.3–(−0.6), apparent cohesion is considerable due to increased water contact between soil particles. Silt/clay slopes maintain sizeable suction values for a wider range of infiltration rates, $\raisebox{1ex}{$q$}\!\left/ \!\raisebox{-1ex}{$k_s$}\right.$ = 0–(−0.7) and slope height, H = 10–30 m.

• Part II

(c)

In the first case study, commercial software programs RS2, Slide 2D, PLAXIS LE, SEEP/W, SLOPE/W, and PLAXIS 2D are used to run seepage and stability analysis of an unsaturated silty slope subjected to rainfall infiltration. The output results are quite similar for most software programs with notable differences mainly due to distinct slip surface search methods, failure mode, and coupling approach in each computer program. The fully coupled model in PLAXIS 2D yields a relatively lower critical $FoS$ (<1.0), sharper wetting front, and shallower failure surface, compared to other software programs after 6 h of rainfall. Distinct failure modes in RS2 and Slide 2D cause differences in the $FoS$ at the initial and end stages despite having the same FE seepage engine. These findings can aid the geotechnical engineer’s judgment on solving similar slope-related problems.

(d)

In the second case study, the computer programs SEEP/W, RS2, and PLAXIS 2D, are used for slope stability analyses. The volumetric water content vs. time, failure time, and failure mechanism are used to validate the obtained results. The displacement vs. time and soil stress paths for PLAXIS 2D and RS2 are plotted and compared. Notable differences in the $FoS$ vs. time and failure surface at the recorded failure time can be attributed to pore water pressure differences arising from the coupling approach and slip surface search criteria in each computer program. Despite the limitation in FE programs to capture progressive failure, the total displacement contours, displacement vs. time, and soil stress paths in PLAXIS 2D, provide a resemblance of the slope movement observed in the experiment. The fully coupled PLAXIS 2D model results are comparatively closer to the experimental results. The partially coupled RS2 model does not sufficiently capture the slope deformation, particularly in terms of total displacement vs. time and stress paths.

(e)

For the partially coupled and fully coupled numerical models in RS2 and PLAXIS 2D, respectively, the $k_s$ is identified as the main controlling parameter of the hydro-mechanical response and failure time of the unsaturated slope under rainfall infiltration. The SEEP/W model is less sensitive to $k_s$. The strength parameter c′ only affects the calculated $FoS$ values. Practically speaking, it is difficult to accurately measure $k_s$. For good geotechnical practice, slopes with an $FoS$ close to 1.0 can be considered to be unsafe and about to collapse.

(f)

Despite the technological advancements in LE-based slope stability programs, the location of the critical slip surface is still dependent on the slip surface search criteria specified by the user. Additionally, the LEM assumes the shear strength is mobilized simultaneously along the entire slip surface which is not realistic. Engineers may bypass such limitations by importing pore pressure contours in form of apparent cohesion to FE software which locates the failure surface based on a calculation of stresses, as demonstrated in case studies 1 and 2.