E ﬀ ects of Two-Phase Flow of Water and Air on Shallow Slope Failures Induced by Rainfall: Insights from Slope Stability Assessment at a Regional Scale

: Over 160 shallow landslides resulted from heavy rainfall that occurred in 26–27 July 2011 at Umyeon Mountain, Seoul, South Korea. To accurately reﬂect the ﬂuid ﬂow mechanism in the void spaces of soils, we considered the two-phase ﬂow of water and air for rainfall inﬁltration analysis using available historical rainfall data, topographic maps, and geotechnical / hydrological properties. Variations in pore water and air pressure from the inﬁltration analysis are used for slope stability assessment. By comparing the results from numerical models applying single- and two-phase ﬂow models, we observed that air ﬂow changes the rate of increase in pore water pressure, inﬂuencing the safety factor on slopes with a low inﬁltration capacity, where ponding is more likely to occur during heavy rainfall. Finally, several slope failure assessments were conducted to evaluate the usefulness of using the two-phase ﬂow model in forecasting slope stability in conditions of increased rainfall sums. We observed that the two-phase ﬂow model reduces the tendency of over-prediction compared to the single-phase model. The results from the two-phase ﬂow model revealed good agreement with actual landslide events.


Introduction
Rainfall-induced landslides have caused many human casualties and severe property damage worldwide [1,2]. The initiation of rainfall-induced landslides is affected by various factors, such as the material properties of soil, topographic features, geological and hydrological characteristics, slope surface cover, and vegetation [3,4]. Slope failures are caused by not only the infiltration of water into the soil of slopes during heavy rainfall but also the groundwater rise or the lateral water flow in the underground of slopes after precipitation, triggering buoyancy forces owing to changes in water pressure near the slope foot [5][6][7][8]. During the rainfall infiltration, the increase in water content causes an increase in the unit weight of the soil [9,10], and the loss of matric suction as wetting bands progress into the slope decreases its shear strength and effective stress [11,12].
Although numerous empirically based approaches have been widely applied to large areas using simple methods, deterministic physics-based approaches are preferable for a better understanding of the mechanisms of shallow landslides [13]. Numerous physically based studies have investigated infiltration behavior during rainfall using various approaches: (1) Adopting assumptions for simplification to solve the differential equation of Richards [14] by applying a single-phase fluid (water) flow [3,10,15,16]; (2) assuming that the air pressure in the void spaces of an unsaturated soil slope is equal to the atmospheric pressure [9,17,18]; (3) applying a two-phase fluid (water and air) flow [19][20][21][22][23]; and (4) considering hydraulic hysteresis, including the effects of capillarity and/or air entrapment, reflecting different hydraulic states under wetting and drying conditions [24,25]. Such physically based numerical studies have mainly been applied at slope or watershed scales.
Slope stability assessments at a regional or local scale are required to simulate landslide events. Simoni et al. [26] simulated shallow landslides using GEOtop-FS [27], which models subsurface flow based on the numerically integrated 3D Richards equation and assesses slope stability depending on a specified topography. Chen et al. [28] utilized Scoops3D [29] to evaluate landslide events in Hualien, Taiwan. However, such simulation in a 3D scheme incurs an excessive amount of computation [16]. Park et al. [30], Jeong et al. [15], Tran et al. [16], and Tran et al. [31] assessed the slope failures that occurred at Umyeon Mountain in Seoul, South Korea in 2011 using SEEP/W and SLOPE/W [32], TRIGRS [33], TiVaSS [34], and Scoops3D [29], respectively. However, air flow was not considered. As the ground consists of materials in three phases (soil particles, water, and air), it is necessary to consider two-phase flow to interpret rainfall infiltration accurately [11]. Several experiments and interpretations have reported that the air flow generated during infiltration changes the rate of the infiltrating water, affecting the safety factor. Therefore, it is apparent that the effect of air on the subsurface flow of fluids is relevant and should be included in infiltration analysis [21,22,35,36].
To efficiently utilize computing resources with application of a two-phase flow model based on numerical methods at a regional scale, each slope based on the geographic information system (GIS)-based topography of Umyeon Mountain located in Seoul, South Korea was simplified to be an infinite slope in a two-dimensional (2D) domain. We then applied the two-phase flow model to the infinite slopes for simulations of the rainfall-induced slope failure. Based on the variations in pore air and water pressures obtained from the simulation of rainfall infiltration, slope failure analyses were performed to compute safety factors at each infinite slope model. The minimum safety factor of elements in the infinite slope was considered to be a safety factor at the corresponding GIS-based topography cell. We evaluated the performance of the analysis for slope failure assessments applying the two-phase flow model. This study is concerned with (1) the investigation of the effects of air flow and pressure on infiltration behavior by comparing the variations in the pore water pressure predicted by the single-phase flow and two-phase flow models, (2) interpretations of the influence of air flow and pressure on safety factors during heavy rain, (3) comparison of the landslide susceptibility maps that were obtained by applying the single-phase flow and two-phase flow models, and (4) assessments of the applicability of single-phase flow and two-phase flow models by comparing areas that are predicted as unstable with the locations of actual landslide events.

