Slope Stability Analysis of Unsaturated Soil Slopes Based on the Site-Specific Characteristics: A Case Study of Hwangryeong Mountain, Busan, Korea

Abstract

Shallow slope failures occur almost every year during the rainy season. Continuous observation of the meteorological parameters and hydrological characteristics is required to more clearly understand the triggering mechanisms of shallow slope failure. In addition, influential factors, such as type of relative permeability models, air flow, and variation of hydraulic conductivity associated with stress–strain behavior of soil, have significant effects on the actual mechanism of rainfall infiltration. Real-time data including hourly rainfall and pore water pressure in response to rainfall was recorded by devices; then, the change in pore pressure from the devices was compared to the results from the infiltration analysis with applications of three relative permeability models, air flow, and the coupled hydro-mechanical analysis to examine an appropriate site-specific approach to a rainfall infiltration analysis. The infiltration and stability analyses based on the site-specific hydrologic characteristics were utilized to create maps of safety factors that depend on the cumulative rainfall. In regions vulnerable to landslides, rainfall forecast information and safety factor maps built by applying various rainfall scenarios can be useful in preparing countermeasures against disasters during the rainy season.

1. Introduction

Slope failures are caused by a combination of factors, such as geological, topographical, and hydrological characteristics, slope surface cover, and climate. The slope failures triggered by rainfall have frequently been reported worldwide [1]. In Korea, mountainous areas occupy approximately 70 percent of the territory, and 70 percent of the annual rainfall higher than 1000 mm is concentrated during a rainy season that occurs between June and September [2]. Numerous slope failures occur almost every year during the rainy season due to localized heavy rains, resulting in loss of life and property damage [3]. In Korea, slope failures, caused by arrival at the critical wetting band depth because of rainfall infiltration, show a tendency to occur on the boundary between a shallow weathered residual soil layer and an underlying rock layer (weathered rock or bedrock) [4]. When water begins to infiltrate the unsaturated soil, matric suction (i.e., negative pore water pressure) starts to decrease, and then loss of matric suction causes a slope failure because of the decrease in the shear strength of the soil [5].

According to Bordoni et al. [6], a physically-based hydrological triggering mechanism concerning the site-specific characteristics of the involved soils can be considered to assess rainfall-induced shallow landslides. Additionally, continuous monitoring of the meteorological parameters and the hydrological characteristics of the unsaturated soils is required to understand the triggering mechanisms of shallow landslides. Therefore, it is necessary to consider data based on real time observations, such as rainfall information and the response of pore water pressure to rainfall, and site-specific characteristics, including physical and hydrological characteristics of soils, to identify infiltration processes adequately.

Not only intensity and duration of rainfall but also hydraulic properties of soils near the surface, that is permeability, unsaturated hydraulic conductivity, and soil-water characteristics, have a great effect on the rainfall infiltration [7,8]. In particular, the unsaturated hydraulic conductivity is an important factor in understanding the hydraulic characteristics of unsaturated soils [7]. The unsaturated hydraulic conductivity function has been estimated by a variety of methods through the use of the soil-water retention curve (SWRC) or the simplified model because of the difficulty of experimental measurements [9,10,11,12,13,14]. When rainfall infiltration characteristics of a specific site are to be understood, it is necessary to evaluate the influence of the unsaturated hydraulic conductivity and to identify an appropriate relative permeability model, which determines the unsaturated hydraulic conductivity for the specific site.

Because the ground consists of three-phase material, including soil particles, water, and air, to strictly interpret rainfall infiltration, it is necessary to consider a fully coupled formulation of water and air flow and the stress–strain behavior of the soil [15]. However, assumptions for simplification have generally been used to solve the differential equation of Richards [16]: (1) the single-phase flow of only water is applied by ignoring air flow and the stress–strain behavior of the soil [17]; (2) infiltration–strain behavior is considered, but the air pressure in the ground is assumed to be equal to the atmospheric pressure, based on the assumption that air flow in the ground is relatively free compared to water flow [18,19,20]; (3) the two-phase flow of water and air is applied, but the stress–strain behavior of the soil is ignored [21,22]; and (4) hysteretic wetting and drying states are not considered, affecting variably saturated hydrologic response and slope failure [20,23,24,25].

The interaction of water and air on the contact surface causes air flow and air compression during rainfall infiltration through voids in unsaturated soils, and it is well known that air flow influences the infiltration of water in unsaturated soils, based on experiments and interpretations in relation to water infiltration [21,26]. In addition, porosity varies depending on changes in the degree of saturation and the stress due to externally applied loads and changes in soil moisture and suction during infiltration-induced deformation [27,28]. Changes in microstructure and permeability are caused by variations in the passage for water and air flow [29]. Hydraulic hysteresis is significant for applying the soil–water–air coupling strictly [19,30,31,32], affecting variation in water content and slope stability [33].

