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

: 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, inﬂuential factors, such as type of relative permeability models, air ﬂow, and variation of hydraulic conductivity associated with stress–strain behavior of soil, have signiﬁcant e ﬀ ects on the actual mechanism of rainfall inﬁltration. 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 inﬁltration analysis with applications of three relative permeability models, air ﬂow, and the coupled hydro-mechanical analysis to examine an appropriate site-speciﬁc approach to a rainfall inﬁltration analysis. The inﬁltration and stability analyses based on the site-speciﬁc 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.


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].

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.
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 = (P C /P 0 ) 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 r w is the residual saturation, θ is the actual volumetric water content, and θ r and θ s are the residual volumetric water content and saturated volumetric water content, respectively.
The matric suction (ψ) is derived by using Equation (3) as follows: 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.
Sustainability 2020, 12, 2839 5 of 21 q g i = −k w ij µ w µ g k g r ∂ ∂x j P g − ρ g g k x k , (6) where k w ij is the saturated mobility coefficient, k r is the relative permeability which is a function of the saturation of the wetting fluid, µ w µ g is the viscosity ratio between the two fluids, ρ w and ρ 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).
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).
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 (=0.02 P 0 ), was introduced to modify the SWRC through tangentially extrapolating from (P 0 , S e0 ) 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.
where λ = where λ = 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 nonwetting fluid.
where is the porosity, is the fluid bulk modulus, is the time, and is the volumetric strain.

Balance of Momentum
The balance of momentum is described as Equation (17).
where is the total stress, is the bulk density, and is the velocity. In unsaturated soil conditions, the bulk density is described as Equation (18). 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.
where n is the porosity, K is the fluid bulk modulus, t is the time, and ε is the volumetric strain.

Balance of Momentum
The balance of momentum is described as Equation (17).
where σ is the total stress, ρ is the bulk density, and .
u is the velocity. In unsaturated soil conditions, the bulk density is described as Equation (18).
where ρ d is the dry density of soil, ρ w is the density of water, and ρ 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).
where τ max is the shear strength, σ b is Bishop's effective stress, φ is the effective friction angle, c is the effective cohesion, andis the effective cohesion, and χ 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 χ with the degree of saturation, Equation (21) was derived from Equation (20). Then, Equation (19) may be expressed as Equation (22).
where S w P w + S g P g is the pore fluid pressure.

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.
where k s0 is the initial hydraulic conductivity, and n 0 is the initial porosity.

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)).
where τ f is the shear strength of soil, τ m is the shear stress on failure surface, β 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.
where z w is the depth to the potential failure surface, γ t is the total unit weight of the soil, γ d is the dry unit weight of the soil, γ w is the unit weight of water, and S w (z) is the saturation based on depth.
Sustainability 2020, 12, x FOR PEER REVIEW 8 of 21 where is the depth to the potential failure surface, is the total unit weight of the soil, is the dry unit weight of the soil, is the unit weight of water, and ( is the saturation based on depth.

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 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 3 with an average of 1220 kg/m 3 , 0.1-6.6 kPa with an average of 2 kPa, 32-39° with an average of 36°, and 1.9 × 10 -5.2 × 10 m/s with an average of 5.8 × 10 m/s, respectively). The dry density (1250 kg/m 3 ), cohesion (0.1 kPa), friction angle (33°), and hydraulic conductivity (1.97 × 10 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. Table 1. Material properties of soil, water, and air.

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 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 3 with an average of 1220 kg/m 3 , 0.1-6.6 kPa with an average of 2 kPa, 32-39 • with an average of 36 • , and 1.9 × 10 −6 -5.2 × 10 −4 m/s with an average of 5.8 × 10 −5 m/s, respectively). The dry density (1250 kg/m 3 ), cohesion (0.1 kPa), friction angle (33 • ), and hydraulic conductivity (1.97 × 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.

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 3 , 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 spflow option, the discretization of the domain, initial conditions, and boundary conditions were

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 3 , 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.

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 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.

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 Kr-sp-FLAC model). The pore water pressures were comparable with the field observations when the mvG-M model and the Kr-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 Kr-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 Kr-sp-FLAC model. During weak rainfall (i.e., 6-15 h), however, the   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 × 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. Residuals between the measured and estimated pore water pressures ( and , respectively) were computed by using an equation, Residual = ln( / , to evaluate the prediction accuracy of the models. Figure 8   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 (Pw mea and Pw est , respectively) were computed by using an equation, Residual = ln(Pw mea /Pw est ), to evaluate the prediction accuracy of the models. Figure 8  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.

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 × 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 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 × 10 −6 to 8 × 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. 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 ( + ) 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.  (Table 1) were applied with the exception of changing only the values of: (a) air entry pressure (from 3 to 9 kPa); (c) van Genuchten SWRC coefficient (from 0.52 to 0.46); (e) saturated hydraulic conductivity (from 1.97 × 10 to 8 × 10 m/s). Soil-water characteristic curves that depended on variations in: (b) air entry pressure (from 3 to 9 kPa); (d) van Genuchten SWRC coefficient (from 0.52 to 0.46).  (Table 1) were applied with the exception of changing only the values of: (a) air entry pressure (from 3 to 9 kPa); (c) van Genuchten SWRC coefficient (from 0.52 to 0.46); (e) saturated hydraulic conductivity (from 1.97 × 10 −6 to 8 × 10 −6 m/s). Soil-water characteristic curves that depended on variations in: (b) air entry pressure (from 3 to 9 kPa); (d) van Genuchten SWRC coefficient (from 0.52 to 0.46). 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 2 /999,650 m 2 ) to 8.31% (83,025 m 2 /999,650 m 2 ), 8.46% (84,575 m 2 /999,650 m 2 ), and 8.66% (86,580 m 2 /999,650 m 2 ), 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]. 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 2 /999,650 m 2 ) to 8.31% (83,025 m 2 /999,650 m 2 ), 8.46% (84,575 m 2 /999,650 m 2 ), and 8.66% (86,580 m 2 /999,650 m 2 ), 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].

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.