Study Area
Umyeon Mountain is in Seoul, South Korea, at 37 • 27 00"-37 • 29 55" N and 126 • 59 02"-127 • 01 41" E. It has an area of approximately 5 km 2 . Figure 1 shows the elevation distribution of the study area based on 1:5000 digital maps produced by the National Geographic Information Institute of Korea [37] using the airborne laser terrain mapping system. The highest elevation is 295 m above sea level. The area around Umyeon Mountain consists mainly of Pre-Cambrian banded biotite gneiss and granitic gneiss [38]. According to Jeong et al. [15], the bedrock is primarily intensely weathered and substantially fractured plagioclase, quartz, biotite, feldspar, and amphibole. Umyeon Mountain is completely surrounded by urban development. The weather station closest to Umyeon Mountain is Namhyun station. It is approximately 2 km away and is included in the automatic warning system (AWS) operated by the Korea Meteorological Administration (KMA) [40]. A precipitation of 1120 mm (i.e., approximately 55% of the normal annual precipitation in 2011) occurred in the area of Umyeon Mountain over July 2011. A maximum hourly rainfall of 114 mm/h occurred from 07:44 to 08:44 on July 27. Figure 2 shows the time series for hourly and cumulative rainfall recorded at the Namhyun station during the two days from July 26 to July 27, 2011. Continuous rainfall of 320 mm occurred during the 8 h prior to the landslide events with strong rainfall intensities (higher than 17.5 mm/h). During the 2 h immediately prior to the landslide events, the rainfall attained 167 mm (i.e., 43% of the cumulative rainfall during the two days from July 26 to July 27, 2011).

Landslide Inventories and Spatial Data
Landslide inventories of information on slope failures are important for validating the applicability of slope stability assessments. We compared satellite images with a resolution of 5 m and aerial photographs with a resolution of 25 cm to map the slope failure sites extracted from the study area before and after the 2011 landslides. A total of 161 slope failure initiation sites were constructed in a GIS format (Figure 3a). We also assessed the suitability of the available landslide inventories by comparing them with field survey-based data [41].
GIS-based topographic information is also required for slope stability assessments at a regional scale. We used a digital elevation model (DEM) with a 10 m-cell size computed from 1:5000 digital maps produced by the National Geographic Information Institute of Korea [37] (see Figure 1). Approximately 160 catastrophic shallow landslides occurred in this area owing to heavy rain on 26-27 July 2011, with a cumulative total of 470 mm over the two days [39]. Debris flows that resulted from slope failures occurred along 20 watersheds from 7:00 to 9:00 local time on 27 July 2011, resulting in 16 deaths and severe structural damage to 10 buildings. The economic loss caused by the landslides was approximately 15 million US dollars [30].
The weather station closest to Umyeon Mountain is Namhyun station. It is approximately 2 km away and is included in the automatic warning system (AWS) operated by the Korea Meteorological Administration (KMA) [40]. A precipitation of 1120 mm (i.e., approximately 55% of the normal annual precipitation in 2011) occurred in the area of Umyeon Mountain over July 2011. A maximum hourly rainfall of 114 mm/h occurred from 07:44 to 08:44 on 27 July. Figure 2 shows the time series for hourly and cumulative rainfall recorded at the Namhyun station during the two days from 26 July to 27 July 2011. Continuous rainfall of 320 mm occurred during the 8 h prior to the landslide events with strong rainfall intensities (higher than 17.5 mm/h). During the 2 h immediately prior to the landslide events, the rainfall attained 167 mm (i.e., 43% of the cumulative rainfall during the two days from 26 July to 27 July 2011).

Landslide Inventories and Spatial Data
Landslide inventories of information on slope failures are important for validating the applicability of slope stability assessments. We compared satellite images with a resolution of 5 m and aerial photographs with a resolution of 25 cm to map the slope failure sites extracted from the study area before and after the 2011 landslides. A total of 161 slope failure initiation sites were constructed in Water 2020, 12, 812 4 of 20 a GIS format (Figure 3a). We also assessed the suitability of the available landslide inventories by comparing them with field survey-based data [41].
GIS-based topographic information is also required for slope stability assessments at a regional scale. We used a digital elevation model (DEM) with a 10 m-cell size computed from 1:5000 digital maps produced by the National Geographic Information Institute of Korea [37] (see Figure 1).
Water 2020, 12, x FOR PEER REVIEW 4 of 19 Figure 2. Hourly rainfall and cumulative rainfall recorded at the Namhyun station during the two days from July 26 to July 27, 2011 (obtained from KMA).