This paper proposes a procedure which facilitates identifying characteristics of hydrologic responses to rainfall adequately in specific sites and constructing an optimal combination of models utilized for rainfall infiltration analysis based on the site-specific characteristics. In this study, based on the soil properties obtained from field/laboratory tests and the measured rainfall data, infiltration analyses were performed to interpret hydrologic responses that depended on applications of different relative permeability models, two-phase flow of water and air, and the fully coupled formulation of the stress–strain behavior of the soil. The accurate site-specific approach was examined that well-modeled a process of infiltration for the study site comparing the measured and simulated pore water pressures. Finally, infiltration and slope stability analyses based on the site-specific hydrologic characteristics were performed in response to the assumed rainfall scenario to generate the maps of the safety factor that depended on cumulative rainfall, which could be useful for predicting and preventing landslide damage and for preparing countermeasures against disasters.

2. Study Area and Observation Devices

Hwangryeong Mountain is located in the center of Busan Metropolitan City and is in a typical downtown area, as shown in Figure 1. It extends from 35°8′00″–37°10′30″ at the north and from 129°4′00″–129°7′00″ at the east. The area of Hwangryeong Mountain (13,020,750 m²) is completely encircled by buildings and roads. The altitude at the top of the mountain is 427 m above sea level. The geology of the Hwangryeong Mountain area consists mainly of sedimentary rocks from the Mesozoic to Cretaceous period and consists of volcanic rocks, such as andesite and volcanic breccia. In the mountain area, various facilities, such as a youth training center, a broadcasting station, a training center in the subordinate of Busan City Hall, and an industrial road are located.

At 11:50 a.m. on 11 September 1999, approximately 140,000 m³ of soil erosion occurred because of slope failures caused by heavy rain from 9 to 11 September 1999. The three-day accumulation was 200 mm, and the maximum hourly rainfall on the day of the accident was 39 mm/h. The slope failures caused one death, injured three, and resulted in approximately 9.9 million US dollars of restoration expenditure. According to the Ministry of Science, ICT and Future Planning [34], the area of Hwangryeong Mountain was classified as an area vulnerable to landslides.

A rainfall gauge and two pore pressure gauges were installed at depths of 0.5 m and 1.0 m from the ground surface, as shown in Figure 1b. The pore pressure gauge (MPS–6 manufactured by Decagon Devices, Inc.) is able to measure negative pore water pressure up to −1500 kPa. The rainfall gauge (WDR–205 manufactured by Wedaen Co.) is able to measure the real-time rainfall and rainfall intensity (mm/h). This is a tipping bucket-type rainfall gauge, in which rainfall is collected in a 20 cm diameter water gauge. The hourly and daily rainfall and the pore water pressure data per hour were measured in real time. They were sequentially transferred to a data logger (CR1000 manufactured by Campbell Scientific, Inc.), which transmits the data to the monitoring server by using a wireless communication method. In this study, rainfall conditions for infiltration analysis were determined based on the observed rainfall data, and the observed pore water pressure was used for comparison with the results obtained from the infiltration analyses to validate models used in the analyses.

3. Materials and Methods

3.1. Rainfall Infiltration Analysis

3.1.1. Unsaturated Flow

In this study, Fast Lagrangian Analysis of Continua (FLAC) Ver. 7.0 [35], which is finite difference code software available for an application of the two-phase flow (tp-flow) of water and air in a 2D package, was used to conduct the rainfall infiltration analysis in unsaturated soil slopes. The tp-flow option within the FLAC program can be applied to model two immiscible fluids through a porous medium. Assuming that the fluids within the tp-flow option are air and water allows numerical modeling of unsaturated soil conditions [36]. The FLAC program contains an in-built programming language called FISH that allows the user to write and apply functions within the code.

The pore spaces within the porous media are assumed to be completely filled by the two fluids, as expressed by Equation (1). One of the fluids is called the wetting fluid and wets the porous media more than another, which is known as the non-wetting fluid. The subscripts ($w$) and ($g$) represent the wetting fluid and the non-wetting fluid, respectively. In unsaturated soils, the wetting fluid and the non-wetting fluid can be replaced by water and air, respectively.

$$S_w+S_g=1,$$

(1)

where $S_w$ is the saturation of the wetting fluid, and $S_g$ is the saturation of the non-wetting fluid.

The difference between the pressures of the two fluids is the capillary pressure, as expressed by Equation (2). Equation (3), proposed by van Genuchten [14], illustrates a unique relationship between the capillary pressure and the effective saturation related to both the actual saturation and the residual saturation.

$$P_C=P_g-P_w$$

(2)

$$S_e={\left[{\left(P_C/P_0\right)}^{\frac{1}{1-a}}+1\right]}^{-a}=\frac{S_w-S_w^r}{1-S_w^r}=\frac{\theta -{\theta }_r}{{\theta }_s-{\theta }_r}$$

(3)

where $P_C$ is the capillary pressure, $P_g$ is the pressure of the non-wetting fluid, $P_w$ is the pressure of the wetting fluid, $S_e$ is the effective saturation, $P_0$ is the parameter that depends approximately on the air-entry value, $a$ is the parameter equivalent to the van Genuchten m constant, $S_w$ is the actual saturation, $S_w^r$ is the residual saturation, $\theta$ is the actual volumetric water content, and ${\theta }_r$ and ${\theta }_s$ are the residual volumetric water content and saturated volumetric water content, respectively.

The matric suction ($\psi$) is derived by using Equation (3) as follows:

$$\psi =P_C=P_0{\left[{\left(S_e\right)}^{-\frac{1}{a}}-1\right]}^{1-a}.$$

(4)

In the tp-flow option, transport of the wetting fluid and the non-wetting fluid is calculated using Darcy’s law. Equations (5) and (6) exhibit the velocities of the wetting fluid and the non-wetting fluid, respectively.

$$q_i^w=-k_{ij}^w{k}_r^w\frac{\partial }{\partial x_j}\left(P_w-{\rho }_w{g}_k{x}_k\right),$$

(5)

$$q_i^g=-k_{ij}^w\frac{{\mu }_w}{{\mu }_g}{k}_r^g\frac{\partial }{\partial x_j}\left(P_g-{\rho }_g{g}_k{x}_k\right),$$

(6)

where $k_{ij}^w$ is the saturated mobility coefficient, $k_r$ is the relative permeability which is a function of the saturation of the wetting fluid, $\frac{{\mu }_w}{{\mu }_g}$ is the viscosity ratio between the two fluids, ${\rho }_w$ and ${\rho }_g$ are the densities of the wetting and non-wetting fluids, respectively, and $g$ is the gravitational acceleration.

The mobility coefficient used in FLAC represents the permeability and can be described using the hydraulic conductivity, commonly used when Darcy’s law is expressed in terms of head, as expressed by Equation (7).

$$k=k_s/(g{\rho }_w),$$

(7)

The relative permeability of the wetting fluid, which is described as a function of saturation, represents a percentage of maximum permeability of the wetting fluid and is used as one component to calculate the unsaturated hydraulic conductivity. The relative permeability used in FLAC was proposed by van Genuchten [14] by integrating the soil-water retention curve (SWRC) and using a closed form equation based on the Mualem [13] capillary model, as expressed by Equation (8). The relative permeability of the non-wetting fluid used in FLAC was proposed by Lenhard and Parker [37] and it is expressed as Equation (9).

$$k_r^w={{(S}_e)}^b{\left\{1-{\left[1-{\left(S_e\right)}^{\frac{1}{a}}\right]}^a\right\}}^2,$$

(8)

$$k_r^g={\left(1-S_e\right)}^c{\left[1-{{(S}_e)}^{\frac{1}{a}}\right]}^{2a},$$

(9)

where $b$ and $c$ are the constants for water and air, respectively, which are both equal to 0.5 [36].

The van Genuchten–Mualem (vG–M) relative permeability model was widely applied to the infiltration analysis in previous research [15,21,24,27,36,38] because it can be easily derived directly using the constant of the van Genuchten SWRC model. However, numerical instability caused by a very steep gradient of the vG–M model in a range of small matric suction values can induce a highly nonlinear solution [39]. To supplement the problem, the original fit of the van Genuchten SWRC can be modified to a non-smooth curve at small matric suction before integrating the SWRC [7,39,40,41]. Oh et al. [7] proposed a modified relative permeability model that was shown to accurately predict the unsaturated behavior for Korean weathered soils using the van Genuchten SWRC modified to a non-smooth curve at a range of small matric suction values. Figure 2 exhibits a conceptual diagram of the modified van Genuchten SWRC suggested by Oh et al. [7]. An arbitrary parameter to replace $P_0$, $P_0^{\prime }$ (=0.02$P_0$), was introduced to modify the SWRC through tangentially extrapolating from ($P_0^{\prime }$, ${S_{e0}}^{\prime }$) to full saturation on the log–P_c scale. The modified SWRC can be expressed as Equations (10)–(12). Substituting these three equations into a closed form equation based on Mualem [13], described as Equation (13), the modified relative permeability is finally derived. In this study, the vG–M model and the modified van Genuchten–Mualem (mvG–M) model were applied to understand the effects on rainfall infiltration depending on the types of relative permeability. In addition, the relative permeability used in the single-phase flow (sp-flow) option in FLAC, expressed as Equation (14), with more gentle gradients in a range of small matric suction values was considered. Because the mvG–M model and the relative permeability used in the sp-flow option are basically not supported in the tp-flow option in FLAC, they were programmed to be applied at every time step using FISH, the in-built programming language.

$$P_C=P_0{\left[{\left(S_e\right)}^{-1/a}-1\right]}^{1-a},\mathrm{if}{S}_e<{S_{e0}}^{\prime },$$

(10)

$$P_C=0.02P_0\exp \left\{\lambda \left(S_e-S_{e0}^{\prime }\right)\right\},\mathrm{if}{S}_{e0}^{\prime }\le S_e<1,$$

(11)

$$0\le P_C<P_s,\mathrm{if}{S}_e=1,$$

(12)

$$k_r^w=S_e^b\frac{{\left({{\displaystyle \int }}_0^{S_{e0}{}^{\prime }}\frac{d{S}_e}{\psi }+{{\displaystyle \int }}_{S_{e0}^{\prime }}^{S_e}\frac{d{S}_e}{\psi }\right)}^2}{{\left({{\displaystyle \int }}_0^{S_{e0}{}^{\prime }}\frac{d{S}_e}{\psi }+{{\displaystyle \int }}_{S_{e0}{}^{\prime }}^1\frac{d{S}_e}{\psi }\right)}^2},$$

(13)