Material Properties
Representative geological units in the area of Umyeon Mountain are PCEbngn (Pre-Cambrian Era banded biotite gneiss), PCEggn (Pre-Cambrian Era granitic gneiss), and Jdgr (Jurassic daebo granite), as shown in Figure 3b. Pre-CambrianEra banded biotite gneiss, the lowest strata, is widely distributed in the study area, and Pre-Cambrian Era granitic gneiss and Jurassic daebo granite are distributed in relatively small areas near the top and western part of Umyeon Mountain, respectively. Jurassic daebo granite formed in the Mesozoic Jurassic period and has intruded into Pre-Cambrian Era banded biotite gneiss. The rocks that form the bedrock are slightly weathered or moderately weathered, and the rocks at the top of the mountain are more weathered than those at the bottom of the mountain [41]. Colluvium underlain by weathered soils is distributed in the entire area of Umyeon Mountain with unconformity, composed of loose silty sand mixed with gravel [42]. Figure  3c shows the boring logs and groundwater levels obtained from three different sites (BL-1, BL-2, and BL-3 shown in Figure 3b). The groundwater tables were distributed underneath the colluvium layer. The Korean Geotechnical Society [42] reported that most of the slope failures occurred at shallow depths with a range from 0.5 to 2 m and that depths of colluvium near the slope failure sites ranged from 0.3 to 3 m. We assumed that slope failure could occur in the colluvium layer with constant geotechnical properties. Infiltration and slope stability analyses were then performed with a focus more on the geotechnical characteristics than on the geological structure.
The Seoul Metropolitan Government [41] carried out detailed geotechnical field investigations at Umyeon Mountain by drilling 58 boreholes. Park et al. [30] used the geotechnical properties available from 13 boreholes and divided Umyeon Mountain into five zones based on the distribution of boreholes, geotechnical properties, and watersheds ( Figure 3b). The cohesion, friction angle, and saturated hydraulic conductivity at borehole sites in zone 1 were in the ranges of 8.4-11.9 kPa with an average of 9.6 kPa, 24.8-26° with an average of 25.3°, and 6.67 × 10 −6 -8.08 × 10 −6 m/s with an average of 7.15 × 10 −6 m/s, respectively, and those at boreholes in zone 5 were in the ranges of 5.1-7.5 kPa with an average of 6.3 kPa, 27.6-30.3° with an average of 28.2°, and 2.08 × 10 −6 -4.45 × 10 −6 m/s with an average of 3.69 × 10 −6 m/s, respectively. Considering that the variations of them were not large and that simplification is required for simulations at a regional scale, we determined that

Material Properties
Representative geological units in the area of Umyeon Mountain are PCEbngn (Pre-Cambrian Era banded biotite gneiss), PCEggn (Pre-Cambrian Era granitic gneiss), and Jdgr (Jurassic daebo granite), as shown in Figure 3b. Pre-CambrianEra banded biotite gneiss, the lowest strata, is widely distributed in the study area, and Pre-Cambrian Era granitic gneiss and Jurassic daebo granite are distributed in relatively small areas near the top and western part of Umyeon Mountain, respectively. Jurassic daebo granite formed in the Mesozoic Jurassic period and has intruded into Pre-Cambrian Era banded biotite gneiss. The rocks that form the bedrock are slightly weathered or moderately weathered, and the rocks at the top of the mountain are more weathered than those at the bottom of the mountain [41]. Colluvium underlain by weathered soils is distributed in the entire area of Umyeon Mountain with unconformity, composed of loose silty sand mixed with gravel [42]. Figure 3c shows the boring logs and groundwater levels obtained from three different sites (BL-1, BL-2, and BL-3 shown in Figure 3b). The groundwater tables were distributed underneath the colluvium layer. The Korean Geotechnical Society [42] reported that most of the slope failures occurred at shallow depths with a range from 0.5 to 2 m and that depths of colluvium near the slope failure sites ranged from 0.3 to 3 m. We assumed that slope failure could occur in the colluvium layer with constant geotechnical properties. Infiltration and slope stability analyses were then performed with a focus more on the geotechnical characteristics than on the geological structure.
The Seoul Metropolitan Government [41] carried out detailed geotechnical field investigations at Umyeon Mountain by drilling 58 boreholes. Park et al. [30] used the geotechnical properties available from 13 boreholes and divided Umyeon Mountain into five zones based on the distribution of boreholes, geotechnical properties, and watersheds ( Figure 3b). The cohesion, friction angle, and saturated hydraulic conductivity at borehole sites in zone 1 were in the ranges of 8.4-11.9 kPa with an average of 9.6 kPa, 24.8-26 • with an average of 25.3 • , and 6.67 × 10 −6 -8.08 × 10 −6 m/s with an average of 7.15 × 10 −6 m/s, respectively, and those at boreholes in zone 5 were in the ranges of 5.1-7.5 kPa with an average of 6.3 kPa, 27.6-30.3 • with an average of 28.2 • , and 2.08 × 10 −6 -4.45 × 10 −6 m/s with an average of 3.69 × 10 −6 m/s, respectively. Considering that the variations of them were not large and that simplification is required for simulations at a regional scale, we determined that the zones, categorized in Park et al. [30], could be appropriately applied in this study. Table 1 summarizes the soil properties (internal friction angle, cohesion, total/dry unit weight of soil, and saturated hydraulic conductivity) obtained from Park et al. [30] and the soil-water characteristics (i.e., soil-water characteristic curve (SWCC) [43] coefficients) obtained from the Seoul Metropolitan Government [41] in which the pressure plate method (GTCS) [44] and the filter paper method were applied. The representative soil properties are determined as the average values of the boreholes located in each zone. We applied two van Genuchten SWCCs, generally used in numerical methods, to the north side (zones 1 and 2) and south side (zones 3, 4, and 5) of Umyeon Mountain (Figure 4a). Figure 4b shows the Mualem-van Genuchten relative permeability curves [43] applied to the north and south areas in this study. The material properties of water and air applied to the infiltration analysis are presented in Table 2. soil-water characteristic curve (SWCC) [43] coefficients) obtained from the Seoul Metropolitan Government [41] in which the pressure plate method (GTCS) [44] and the filter paper method were applied. The representative soil properties are determined as the average values of the boreholes located in each zone. We applied two van Genuchten SWCCs, generally used in numerical methods, to the north side (zones 1 and 2) and south side (zones 3, 4, and 5) of Umyeon Mountain (Figure 4a). Figure 4b shows the Mualem-van Genuchten relative permeability curves [43] applied to the north and south areas in this study. The material properties of water and air applied to the infiltration analysis are presented in Table 2. (c) Boring logs and groundwater levels obtained at three sites shown in (b) (modified from the Seoul Metropolitan Government [41]). (d) Locations of the sites' selected borehole for the investigation conducted by the Seoul Metropolitan Government [41] and the distribution of the five zones proposed by Park et al. [30].   [30] and soil-water characteristics [41] for each zone used in the infiltration and slope failure analyses.   Figure 4a. The constants of the two van Genuchten SWCCs used to compute Mualem-van Genuchten relative permeability curves [43] are presented in Table 1. Table 2. Summary of the material properties of water and air used in the infiltration analysis.