$$k_r^w={\left(S_e\right)}^2(3-2S_e),$$

(14)

where $\mathsf{\lambda }=\frac{{\left\{1+{\left(0.02\right)}^{1/(1-a)}\right\}}^{a+1}}{-a{\left(1-a\right)}^{-1}{\left(0.02\right)}^{1/(1-a)}}$, and $P_s=0.02P_0\exp \left\{\lambda \left(1-S_{e0}^{\prime }\right)\right\}$.

Balance equations of the two fluids are described as Equations (15) and (16), that were derived combining the fluid balance laws with the fluid constitutive laws for the wetting fluid and the non-wetting fluid.

$$n\left[\frac{S_w}{K_w}\frac{\partial P_w}{\partial t}+\frac{\partial S_w}{\partial t}\right]=-\left[\frac{\partial q_i^w}{\partial x_i}+S_w\frac{\partial \epsilon }{\partial t}\right],$$

(15)

$$n\left[\frac{S_g}{K_g}\frac{\partial P_g}{\partial t}+\frac{\partial S_g}{\partial t}\right]=-\left[\frac{\partial q_i^g}{\partial x_i}+S_g\frac{\partial \epsilon }{\partial t}\right],$$

(16)

where $n$ is the porosity, $K$ is the fluid bulk modulus, $t$ is the time, and $\epsilon$ is the volumetric strain.

3.1.2. Balance of Momentum

The balance of momentum is described as Equation (17).

$$\frac{\partial {\sigma }_{ij}}{\partial x_j}+\rho g_i=\rho \frac{d\dot{u_i}}{dt},$$

(17)

where $\sigma$ is the total stress, $\rho$ is the bulk density, and $\dot{u}$ is the velocity. In unsaturated soil conditions, the bulk density is described as Equation (18).

$$\rho ={\rho }_d+n(S_w{\rho }_w+S_g{\rho }_g),$$

(18)

where ${\rho }_d$ is the dry density of soil, ${\rho }_w$ is the density of water, and ${\rho }_g$ is the density of air.

Bishop’s effective stress [42] was used to analyze the failure of unsaturated soil in the Mohr–Coulomb constitutive model. The yield criterion of the Mohr–Coulomb constitutive model is described as Equation (19). Bishop’s effective stress is defined as Equation (20).

$${\tau }_{max}={{\sigma }^{\prime }}_b\tan {\varphi }^{\prime }+c^{\prime },$$

(19)

$${{\sigma }^{\prime }}_b=\left(\sigma -P_g\right)+\chi (P_g-P_w),$$

(20)

where ${\tau }_{max}$ is the shear strength, ${{\sigma }^{\prime }}_b$ is Bishop’s effective stress, ${\varphi }^{\prime }$ is the effective friction angle, $c^{\prime }$ is the effective cohesion, andis the effective cohesion, and $\chi$ is the matric suction coefficient, which varies between 0 (dry condition) and 1 (fully saturated condition).

The matric suction coefficient can be approximately replaced by the degree of saturation [24,38,43,44,45]. By replacing $\chi$ with the degree of saturation, Equation (21) was derived from Equation (20). Then, Equation (19) may be expressed as Equation (22).

$${{\sigma }^{\prime }}_b=\sigma -(S_w{P}_w+S_g{P}_g),$$

(21)

$${\tau }_{max}=\left[\sigma -\left(S_w{P}_w+S_g{P}_g\right)\right]\tan {\varphi }^{\prime }+c^{\prime },$$

(22)

where $S_w{P}_w+S_g{P}_g$ is the pore fluid pressure.

3.1.3. Coupled Hydro-Mechanical Analysis

Porosity varies with the changing degree of saturation and stress due to externally applied loads and changes in soil moisture and suction during infiltration-induced deformation [27,28]. Changes in microstructure and permeability are caused by variations in void space [29]. Hydraulic conductivity depending on changes in porosity can be derived from the Kozeny–Carman equation [46] (e.g., [47,48,49]) and is expressed as Equation (23). Because changes of porosity and hydraulic conductivity based on the slope deformation during rainfall infiltration cannot be considered in FLAC, the equation was programmed to be applied at every time step using FISH. An application of changes in porosity and hydraulic conductivity greatly affects the time required for computation because it affects the critical time interval and the matrix of elements.

$$k_s=k_{s0}{\left(\frac{n}{n_0}\right)}^3{\left(\frac{1-n_0}{1-n}\right)}^2,$$

(23)

where $k_{s0}$ is the initial hydraulic conductivity, and $n_0$ is the initial porosity.

3.2. Slope Stability Analysis

Shallow slope failures caused by rainfall can be interpreted by the infinite slope analysis method because the failure surface is usually parallel to the slope surface [4]. In this study, the limit equilibrium method was applied to the infinite slope to calculate the safety factor considering that the surface of shallow slope failure tends to be parallel to the slope surface in Korea [50]. Figure 3 exhibits a conceptual diagram of an infinite slope with a shallow impermeable layer. The safety factor of the infinite slope is described as Equation (24) using the Mohr–Coulomb failure criterion (Equation (22)) with the shear strength of soil and Bishop’s effective stress (Equation (21)).

$$F_s=\frac{{\tau }_f}{{\tau }_m}=\frac{{\sigma }_b^{\prime }\tan {\varphi }^{\prime }+c^{\prime }}{W\sin \beta \cos \beta }=\frac{\left[W{\cos }^2\beta -\left(S_w{P}_w+S_g{P}_g\right)\right]\tan {\varphi }^{\prime }+c^{\prime }}{W\sin \beta \cos \beta },$$

(24)

where ${\tau }_f$ is the shear strength of soil, ${\tau }_m$ is the shear stress on failure surface, $\beta$ is the slope angle, and $W$ is the weight of a soil slice.

The total unit weight of the soil varies according to changes in saturation because of rainfall infiltration; therefore, the weight of a soil slice per unit area on the potential failure surface can be expressed as Equation (25) to account for an increase in the unit weight of the soil during rainfall.

$$W=\underset{0}{\overset{z_w}{{\displaystyle \int }}}{\gamma }_{t}dz=\underset{0}{\overset{z_w}{{\displaystyle \int }}}\left[{\gamma }_d+n{S}_w\left(z\right){\gamma }_w\right]dz,$$

(25)

where $z_w$ is the depth to the potential failure surface, ${\gamma }_t$ is the total unit weight of the soil, ${\gamma }_d$ is the dry unit weight of the soil, ${\gamma }_w$ is the unit weight of water, and $S_w\left(z\right)$ is the saturation based on depth.

3.3. Material Properties

Material properties of the soil, water, and air are required for modeling a slope to conduct infiltration analysis. Undisturbed soil samples were taken at a depth of approximately 1 m near the location of the pore pressure gauge shown in Figure 1b, and field/laboratory tests, such as a field permeability test, a direct shear test, and a soil-water characteristic curve (SWCC) test, were conducted to obtain material properties of the soil. Material properties of water and air were acquired by referring to Cho [24].

Table 1. Material properties of soil, water, and air.

Parameter		Value
Dry density of soil, ${\rho }_d$		1250 kg/m3
Initial porosity, no		0.46
Soil classification (USCS)		SP–SM
Cohesion, $c^{\prime }$		0.1 × 103 Pa
Friction angle, ${\varphi }^{\prime }$		33°
Dilation angle, $\Psi$		0°
Bulk modulus, K		33.3 × 106 Pa
Shear modulus, G		11.1 × 106 Pa
Hydraulic conductivity, ks		1.97 × 10−6 m/s
Soil-water characteristics (wetting)	Residual saturation, Sr	0
	P0	3000 Pa
	van Genuchten parameter, a	0.52
Constant in relative permeability for water, b		0.5
Constant in relative permeability for air, c		0.5
Viscosity ratio, ${\mu }_w/{\mu }_a$		56
Water density, ${\rho }_w$		1000 kg/m3
Air density, ${\rho }_a$		1.25 kg/m3
Bulk modulus of water, Kw		2 × 109 Pa
Bulk modulus of air, Ka		1 × 105 Pa

exhibits material properties of soil (dry density, initial porosity, cohesion, friction angle, bulk modulus, shear modulus, hydraulic conductivity, soil-water characteristics for van Genuchten [14] SWRC), water, and air that were applied in this study. Kim et al. [51] conducted laboratory tests using undisturbed soil samples that were taken from 83 sites around the Hwangryeong Mountain. They reported ranges of the dry density, cohesion, friction angle, and hydraulic conductivity obtained from the tests (e.g., 1100–1500 kg/m³ with an average of 1220 kg/m³, 0.1–6.6 kPa with an average of 2 kPa, 32–39° with an average of 36°, and $1.9\times {10}^{-6}$–$5.2\times {10}^{-4}$ m/s with an average of $5.8\times {10}^{-5}$ m/s, respectively). The dry density (1250 kg/m³), cohesion (0.1 kPa), friction angle (33°), and hydraulic conductivity ($1.97\times {10}^{-6}$ m/s) applied in this study are in the ranges reported by Kim et al. [51]. Therefore, it was assumed that material properties of the soil samples were reasonable to represent characteristics of soils adjacent to the pore pressure gauge.

Figure 4 exhibits hydraulic properties used for infiltration analyses in this study. Experimental tests were performed through applying the pressure plate method [52] for measuring volumetric water content that depends on matric suction under both wetting and drying conditions to consider the water retention curve hysteresis. Figure 4a shows wetting and drying SWCCs fitted with the van Genuchten [14] SWRC model (i.e., coefficients of the model are given in Table 1). Using such coefficients, the vG–M relative permeability curves [14] for water and air, the mvG–M relative permeability curve [7], and the relative permeability used in the sp-flow option in FLAC (K_r–sp–FLAC) were computed, as shown in Figure 4b. Effects of the three relative permeability curves for water on the infiltration rate were examined in this study, whereas the vG–M relative permeability curve for air was consistently applied. The vG–M relative permeability curve for water was steepest near saturation among the three relationships, and it has the smallest values in the entire range of effective saturation. On the other hand, the K_r–sp–FLAC has the largest values.

3.4. Slope Geometry and Boundary and Initial Conditions