Rainfall Infiltration Analysis
We used FLAC 2D Ver. 7.0 [46] to apply the two-phase flow of water and air for the infiltration analysis of unsaturated soil slopes. This simulates the flow of two immiscible fluids in a porous medium. The pore spaces are assumed to be completely filled with the two fluids. In unsaturated soils, water and air can be substituted for the wetting fluid (w) and non-wetting fluid (g), respectively. van Genuchten [43] proposed Equation (1), which describes a relationship between effective saturation (S e ) and matric suction (ψ = non-wetting fluid pressure (P g )-wetting fluid pressure (P w )) as follows: where S w is the actual saturation of the wetting fluid, and S r w is the residual saturation. Darcy's law allows calculation of how the wetting and non-wetting fluids are transported in the two-phase flow option. Equations (2) and (3) represent the velocities of the wetting and non-wetting fluids (q w i and q g i ), respectively: where k w ij is the saturated mobility coefficient (= k s /(gρ w )), k s is the saturated hydraulic conductivity, g is the gravitational acceleration, ρ w is the density of the wetting fluid, k r is the relative permeability, µ w / µ g is the viscosity ratio between the wetting and non-wetting fluids, and ρ g is the density of the non-wetting fluid.
The relative permeability of a wetting fluid is required to calculate the unsaturated hydraulic conductivity. van Genuchten [43] proposed the relative permeability (k w r ) based on the capillary model of Mualem [47] (Equation (4)). Lenhard and Parker [48] calculated the relative permeability of a non-wetting fluid (k g r ), expressed as in Equation (5): where b and c are constants. Both b and c are equal to 0.5 for water and air [11]. Equations (6) and (7) are the balance equations for the wetting and non-wetting fluids, respectively. These equations can be derived by combining the fluid balance laws with the fluid constitutive laws for the two fluids: where n is the porosity, K is the fluid bulk modulus, t is the time, and ε is the volumetric strain. The balance of momentum is described as in Equation (8): where σ is the total stress, ρ is the bulk density, and .
u is the velocity of a soil particle. The bulk density is described by Equation (9) for unsaturated soil conditions: where ρ d is the dry density of soil. We used Bishop's effective stress [49] to analyze the failure of unsaturated soil in the Mohr-Coulomb constitutive model. Equation (10) presents the yield criterion of the Mohr-Coulomb constitutive model. Bishop's effective stress (σ b ) is expressed by Equation (11) as below: where τ max is the shear strength, φ is the effective friction angle, c is the effective cohesion, and χ is the matric suction coefficient that varies from zero (i.e., dry conditions) to one (i.e., fully saturated conditions). The degree of saturation can be approximately substituted for the matric suction coefficient [11,19,[50][51][52]. Equation (12) was derived from Equation (11) by substituting the degree of saturation for χ. Equation (10) is then expressed as Equation (13): where S w P w + S g P g is the pore fluid pressure. The two-phase flow model is composed of the fluid flow and mechanical loops. The fluid flow loop computes fluid flows based on pressure gradients and variations in saturation and pore water and air pressures induced by unbalanced flows, utilizing from Equation (2) to Equation (7). The mechanical loop computes total stress, which depends on the coordinates, velocities, and pore pressure generation induced by mechanical volume strain, utilizing from Equation (8) to Equation (13). The stress state from the mechanical loop is used for the next time step of the fluid flow loop.

Slope Stability Analysis
The infinite slope analysis method can be used to interpret rainfall-induced shallow slope failure in Korea. This is because the failure surface is generally parallel to the slope surface in this region [3]. We applied the limit equilibrium method to the infinite slope. The safety factor of the infinite slope (F s ) can be expressed as in Equation (14), given the Mohr-Coulomb failure criterion (Equations (10) and (13)) and Bishop's effective stress (Equations (11) and (12)): where τ f is the shear strength of the soil, τ m is the shear stress on the failure surface, β is the slope angle, and W is the weight of a soil slice. The total unit weight of soil (W) depends on variations in saturation caused by rainfall infiltration. To consider an increase in the unit weight of soil during rainfall, the weight of a soil slice per unit area above a potential failure surface can be described as in Equation (15): where z w . is the depth from the slope surface to the potential failure surface, γ t is the total unit weight of soil, γ d is the dry unit weight of soil, γ w is the unit weight of water, and S w (z) is the saturation, which depends on depth.

Slope Geometry and Boundary/Initial Conditions
In Korea, the depth at which shallow slope failure surfaces are present generally ranges from 1 to 3 m [53], and such failure surfaces are parallel to the ground surface. Therefore, to assess slope failure, Park et al. [30] assumed a uniform soil depth of 2 m for Umyeon Mountain by considering an investigation of the distribution of soil depths by the Korean Geotechnical Korean Geotechnical Society [42]. We also modeled 2D infinite slopes with a uniform soil depth of 2 m and a length of 10 m, and we assumed the ground to be isotropic and homogeneous considering the engineering soil properties summarized in Table 1. Figure 5 illustrates an example of an infinite slope model with an angle of 31 • , which is approximately equal to the mean slope angle of the slope failure sites in the study area (i.e., 31.2 • ). We applied the two-phase flow model for simulating rainfall infiltration to the independent infinite slope models, where slope angles were considered based on different cells of the GIS-based slope raster in the study area. Variations in pore water and air pressures (P w and P a , respectively) and saturations with the time step could then be computed at the infinite slopes for 17 h during rainfall. The slope stability on the infinite slopes for the 17 h was calculated by applying the variations in P w , P a , and saturations computed from the two-phase flow model to Equation (14). Finally, we determined the safety factor of each cell of the GIS-based topography in the study area to be the minimum safety factor of the corresponding infinite slope model.
Finally, we determined the safety factor of each cell of the GIS-based topography in the study area to be the minimum safety factor of the corresponding infinite slope model.
Water flux was set as equal to the hourly rainfall shown in Figure 2 and was applied at the ground surface as a boundary condition. The bottom surface was set to be completely constrained, at which displacements in both horizontal and vertical directions were fixed, and vertical fluid flows were unaccepted. Displacements in the horizontal direction were fixed at the left and right sides. The slope surface was subjected to the atmospheric boundary conditions. We set the initial on the slope as equal to atmospheric pressure. The other boundaries were assumed to be airtight. The applied in the single-phase flow option was set to zero.
The matric suction measured in the study area during a rainfall season in 2012 was maintained at 10-20 kPa prior to rainfall [41]. Therefore, we set the initial slope conditions to be steady-state with an initial of −20 kPa (a constant value by depth). The initial can be set as constant near the surface grounds where the depth of the groundwater table is known to be deep [3,11].

Performance Evaluation of Slope Stability Assessment
To evaluate the quantitative prediction performance for the slope failure areas, a confusion matrix [54] can be utilized. We also applied the modified success rate (MSR) proposed by Huang and Kao [55]. The modified success rate (MSR) defined by Equation (16) can consider the prediction performance in both stable and unstable areas. Moreover, it exhibits the advantage of preventing over-and under-prediction. MSR ranges between zero to one. MSR = 0.5 × SR + 0.5 × successfully predicted stable cells total number of actual stable cells , where SR is the former success rate (= ). Water flux was set as equal to the hourly rainfall shown in Figure 2 and was applied at the ground surface as a boundary condition. The bottom surface was set to be completely constrained, at which displacements in both horizontal and vertical directions were fixed, and vertical fluid flows were unaccepted. Displacements in the horizontal direction were fixed at the left and right sides. The slope surface was subjected to the atmospheric boundary conditions. We set the initial P a on the slope as equal to atmospheric pressure. The other boundaries were assumed to be airtight. The P a applied in the single-phase flow option was set to zero.
The matric suction measured in the study area during a rainfall season in 2012 was maintained at 10-20 kPa prior to rainfall [41]. Therefore, we set the initial slope conditions to be steady-state with an initial P w of −20 kPa (a constant value by depth). The initial P w . can be set as constant near the surface grounds where the depth of the groundwater table is known to be deep [3,11].

Performance Evaluation of Slope Stability Assessment
To evaluate the quantitative prediction performance for the slope failure areas, a confusion matrix [54] can be utilized. We also applied the modified success rate (MSR) proposed by Huang and Kao [55]. The modified success rate (MSR) defined by Equation (16) can consider the prediction performance in both stable and unstable areas. Moreover, it exhibits the advantage of preventing overand under-prediction. MSR ranges between zero to one.
successfully predicted stable cells total number of actual stable cells , where SR is the former success rate (= number of successfully predicted slope failures total number of actual slope failures ).

Validation with Comparison to Field Measurements
To validate the two-phase flow model utilized for infiltration simulations, we compared the matric suction measured by the Seoul Metropolitan Government [41] in the study area (zone 1) from 0:00 on 29 June 2012 to 20:00 5 July 2012 and the suction simulated from the two-phase flow model ( Figure 6). We modeled an infinite slope with a depth of 2 m and an angle of 25 • , applying the engineering soil properties of zone 1 summarized in Table 1. Water flux was set as equal to the hourly rainfall shown in Figure 6. Other boundary conditions were set identical to those described in the section "slope geometry and boundary/initial conditions". Considering the initial matric suctions measured at depths of 30, 70, and 110 cm (as shown in Figure 6), we set the initial slope conditions to be hydrostatic, with an initial P w of −76 kPa (a constant value by depth). It can be observed from Figure 6 that the simulated matric suction at depths of 30, 70, and 110 cm were in good agreement with the field measurements. Therefore, we concluded that the two-phase flow model could be appropriately used to simulate rainfall infiltration in the study area.
the single-phase flow model. The Pw in both the models decreased at 9 h from the start owing to the negligible rainfall between 7 and 9 h (see Figure 2). After 12 h, the infiltration rates increased sharply in the single-phase flow model during heavy rainfall whereas those from the two-phase flow model increased relatively gradually during this period. This is because air pressure interrupted the water flow into the void spaces and was accompanied by ponding ( [20,35,56]). This shows clearly that air pressure impacts the infiltration rate particularly in the low permeable soils where ponding effects occur frequently.

Infiltration Analysis
We conducted infiltration analyses for 17 h (from 16:00 on July 26 to 9:00 on July 27), applying the single-phase flow and two-phase flow models to 2D infinite slopes. The representative material properties of the soil at the five zones (Figure 3b) are summarized in Table 1. Figure 7 presents an example of variations in P w and matric suction with time at the five zones when the single-phase flow model was applied to the slope of 30 • , which is approximately equal to the mean slope angle of the slope failure sites in the study area (i.e., 31.2 • ). The variations in P w , P a , and matric suction from the two-phase flow model are shown in Figure 8. In the early stages (less than 9 h), variations in P w in the single-phase flow and two-phase flow models are similar at zones 1 and 4; ponding effects are negligible because of a high infiltration capacity with relatively high saturated hydraulic conductivities (see Table 1). However, at zones 2, 3, and 5, the ponding caused by a rainfall intensity higher than the infiltration capacity resulted in an increase in P a under the two-phase flow model. The increases in P a caused P w to increase even at certain depths under the slope where the rainfall had not reached. Therefore, increases in P w in the two-phase flow model became higher than those in the single-phase flow model. The P w in both the models decreased at 9 h from the start owing to the negligible rainfall between 7 and 9 h (see Figure 2). After 12 h, the infiltration rates increased sharply in the single-phase flow model during heavy rainfall whereas those from the two-phase flow model increased relatively gradually during this period. This is because air pressure interrupted the water flow into the void spaces and was accompanied by ponding ( [20,35,56]). This shows clearly that air pressure impacts the infiltration rate particularly in the low permeable soils where ponding effects occur frequently.  Variations in pore water pressure and matric suction with time when the single-phase flow model was applied to a slope of 30°. The vertical axis of each graph represents the depth from the ground surface. The first and second columns represent the ranges of the pore water pressure and matric suction, respectively.
In the single-phase flow model, the matric suction sequentially dissipated from the surface without the inclusion of Pa (see Figure 7). Meanwhile, in the two-phase flow model, the matric suction was marginal at shallow depths. However, it did not dissipate owing to Pa being higher than Pw (see Figure 8).  Because the relative permeabilities at the north side were too marginal at matric suction below 20 kPa (Figure 4b), increases in Pw at the north side (zones 1 and 2) were slower than those at the south side (zones 3, 4, and 5) in both the models. An exception was zone 3, where the saturated hydraulic conductivity was the lowest (1.80 × 10 −6 m/s). The infiltration rate was the highest at zone 4 owing to the largest relative permeability and saturated hydraulic conductivity (9.7 × 10 -6 m/s). The ground here was almost completely saturated after 15 h (7:00 on 27 July, 2 h prior to the landslide Figure 8. Variations in pore water/air pressures and matric suction with time when the two-phase flow model was applied to a slope of 30 • . The vertical axis of each graph represents the depth from the ground surface. The columns represent the ranges of the pore water pressure, pore air pressure, and matric suction from left to right. In the single-phase flow model, the matric suction sequentially dissipated from the surface without the inclusion of P a (see Figure 7). Meanwhile, in the two-phase flow model, the matric suction was marginal at shallow depths. However, it did not dissipate owing to P a being higher than P w (see Figure 8).
Because the relative permeabilities at the north side were too marginal at matric suction below 20 kPa (Figure 4b), increases in P w at the north side (zones 1 and 2) were slower than those at the south side (zones 3, 4, and 5) in both the models. An exception was zone 3, where the saturated hydraulic conductivity was the lowest (1.80 × 10 −6 m/s). The infiltration rate was the highest at zone 4 owing to the largest relative permeability and saturated hydraulic conductivity (9.7 × 10 -6 m/s). The ground here was almost completely saturated after 15 h (7:00 on 27 July, 2 h prior to the landslide events). The infiltration rate was the lowest at zone 2 owing to the marginal relative permeability and saturated hydraulic conductivity.
At the south side, P w increased toward the surface and larger depths, whereas at the north side, it increased first at shallow depths and subsequently at larger depths as the wetting bands progressed gradually into the depths because of the marginal relative permeability. In the two-phase flow model, positive P w occurred at shallow depths because of ponding effects and heavy rainfall (after 12 h). However, the matric suction was not dissipated completely in the shallow areas because here, P a is higher than P w (i.e., not fully saturated).

Slope Stability Assessment
The safety factor is dependent on shear strength parameters (i.e., internal friction angle and cohesion) and variations in P w and P a . Figure 9 shows the variations in the minimum safety factors on infinite slopes whose angles are the maximum, average, and minimum of the actual slope failure sites at each zone. The initial safety factors at zone 3 are larger than those at the other zones, with values higher than 3 (Figure 9c) owing to the large internal friction angle (37.6 • ) and high cohesion (7.7 kPa). Initial safety factors on slopes with the average angles at zones 1, 2, and 5, where the internal friction angles and cohesions are marginal (see Table 1), are smaller than two (Figure 9a,b,e). At the north side (zones 1 and 2), the safety factors decreased gradually owing to gradual infiltration rates and reduced abruptly after 13 h. Meanwhile, at the south side (zones 3, 4, and 5), the safety factors reduced gradually from the initial values and decreased to below one more rapidly than those at the north side. This is owing to the rapid infiltration (with an exception of zone 3 because of the smallest saturated hydraulic conductivity (1.8 × 10 −6 m/s)). Therefore, hydrological properties (i.e., saturated hydraulic conductivity and relative permeability) strongly affect variations in the safety factor, thereby influencing infiltration rates.
The safety factors from both the single-and two-phase flow models were similar during the early stages. However, the decrease in the safety factor in the two-phase flow model became higher than that in the single-phase flow model during the middle stages owing to the effects of ponding. Whereas the decrease in the effective stress resulted solely from the increase in P w in the single-phase flow model, the increases in both P w and P a in the two-phase flow model resulted in a higher decrease in the safety factor. This is because the shear strength decreased as P w and P a increased (referring to Equation (13)). However, the safety factor in the single-phase flow model became lower than that in the two-phase flow model during the final stage. This is because the slopes in the single-phase flow model became saturated more rapidly than those in the two-phase flow model (as shown in Figures 7 and 8). These observations are comparable with the observations of Cho [11].
We generated landslide susceptibility maps for Umyeon Mountain based on topographic maps (i.e., DEM and slope map) with a cell size of 10 m, the geotechnical properties summarized in Table 1, and the distributions of P w and P a obtained from infiltration analyses by applying the single-phase flow and two-phase flow models in FLAC (Figure 10a,b). The distribution patterns for the areas of predicted slope failure are similar because the same topographic maps and geotechnical properties were applied to generate both maps. However, applying different models resulted in differences in the total area predicted for slope failure, as shown in Figure 10c,d. On the north side (zones 1 and 2), the total area of potential slope failure in the two-phase flow model (234,600 m 2 ) is smaller than that in the single-phase flow model (405,525 m 2 ). This is because the rainfall infiltration at the north side is slower in the two-phase flow model than it is in the single-phase flow model. At zones 4 and 5, the slope failure areas predicted from the two-and single-phase flow models (1,018,925 and 1,189,575 m 2 , respectively) are similarly distributed because the ground approached saturated conditions owing to rapid infiltration with large relative permeabilities. However, under the application of the two-phase flow model, the potential slope failure areas are markedly marginal in Zone 3 because of the considerable matric suction there (as shown in Figure 8). This is comparable with the fact that there was no actual slope failure in zone 3 with a large representative internal friction angle of 37.6 • .
Water 2020, 12, x FOR PEER REVIEW 14 of 19 Figure 9. Variations in minimum safety factors with time when the two-and single-phase flow models were applied to slopes with maximum, average, and minimum slope angles for actual failure slopes: (a)-(e) for zones 1-5. We used the maximum, average, and minimum slope angles for the entire actual slope failure sites in the study area only for zone 3, where an actual slope failure did not occur. We then applied the modified success rate (MSR) to evaluate the quantitative prediction performance for the slope failure areas. In total, 161 actual slope failures occurred in the study area, and the total number of slope failures that were successfully predicted by the two-and single-phase flow models was 148 and 150, respectively. Therefore, the performance values produced by the former success rate (SR) for the three models are 0.92 (= 148/161) and 0.93 (= 150/161), respectively. The MSRs for both models are 0.86 and 0.84, respectively, and the stable cell coverage over the total watershed area is 80.1 and 74.7, respectively. Huang and Kao [55] reported that the best simulation could be achieved by MSR, where the stable cell coverage in a total watershed area ranges from 80% to 90%. Landslide over-prediction is likely to occur when the stable cell coverage is less than 80%, and under-prediction is likely when it is larger than 90%. We then applied the modified success rate (MSR) to evaluate the quantitative prediction performance for the slope failure areas. In total, 161 actual slope failures occurred in the study area, and the total number of slope failures that were successfully predicted by the two-and single-phase flow models was 148 and 150, respectively. Therefore, the performance values produced by the former success rate (SR) for the three models are 0.92 (= 148/161) and 0.93 (= 150/161), respectively. The MSRs for both models are 0.86 and 0.84, respectively, and the stable cell coverage over the total watershed area is 80.1 and 74.7, respectively. Huang and Kao [55] reported that the best simulation could be achieved by MSR, where the stable cell coverage in a total watershed area ranges from 80% to 90%. Landslide over-prediction is likely to occur when the stable cell coverage is less than 80%, and under-prediction is likely when it is larger than 90%. The largest MSR was derived from the two-phase flow model, where an adequate range of stable cell coverage was used for the best simulation (i.e., 80%-90%). The largest SR was derived from the single-phase flow model. However, the MSR-derived value was the second largest in this case, owing to the lesser stable cell coverage that probably causes an over-prediction of landslides. Table 3 summarizes the results of the performance evaluation of the two-and single-phase flow models utilizing the confusion matrix. True positive (TP), true negative (TN), false positive (FP), and false negative (FN) for the two-phase flow model are respectively 0.23%, 80.11%, 19.64%, and 0.02%, and those for the single-phase flow model are 0.24%, 74.69%, 25.05%, and 0.02%, respectively. Efficiency, positive predictive power, sensitivity, and specificity were computed from the confusion matrix to quantitatively evaluate the estimation performance. The efficiency means the proportion of correctly predicted class and is most frequently used in the literature [54]. The efficiency of the twophase flow model (80.19%) was greater than that of the single-phase flow model (74.79%). The positive predictive power means the proportion of TP in a total of positive predicted class. Low values of TP resulted from both the two-and single-phase flow models (0.23% and 0.24%, respectively). This is because a number of actual slope failure sites (161 cells) were quite small compared to the mountainous area considered in this study (63,157 cells). Even if all the actual slope failure sites are correctly predicted, the maximum value of TP is 0.25% (161/63,157). Therefore, regardless of the method of analysis, the value of the positive predictive power is consequently low. The positive predictive power of the two-phase flow model (1.17%) was slightly larger than that of the single-phase flow model (0.93%) due to the large FP, resulting in an over-prediction of landslides. Only for the sensitivity that is the same as the SR used to compute the MSR, the result from the single- The largest MSR was derived from the two-phase flow model, where an adequate range of stable cell coverage was used for the best simulation (i.e., 80%-90%). The largest SR was derived from the single-phase flow model. However, the MSR-derived value was the second largest in this case, owing to the lesser stable cell coverage that probably causes an over-prediction of landslides. Table 3 summarizes the results of the performance evaluation of the two-and single-phase flow models utilizing the confusion matrix. True positive (TP), true negative (TN), false positive (FP), and false negative (FN) for the two-phase flow model are respectively 0.23%, 80.11%, 19.64%, and 0.02%, and those for the single-phase flow model are 0.24%, 74.69%, 25.05%, and 0.02%, respectively. Efficiency, positive predictive power, sensitivity, and specificity were computed from the confusion matrix to quantitatively evaluate the estimation performance. The efficiency means the proportion of correctly predicted class and is most frequently used in the literature [54]. The efficiency of the two-phase flow model (80.19%) was greater than that of the single-phase flow model (74.79%). The positive predictive power means the proportion of TP in a total of positive predicted class. Low values of TP resulted from both the two-and single-phase flow models (0.23% and 0.24%, respectively). This is because a number of actual slope failure sites (161 cells) were quite small compared to the mountainous area considered in this study (63,157 cells). Even if all the actual slope failure sites are correctly predicted, the maximum value of TP is 0.25% (161/63,157). Therefore, regardless of the method of analysis, the value of the positive predictive power is consequently low. The positive predictive power of the two-phase flow model (1.17%) was slightly larger than that of the single-phase flow model (0.93%) due to the large FP, resulting in an over-prediction of landslides. Only for the sensitivity that is the same as the SR used to compute the MSR, the result from the single-phase flow model (93.17%) was better than that from the two-phase flow model (91.93%). As for the specificity, which is the proportion of negative observations accurately predicted, the two-phase flow model resulted in a greater value (80.16%) than the single-phase flow model (74.74%). Given the performance values produced by the MSR and confusion matrix, the two-phase flow model exhibits good agreement with the actual landslide events at Umyeon Mountain.

Conclusions
Although air flow changes the infiltration rate of rainfall and variations in the safety factor of unsaturated soils, it has rarely been considered. We applied the two-phase (water and air) flow model to infiltration and slope stability analysis and assessed the model applicability for the Umyeon Mountain landslides that occurred in 2011. We used the available historical rainfall data, topographic maps, and geotechnical/hydrological properties to construct the model. Location information of the actual slope failures, the MSR, and the confusion matrix were used to evaluate the model performance.
The following are interpretations of the results of the infiltration and slope stability analyses:

1.
Considering that air flow changes the rate of increase in P w on slopes with a low infiltration capacity when ponding occurs during heavy rainfall, and that the saturated hydraulic conductivity is relatively marginal in the mountainous areas of Korea (i.e., less than 2.78 × 10 −5 m/s according to the National Disaster Management Institute [57]), it is necessary to apply the two-phase flow model to accurately interpret rainfall infiltration.

2.
The initial safety factor prior to rainfall depends on the shear strength parameters (i.e., internal friction angle and cohesion); however, variations in the safety factor are strongly dependent on an increase in the P w rate. Slopes in the single-phase flow model rapidly become saturated during heavy rainfall because of the high rainfall infiltration rates without the interruption of air. Thus, safety factors decrease rapidly compared to those in the two-phase flow model.

3.
Landslide susceptibility maps change depending on an increase in P w , which depends on relative permeability. It is also sensitive to the model type. The MSR and confusion matrix yield the highest performance for the two-phase flow model with an appropriate range of stable cell coverage for the best simulation. Thus, it is concluded that infiltration and slope stability analyses using the two-phase flow model have good applicability in evaluating landslide events in Umyeon Mountain. However, it is necessary to additionally evaluate the applicability of this model for other landslide cases.

4.
We performed infiltration and slope stability analyses focusing on the geotechnical characteristics of unsaturated soils in the upper layer. Besides, we excluded geological characteristics given that most of the slope failure in the study area occurred at the colluvium layer with shallow depths up to 2 m. The two-phase flow model can be usefully applied to the region where shallow slope failure occurs primarily. However, further study is required to examine the effects of the geological structure for improving the applicability of the two-phase flow model to the region where deep slope failure occurs primarily, which is affected not only by geotechnical characteristics but also by geological structure. 5.
The performance of slope stability assessments at a regional scale is greatly affected by the uncertainty and variability of geotechnical and hydrological input parameters when physically based models are applied. Geotechnical and hydrological properties have probabilistically or statistically been characterized based on types of soil and lithology to properly consider the uncertainty and variability of them [58,59], whereas we used representative constant soil properties at each of the five zones as averaged from field measurements. The applicability of the two-phase flow model to a regional scale can be improved by further study applying the probabilistic approach for characterizing the uncertainty and variability of geotechnical and hydrological input parameters.