In this study, the infiltration analysis was performed on the unsaturated soil slope with a depth of 3 m on the layer of weathered rocks, as shown in Figure 5, considering the depths from ground surfaces to weathered rocks observed during the field investigation. Points (a) and (b) (section C–H) at vertical depths of 0.5 m and 1.0 m from the surface, respectively, indicate observation points to check changes in pore water pressure. The slope was determined to be approximately 34° considering the actual slope angle at the location where the pore pressure gauges were installed. Material properties of the soil slope applied in the analysis are listed in Table 1. Material properties of the layer of weathered rocks (e.g., dry density of 2750 kg/m³, cohesion of 43.4 kPa, and friction angle of 43.1°) were applied referring to La et al. [53] who conducted in-situ investigations and laboratory tests using rock samples taken from the area of Hwangryeong Mountain. It was assumed that soils and weathered rocks of the slope are homogeneous, and material properties are isotropic.

Water flux was applied at the ground surface (ABCDE) as a boundary condition to simulate infiltration of rainfall. A displacement into the parallel direction was fixed at the left and right sides (AK and EL), and the bottom surface (KL) was fully constrained. The initial air pressure on the slope was assumed to be atmospheric pressure. Additionally, the slope surface (ABCDE) was subjected to the atmospheric boundary condition and other boundaries were assumed to be airtight. For the sp-flow option, the discretization of the domain, initial conditions, and boundary conditions were similar to those for the tp-flow option, except that air pressure was set to zero.

4. Results and Discussion

4.1. Site-Specific Rainfall Infiltration Analysis

Rainfall infiltration was simulated based on the slope geometry as shown in Figure 5, material properties summarized in Table 1, and rainfall data measured from 20:00 on 27 September 2016 to 23:00 on 28 September 2016, as shown in Figure 6. It was assumed that initial pore water pressures were consistently −32.95 kPa at the soil layer (steady-state condition), which was the average value of initial pore water pressures at the beginning of the rainfall measured at depths of 0.5 and 1.0 m from the ground surface (–35.2 and –30.7 kPa, respectively). To investigate an accurate approach to the site-specific infiltration analysis, the influence of two types of fluid flow models (sp-flow model and coupled hydro-mechanical (CHM) model considering tp-flow of water and air) and three relative permeability models on hydrologic responses was examined through comparing changes in pore water pressure at depths of 0.5 and 1.0 m both from observation devices and from simulations. It is necessary to observe pore water pressures at various depths to precisely investigate site-specific infiltration characteristics, whereas several previous studies have observed pore water pressure using few (3–4) gauges at vertical profiles [6,25,54,55,56]. Considering that the soil layer of the slope applied in this study was much shallower than that applied in the previous studies, it was assumed that hydrological behaviors of the soil layer could be evaluated based on the pore water pressure observed at depths of 0.5 and 1.0 m.

Figure 7 shows changes in pore water pressure at two observation points, (a) and (b) shown in Figure 5, obtained from the sp-flow and CHM models with three relative permeability models (e.g., vG–M model, mvG–M model, and K_r–sp–FLAC model). The pore water pressures were comparable with the field observations when the mvG–M model and the K_r–sp–FLAC model were applied. However, increases in the pore water pressure at both depths of 0.5 and 1.0 m were consistently slowest when the vG–M model was applied because the unsaturated hydraulic conductivity from the vG–M model was smallest, as shown in Figure 4b. It can be observed from Figure 7a,c that the increase rate of pore water pressure from the K_r–sp–FLAC model was faster than those from the mvG–M model at a depth of 0.5 m during rainfall with the intensity greater than 0.7 mm/h (i.e., 0–5 h and 16–24 h). This is because the unsaturated hydraulic conductivity from the mvG–M model was smaller than that from the K_r–sp–FLAC model. During weak rainfall (i.e., 6–15 h), however, the increase rate of pore water pressure from the K_r–sp–FLAC model became slower than that from the mvG–M model. The effective saturation gradually increased during rainfall after the infiltrated rainfall reached a depth of 0.5 m, and then the relative permeability also increased significantly when the K_r–sp–FLAC model was applied. As the rainfall intensity became small, the inflow rate of rainfall into the point at a depth of 0.5 m decreased, whereas the water passed through the point more quickly because of large relative permeability. Therefore, the pore water pressure increased slowly during the rainfall with low intensity when the K_r–sp–FLAC model was applied. At a depth of 1.0 m, relative permeability increased quite slowly because rainfall arrived late due to the small saturated hydraulic conductivity ($1.97\times {10}^{-6}$ m/s), and the increase rate of pore water pressure depended on the relative permeability model without effects of rainfall intensity. Therefore, the pore water pressure from the K_r–sp–FLAC model increased faster than that from the mvG–M model at a depth of 1.0 m, as shown in Figure 7b,d.

The pore water pressures simulated from the sp-flow model, as shown in Figure 7a,b, were slightly greater than those simulated from the CHM model considering tp-flow of water and air, as shown in Figure 7c,d. These results could be attributed to effects of air flow and air pressure in the slope, interrupting water flow and propagation of the wetting front, and those are in agreement with the findings by Hu et al. [29] and Cho [24]. In addition, the decrease in hydraulic conductivity caused by the reduction of porosity during the deformation of the slopes could result in low increase rates of the pore water pressure.

Residuals between the measured and estimated pore water pressures ($P{w}_{\mathrm{mea}}$ and $P{w}_{\mathrm{est}}$, respectively) were computed by using an equation, $\mathrm{Residual}=\ln \left(P{w}_{\mathrm{mea}}/P{w}_{\mathrm{est}}\right)$, to evaluate the prediction accuracy of the models. Figure 8 shows the distribution of residuals at depths of 0.5 and 1.0 m against time variation that depended on the types of fluid flow model and relative permeability model. Regardless of the type of fluid flow model applied, the residuals from the mvG–M model and the K_r–sp–FLAC model were consistently close to zero at both depths from the start to the end with the exception that the residuals from the K_r–sp–FLAC model increased rapidly from 20 h. Whereas the residuals from the vG–M model remained close to zero from the start to 18 h and then decreased sharply at depths of 0.5 and 1.0 m from 18 and 20 h, respectively. This indicates that the relative permeability model had a greater effect on the pore pressure change than the fluid flow model during the infiltration analysis. Values for mean and standard deviation of the residuals are summarized in

Table 2. Values for mean and standard deviation of residuals that depended on the two fluid flow models (sp-flow model and CHM model) and three relative permeability models (vG–M model, mvG–M model, and K_r–sp–FLAC model) at depths of 0.5 and 1.0 m.

Fluid Flow Model	Relative Permeability Model	Mean		Standard Deviation
		Depth = 0.5 m	Depth = 1.0 m	Depth = 0.5 m	Depth = 1.0 m
sp-flow	vG–M	−0.224	−0.164	0.409	0.115
	mvG–M	0.037	−0.071	0.057	0.027
	Kr–sp–FLAC	0.041	−0.055	0.065	0.11
CHM	vG–M	−0.223	−0.164	0.41	0.114
	mvG–M	0.03	−0.068	0.055	0.026
	Kr–sp–FLAC	0.036	−0.063	0.061	0.1

. If both the mean and standard deviation of the residuals are zero, the simulation results and the observed values are completely in agreement. The closer the mean and standard deviation of the residuals to zero, the better the prediction performance of a simulation. The application of the CHM model with the mvG–M model resulted in the smallest values of the mean (0.03) and the standard deviation (0.055) at a depth of 0.5 m. When the sp-flow model was applied with the K_r–sp–FLAC model, the mean value (−0.055) was smallest at a depth of 1.0 m, but the standard deviation (0.11) was quite larger than the smallest standard deviation (0.026) obtained from the CHM model with the mvG–M model. The large standard deviations from the K_r–sp–FLAC model resulted from the relatively large differences between the measured and simulated pore water pressures compared to the results from the mvG–M model, as shown in Figure 7b,d and Figure 8b. Therefore, it was concluded that the CHM model with the mvG–M model could be used for simulating a process of infiltration for the study site appropriately.

4.2. Slope Stability Assessment

The rainfall infiltration was simulated with applications of the CHM model and the mvG–M model to identify characteristics of the infiltration behavior of the study site that depended on rainfall intensity. Three rainfall scenarios up to a cumulative rainfall of 400 mm were applied: with a constant rainfall intensity of 5 mm/h that was smaller than the saturated hydraulic conductivity ($k_s=1.97\times {10}^{-6}$ m/s); with a constant rainfall intensity of 7.1 mm/h that was similar to $k_s$; and with a constant rainfall intensity of 10 mm/h that was larger than $k_s$. Figure 9a,c,e shows vertical profiles of the matric suction against time variation at section C–H in the slope when rainfall intensities of 5, 7.1, and 10 mm/h were applied, respectively. When the rainfall intensity is greater than the saturated hydraulic conductivity, air flow disturbs rainfall infiltration [21]. Therefore, the greater the rainfall intensity, the smaller the reduction of matric suction at the lower part of the slope. However, the decrease of matric suction at the upper part of the slope was large when the large rainfall intensity was applied because hydrologic responses to rainfall are relatively fast at the upper part of the slope. Influences of the rainfall intensity and the slope angle were examined by changes in the minimum safety factor at section C–H against cumulative rainfall variations (up to 400 mm) varying the angle of the slope to 25°, 30°, and 34°, as shown in Figure 9b,d,f. It is observed that rainfall intensity and slope angle have little effect on the decrease rate in the safety factor during the rainfall with a cumulative rainfall of 400 mm. When the large rainfall intensity and the small slope angle were applied, the reduction rate of the safety factor increased only marginally.

The representative hydraulic properties of soils utilized in the study site were determined based on a limited number of soil samples, although material properties vary depending on complex geological features and the structure of the strata. Influences of the air entry pressure ($P_0$), van Genuchten SWRC coefficient ($a$), and saturated hydraulic conductivity ($k_s$) on the reduction in the safety factor were examined when rainfall was applied. Figure 10 exhibits changes in the minimum safety factor at section C–H against time variations when the material properties summarized in Table 1 were consistently applied, with the exception of changing only the values of air entry pressure or van Genuchten SWRC coefficient or saturated hydraulic conductivity. In general, safety factors decreased slowly because of the small saturated hydraulic conductivity (from $1.97\times {10}^{-6}$ to $8\times {10}^{-6}$ m/s). When the air entry pressure increased from 3 to 9 kPa, the initial safety factor increased from 1.058 to 1.236, as shown in Figure 10a. The increase in air entry pressure resulted in an increase in the saturation of water and the volumetric water content at the same matric suction, as shown in Figure 10b, and the pore fluid pressure ($S_w{P}_w+S_g{P}_g$) decreased sequentially because of negative values of the initial pore water pressure, causing the increase in the safety factor. The larger the value of air entry pressure, the larger the reduction rate of the safety factor because infiltration rates increased due to the saturation of water increased. When the SWRC coefficient decreased (from 0.52 to 0.46), the initial safety factor and the reduction rate of the safety factor increased, as shown in Figure 10c. This is because effects of decreasing the van Genuchten SWRC coefficient on the form of SWRC are similar to effects of increasing the air entry pressure on that, as shown in Figure 10d. It can be observed from Figure 10e that when the large saturated hydraulic conductivity was applied, the increased rate of infiltration resulted in the fast reduction in the safety factor.

Slope stability assessments were performed for surroundings of the study site based on the GIS-based topography utilizing the digital elevation model (DEM) with a cell size of 10 m modeled from 1:5000 digital maps that were provided by the National Geographic Information Institute of Korea. It was assumed that rainfall was applied to the infiltration simulation with a constant rainfall intensity of 10 mm/h, and the CHM model and the mvG–M model were applied with the material properties given in Table 1 to the infinite slope models with varying slope angles considering the range of slope angles from the GIS-based topography. The pore pressure profiles against variations in cumulative rainfall obtained from the infiltration simulation were used for slope stability analyses. The minimum safety factors were computed for the infinite slope models that were determined to be safety factors of the cells of the GIS-based topography slope angles, both of which and the infinite slope models were the same. Figure 11 displays the distribution of safety factors when cumulative rainfall reached 100, 200, 300, and 400 mm. As the cumulative rainfall increased from 100 mm to 200, 300, and 400 mm (from 10 h to 20, 30, and 40 h with a rainfall intensity of 10 mm/h), the ratio of the slope failure prediction area (i.e., safety factor < 1) to a total area of the study site increased from 8.11% (81,025 m²/999,650 m²) to 8.31% (83,025 m²/999,650 m²), 8.46% (84,575 m²/999,650 m²), and 8.66% (86,580 m²/999,650 m²), respectively. Although additional rainfall of 300 mm was applied, the slope failure prediction area increased slightly due to the hydrological characteristics of low infiltration rates in the study site. It can be concluded that the study site is relatively invulnerable to slope failures, and this is comparable to the fact that no slope failure occurred in the study site on July 16, 2009 during the rainfall of approximately 400 mm, whereas slope failures occurred at 30 sites in the other areas of Hwangryeong Mountain [57].

5. Conclusions

In this study, an appropriate approach to a rainfall infiltration analysis for the study site, located in Hwangryeong Mountain, Busan, South Korea, was investigated. The change in pore water pressure in response to rainfall was measured, and material properties of soils in the study site were acquired through field investigations and lab tests. Variations in pore water pressure were simulated through the infiltration analyses in which three different types of relative permeability models and two different fluid flow models were applied. The simulated pore water pressures were compared with the measured pore water pressures to determine the appropriate site-specific approach for modeling rainfall infiltration in the study site. Compared to the fluid flow model, the relative permeability model significantly affected the increase rate of pore water pressure. In particular, the mvG–M model outperformed the other relative permeability models during infiltration analysis for the study site. As to the results, an application of the mvG–M model and the CHM model resulted in the best simulation with small residuals.

Based on the site-specific approach and characteristics, the effects of rainfall intensity, slope angle, factors in relation to SWRC (i.e., air entry pressure and van Genuchten SWRC coefficient), and saturated hydraulic conductivity on the slope stability were examined when the rainfall with constant intensities was applied. Rainfall intensity and slope angle did not significantly affect the change in the safety factor because of quite slow infiltration rates that resulted from specific hydraulic properties of the study site, whereas the factors concerning SWRC affecting the form of SWRC and saturated hydraulic conductivity influenced the initial safety factor and the reduction rate of the safety factor. It is necessary to apply appropriate site-specific hydraulic properties to perform the slope stability analysis adequately.

Maps of safety factors around the study site were generated that depended on cumulative rainfall based on the site-specific characteristics with the application of rainfall with a constant intensity. Although real-time monitoring of slope failure can be helpful in quickly and accurately identifying where a slope failure will occur, it is difficult to prepare countermeasures for disaster prevention against the slope failure in advance because of the absence of information on the slope failure prediction area. Safety factor maps in response to numerous cumulative rainfalls can be built in regions vulnerable to landslides for various rainfall scenarios through performing slope stability analyses based on site-specific characteristics appropriate to the area. Given that rainfall forecast information (from short- to mid-term) can be easily obtained from the Meteorological Administration, the rainfall-based safety factor map can be useful for managing the regions vulnerable to landslides and for preparing countermeasures for disaster prevention. Particularly in urban areas where mountains vulnerable to landslides are located, the procedure proposed in this study, which helps to improve prevention and response to landslide disasters, can be usefully applied for enhancing local sustainability.