Abstract
In this study, the role of water infiltration on a rainfall-induced landslide in mudstone is investigated. Finite element analysis (FEA) and discrete element analysis (DEA) were employed to explore the driving mechanism in the prefailure regime and the dynamic runout process in the postfailure regime, respectively. The driving mechanism was revealed in terms of the pore water pressures, saturation, and displacement of the sliding zone. To account for water infiltration associated with the dynamic runout process of the landsliding, the sliding friction coefficient in DEA was examined. Based on the results, satisfactory agreement between the numerical analysis and landslide behavior was realized. In prefailure regime, the onset of landslide initiation was found through assessment of rapid change of source displacement (RCSD). The estimated seepage force was about 0.5 N at landslide initiation. In postfailure regime, the results demonstrate that water infiltration and transition in steepness play significant roles in the behavior of the dynamic runout process. The landsliding exhibited a maximum speed of 4.41 m/s and decelerated as it reached a more gentle slope. Overall, the study indicated that the coupling of FEA and DEA can be used to investigate the role of water infiltration and provide useful insights. The approach can be further applied for studies on many other rainfall-induced geohazards to disclose the role of water infiltration with the prefailure and postfailure characteristics.
1. Introduction
Rainfall-induced landslides are a common geohazard worldwide in tropical and subtropical regions [1–3], and studies have clearly indicated that these hazards are destructive to properties and lives [4–6]. Situated within a mountainous subtropical region with frequent typhoons, rainfall-induced landslides have been a critical issue for Taiwan [7]. To facilitate accurate prediction and provide countermeasures for reducing landslide damage, failure mechanisms of rainfall-induced landslides have been studied in recent years [8–10]. These studies indicate that water infiltration may play a critical role in landslide occurrence through a decrease in suction and increase in saturation and unit weight, making slopes unstable. Furthermore, to predict areas prone to rainfall-induced landslides, the dynamic runout behaviors of rainfall-induced landslides have been investigated [11, 12]. Although the above studies highlighted the role of water in rainfall-induced landslides, they are rather limited to areas where sandstone and shale prevail. Note that it is also essential that a practical approach towards connecting failure mechanism and runout behavior be established to enhance understanding of landsliding.
On September 28, 2016, typhoon Megi slammed Taiwan and forced more than 14,000 people to evacuate from mountainous regions near the coast to government shelters. Among the damage caused, a mudstone landslide in Yanchao district, Kaohsiung, was induced by the rainfall of typhoon Megi. Of particular note is that, although the affected area is considered small-scale (<10 ha) [13], it damaged a two-story residential building, led to three deaths, and raised public attention to landslides in the region. Despite the fact that mudstone is common in southwestern Taiwan and has a proclivity to fail under heavy rain [14, 15], the failure mechanism and dynamic process in these site conditions have not been extensively studied. Consequently, the 2016 Yanchao mudstone landslide provides a unique opportunity to examine the driving mechanism and dynamic process of these slope failures.
With rapid advancement in computing capabilities, numerical simulations have become a more prevalent tool to reconstruct and analyze the progress of slope failures, which facilitates better understanding of the driving mechanism and dynamic process of landsliding [16–18]. Of particular utility for the analysis of landslides is finite element analysis (FEA) and discrete element analysis (DEA). FEA has been utilized in investigating the role of water infiltration such as the change of pore water pressure and the strain of geomaterial interacting water. Leshchinsky et al. [19] used coupled transient unsaturated seepage-stress FEA to investigate the complex behaviors of the progressive slope failure in Yumokjeong, Korea. Yang et al. [10] coupled FEA and unsaturated seepage concept to investigate the failure and deformation mechanisms of slopes in Taipei, Taiwan. Both studies demonstrated that FEA is a useful approach for reviewing the role of water infiltration in prefailure regimes. Recent studies have also suggested the FEA can be used to the initiation time of landsliding and defines the initiation time as the time when a rapid change of source displacement (RCSD) surfaced [20, 21].
In addition to FEA, DEA is an effective approach for studying dynamic process of landslides that deaggregate or flow after failure, given that it treats geomaterials as discontinuous bodies [22–26]. The approach was first utilized to simulate problems relating to rock mechanics [22] and other geomaterials [23] and has been applied to simulate landslide problems in recent decades, as it is capable of simulating the kinematic process and runout dynamics that occur after failure [24–26].
In this study, a coupled approach of finite element analysis (FEA) and discrete element analysis (DEA) is employed to analyze the rainfall-induced landslide in Yanchao mudstone, a geomaterial common in the mountainous areas of southwestern Taiwan. Analysis of the landslide was divided into the prefailure and postfailure regimes, involving evaluation of the driving mechanism and dynamic runout process. The influence of partial saturation was considered through coupling of the soil water characteristic curve (SWCC), hydraulic conductivity function (HCF), and unsaturated shear strength for the modeled geomaterials [10, 19, 27–29], used to investigate the onset of failure. The variations of pore water pressure, saturation, and rate of source displacement in the sliding zone are compared and discussed. The dynamic runout process was simulated based on the commercial particle flow code [30], obeying Newton’s second law. The seepage forces and reduced sliding friction coefficient were considered within the discrete element method analyses to capture the impacts of water infiltration [11, 12]. The runout behavior of the landslide is characterized by discrete particles and walls, and six phases were presented to enable better understanding of the dynamic process of landslide movements after failure. This study provides new insight into the influence of precipitation on the onset of failure and the runout behavior of landslides that occur in mudstone using coupled numerical modeling.
2. Study Area
The landslide occurred in the Yanchao district of Kaohsiung city, Taiwan, in the early morning of September 28, 2016. The failure occurred as a result of rains brought by the typhoon Megi. The nearest rain station is the Agongdian rain station which is located approximately 3 km northwest of the landslide (Figure 1, 120.3552 long. 22.8041 lat., 56 m asl.), where rainfall totals were measured hourly. The rain data showed that the total accumulated rainfall was 408.5 mm from the period of September 27, 12 a.m., to September 28, 12 p.m. According to the Central Weather Bureau (CWB), the event was classified as an extreme torrential rain event. Interestingly, very few landslides had been observed in the areas comprised of mudstone, resulting in a newfound awareness of a hazard that was not considered previously. This event resulted in the revision of local landslide hazard maps, produced by the government, to include areas of mudstone that were potentially unstable.
In southwestern Taiwan, the mudstone formation was formed in the end of the Tertiary Period to the early Quaternary Period [15]. According to ISRM classification [31], mudstone is classified as very weak to extremely weak rock in a high water content. Based on its characteristic, the behavior of a mudstone landslide would progress like a flow-like landslide. Understanding the interplay of rainfall in a mudstone landsliding would help reveal the prefailure mechanism. Figure 2 shows the photos taken from the site investigation after the landslide, where cracks can be seen on the surface of the surrounding geomaterials. Figure 3 shows the cross section (line A-A′) before and after the landslide, composed of the mudstone (Ku-ting-keng formation) and weathered formations [32]. The resolution of the cross section was 16 cm and was created using a digital elevation model (DEM) with the aid of unmanned aerial vehicle (UAV) photos, in which the source area was depicted. Note that steeper and gentler slope angles of 27 and 15 degrees, respectively, are delineated based on the DEM. The sliding mass occupied an area of about 4350 m2, and the total estimated volume was about 8700 m3. The height, length, and width of the sliding zone (source area) are approximately 45, 90, and 30 m, respectively. Of particular interest is the observation that the damaged building at its toe displaced about 12–14 m from its original position. To acquire the permeability of geomaterials, the field test proposed by Hoek and Bray [33] was conducted and a saturated hydraulic conductivity of was obtained. Other laboratory tests, including grain size distributions, unit weight tests, and specific gravity tests, were also conducted to get the basic properties of sliding mass.
3. Methodology
The principles that govern FEA and DEA are briefly introduced to provide context to the numerical modeling performed in this study. To address the role of water infiltration, FEA incorporated the flow equation controlling the unsaturated filtration and the coupled stress-pore pressure formulation, which requires the integration of the concepts including soil water characteristic curve (SWCC), hydraulic conductivity function (HCF), and unsaturated shear strength [10, 19, 27–29]. The algorithm in DEA is based on particle flow code (PFC), seepage force, and sliding friction coefficients [11, 12, 26, 30, 34].
3.1. Finite Element Analysis (FEA)
In FEA, the failure mass of the landslide can be represented by a large number of small elements, enabling quantitative assessment of the prefailure behavior of the landslide, evaluated through time-dependent iterations of the coupled stress-pore pressure formulation until a convergence is achieved.
The movement of the water in unsaturated soil was derived from Richards equation [35]: where and represent the hydraulic conductivities in the and directions, respectively, is the total hydraulic head of flow, is water content, is the coefficient of water volume change, is the unit weight of water, is gravity, and t is the time.
The coupled stress-pore pressure formulation is defined as where is volume, is true stress, is the virtual rate of deformation, is the virtual velocity field, is surface tractions per unit area, is the body force barring the weight of the liquid (water), is porosity, is the volume of trapped wetting liquid per unit of current volume, is the degree of saturation, is the unit weight of water, and is gravity [19, 29]. The last term in (1) represents the body forces of trapped wetting liquid in the voids and the volumetric flow rate of water to capture both the weight of the water and the inertial effects of flow, known as the Forchheimer term [36]. The principle of virtual work represented in (1) is evaluated for equilibrium in the volume under consideration in its current configuration at time under transient seepage induced by boundary conditions representative of rainfall.
The SWCC describes the relationship between matric suction and the normalized volumetric water content while the HCF evaluates changes in hydraulic conductivity with the normalized volumetric water content. SWCC and HCF were utilized in the study following van Genuchten-Mualem’s model [27, 28]. Then, since the normalized degree of saturation is equal to the normalized volumetric water content, the SWCC and the HCF can be expressed as where is the normalized degree of saturation; is the water content; is the saturated water content; is the residual water content; , , and are constants; is suction (where and are the pore pressure and pore water pressure, respectively); is the relative hydraulic conductivity; is the unsaturated hydraulic conductivity; and is the saturated hydraulic conductivity.
To evaluate the unsaturated shear strength, the study uses the equation proposed by Vanapalli et al. [37], which is defined as where is the soil shear strength, is the normal stress on the failure surface at failure, is suction on the failure surface at failure, is the effective friction angle, and is effective cohesion.
3.2. Discrete Element Analysis (DEA)
The commercial program particle flow code (PFC), based on the concept of DEA, was utilized in this study. The program, featuring particles and boundary conditions, is suitable for analyzing the postfailure kinematic process of landslide runout [24, 34]. Through a continuous calculation cycle, the force-displacement law is obeyed on the contacts between particles and Newton’s second law is maintained to compute the new velocity and position of particles.
The effect of water infiltration on landsliding has been addressed, in which the importance of water head difference (seepage force) and the decreased interparticle sliding friction coefficient are determined [11, 12, 26]. To address the effect of water infiltration and realistically simulate the landslide, the study considered seepage forces and the decreased friction coefficient stemming from precipitation. The seepage force is computed by the following equation: where is the seepage force, is the unit weight of water, is the hydraulic gradient, and is the volume of each particle [11, 38]. Due to the fact that the damaged building (the reference for runout behavior) moved about 12–14 m in the Yanchao mudstone landsliding, the decreased friction coefficient could be explored through sensitivity analysis to provide insights on the effect of water infiltration [39].
4. Numerical Simulation
The setting of FEA and DEA is introduced in this section. The setting in the FEA includes the modeling, boundary conditions, and material properties, while the setting in DEA includes transformation and calibration of parameters and modeling.
4.1. FEA Model and Boundary Conditions
In FEA modeling, the landslide mass was discretized differently than the area outside the failure plane. The bottom boundary was restrained in the vertical and horizontal directions, and the two lateral boundaries were restrained in the horizontal direction. The two lateral hydraulic boundaries were set to permit change in the ground water level (GWL) accounting for water infiltration while the bottom hydraulic boundary was set as impervious. The slope geometry and boundary conditions are illustrated in Figure 4.
To calculate the change of the pore water pressure and saturation, eight-node pore-stress, reduced-integration quadrilateral elements (CPE8RP) were used to discretize the model. The model is comprised of 1933 elements, as shown in Figure 5. Thirteen monitoring points were prescribed in the sliding zone to measure the change of the pore water pressure, saturation, and displacements throughout the landslide mass and were divided into three sections including the toe, middle, and top (Figure 4).
The FEA model was divided into three steps including the gravity, antecedent rainfall, and typhoon rainfall. In the gravity step, the initial pore stress, based on unit weight of geomaterial, was computed and the distribution of the initial water pressure was established by hydraulic boundaries. The distribution of the initial water pressure was determined by achieving hydraulic equilibrium based on the development feasibility report by China steel [40]. In general, it is difficult to discern groundwater flow regimes with any great certainty in the field. However, at the assumed domain, the lateral and bottom boundaries were considered to be under “far-field” constant flow conditions (i.e., zero fluid flux), and the regimes were specified at nodes located along these boundaries.
The effect of the antecedent rainfall was considered in the study because previous studies have revealed that antecedent rainfalls have an influence on slope stability [18, 41]. According to Rahardjo et al. [42], cumulative rainfall in the week prior to a landslide event is important in considering susceptibility of slope failures. Recorded rainfall data for the week preceding the Yanchao mudstone landslide is considered for the purpose of establishing realistic initial hydrologic conditions in the study, as shown in Figure 6(a). Finally, the typhoon Megi rainfall data obtained from the Agongdian rain station was applied in the typhoon rainfall step (Figure 6(b)) to incite a slope failure.
(a)
(b)
4.2. FEA Geomaterial Properties
The geomaterial parameters of the sliding zone and nonsliding zone were established from laboratory tests [15, 43, 44] and were characterized by an elastoplastic Mohr-Coulomb model widely used in the slope analysis [19, 45, 46]. Table 1 summarizes the input parameters of geomaterials used in the present study. The unsaturated hydraulic properties including SWCC and HCF were evaluated by the RETC program combining the grain size distribution and van Genuchten-Mualem’s model [27, 28], as shown in Figure 7.
(a)
(b)
(c)
5. Transformation and Calibration of Parameters
In DEA, microparameters such as stiffness and parallel bond stiffness were characterized to represent interactions between particles; however, microparameters cannot be obtained from physical experimentation. Through transformation and calibration, the microparameters can be determined. The transformation procedure used in this study was suggested by Potyondy and Cundall [47]. The microparameters were calibrated through uniaxial compression test to reproduce macroparameters such as Young’s modulus, strain behaviors, and uniaxial compressive strength [24, 25, 48]. The calibrated microparameters and macroparameters are presented in Tables 2 and 3, respectively.
5.1. DEA Model
In the DEA model, the profile of the Yanchao mudstone landslide was created for subsequent release. When simulating runout problems, the stipulated distribution of particle sizes should capture a balance between model resolution and practical computational time [49]. When the particle sizes are adequately small, it sufficiently simulates the flow-like landslide behavior after failure, enabling movement of geomaterials in both vertical and longitudinal directions [26]. To capture this behavior, the sliding mass was constructed using 16,200 particles with radiuses of 0.08–0.09 m, and the nonsliding zone was represented by no displacement boundary condition.
To investigate the effect to water infiltration on landslide runout, a sensitivity analysis was performed against the movement of a building displaced by the landslide at the base of the slope, discretized by rigid particles based on Design Specifications for Concrete Structures Taiwan [50], and the parameters are shown in Table 3. The seepage force could be estimated by coupling the water head difference obtained from the result of FEA and (5).
6. Result and Discussion
6.1. FEA Results and Discussion
Figures 8, 9, and 10 illustrate the variation of saturation and pore water pressure at the monitoring points in three sections during typhoon Megi. As shown in Figures 8, 9, and 10, when the rainfall intensity is about 0.5 mm/hr at the initial period of precipitation duration (10–12 hr), the variations of saturation and pore water pressure were almost constant in all monitoring points. At , there is a notable increase in saturation and pore water pressures near the surface. This observation is reasonable as rainfall must first pass through the monitoring points near the surface (mp1, mp6, and mp10), called a wetting front. With continuous water infiltration, saturations and pore water pressures increased gradually. At , the rainfall intensity increased to 48 mm/hr, creating saturations near 100% and positive pore water pressures at the monitoring points near the surface. At , when the water infiltration caused the perceivable increase in the ground water level, pore water pressure near the ground water level (mp3) turned positive. At , the negative pore water pressure (suction) almost dissipates in the toe section and decreased the shear strength and destabilized the slope. In addition, Figure 11 illustrates the displacement for monitoring points near the failure (mp5, mp9, and mp13) of the three sections and the FEA contours during the typhoon rainfall. It is observed that rapid change of source displacement (RCSD) occurred at 28 hr. These results indicate that the likely failure timing was at approximately 28 hr from the beginning of the typhoon rainfall. The results observed in the three monitoring sections also indicated that the effect of the water infiltration resulted in variation of pore water pressure, saturation, and ground water level, which was consistent with the previous studies [51, 52].
(a)
(b)
(a)
(b)
(a)
(b)
(a) Precipitation duration (hr)
(b) Pore water pressure
(c) Deformation
(d) Stress
Figure 12 shows the variation of pore water pressure at the monitoring points near the failure of the three sections during the typhoon rainfall. In Figure 12, the response timing of the pore water pressure at mp5, approximately 17 hr, was faster than that at mp9 and mp13, which were approximately 20 hr and 23 hr, respectively. Such an observation is reasonable given the monitoring point, as higher permeability causes a shorter response timing, corroborating the lag observed in saturation levels at different monitoring points, each with varying conductivities. The variation of the pore water pressure at mp5 is small because it is located under the ground water level (GWL). The head difference between two hydraulic boundaries at was 9 m. Through coupling the water head difference obtained from the result of FEA and (5), the seepage force was determined to be 0.5 N.
6.2. DEA Results and Discussion
The postfailure behavior of the Yanchao mudstone landslide is considered through DEA modeling that considers the seepage force and reduced friction stemming from precipitation. Through FEA, the failure timing of the landslide was determined to occur at 28 hr after the typhoon started, whereas the seepage force at this time was 0.5 N. The seepage force determined from this analysis was applied in DEA to account for more realistic initial conditions on landslide runout behavior. Although the friction coefficient required in the DEA cannot be obtained from the physical experimentation, it could be determined through iterative study and comparison of the displacement of the reference building (Figure 13).
Figure 13 shows the relationship between the friction coefficient and the displacement of the building obtained through DEA. A displacement of 13 m was observed with a friction coefficient of 0.05, agreeing to the result observed in the field. It is revealed that a minor difference of the friction coefficient would lead to a significant difference in the displacement of the reference building, as shown in Figure 13. The result also suggested that water infiltration plays a critical role in the kinematic process of the landslide runout as well as the affected areas.
To better understand the kinematic process of the landslide travel, the failure process is divided into six phases (Figure 14). The six phases represented (a) the initial phase, (b) acceleration phase, (c) blockage phase, (d) collision phase, (e) interaction phase, and (f) the final phase. In general, the results and discussion on each stage are as follows: (i)Stage (a): at the beginning of the postfailure regime, the seepage force was applied, coupling the result of FEA.(ii)Stage (b): the sliding mass accelerated along the steeper section of the slope (27°), indicating that the front of the sliding mass was dominated by free falls attributed to the terrain. The maximum average velocity of the sliding mass was reached at this stage, at approximately 4.41 m/s.(iii)Stage (c): the acceleration of the front of the sliding mass began to drop noticeably as it reached a gentler slope while the back of the sliding mass accelerated continuously. At this stage, the sliding mass had almost blocked the road.(iv)Stage (d): the sliding mass moved continuously and began to collide with the reference building. Due to the sliding mass colliding the reference building, the velocity decreased slightly once again.(v)Stage (e): the interaction between sliding mass and reference building was observed. In this stage, the kinematic process did not yet reach a force balance.(vi)Stage (f): in the final phase, deposition of the sliding mass was completed. The final displacement of the reference building was 13 m, satisfying the result observed in the field.
(a)
(b)
(c)
(d)
(e)
(f)
Figure 15 presents the relationship between average velocity and time obtained through DEA. As shown in Figure 15, the average velocity of the landslide before colliding with the reference building was 1.97 m/s ranging from very rapid to extremely rapid velocity classes according to Cruden and Varnes [53], leaving only 16 s to respond from the onset of the landslide. Besides, the velocity fluctuates as a result of particle collisions during the sliding process [54]. Based on the result, it is revealed that the variation of velocity was dominated by the variation of terrain.
6.3. General Discussion
The focus of this paper is on the application of a coupled finite and discrete element on analyzing a rainfall-induced mudstone landslide. The results show that the effects of water infiltration including increases of saturation, pore water pressure, and displacement could be observed during the typhoon rainfall. The variation of saturation and pore water pressure near the surface increased remarkably at , which was earlier than other locations, because rainfall first passes through the monitoring points near the surface.
In general, there are three types of rheology commonly used for simulations of landslides with rapid movements: (i) frictional rheology, (ii) the Voellmy rheology, and (iii) the Pouliquen rheology [55, 56]. Pirulli and Mangeney [56] conducted a comparison between the three rheologies by back analysis on a rock avalanche. They concluded that the three rheologies assumed were capable of producing the runout area and the final deposit with satisfactory accuracy. The main advantage of the frictional rheology is that only one parameter, the dynamic friction angle, needs to be calibrated. This advantage is particularly important for the sensitivity analysis performed in the study. For this reason, the frictional rheology was selected with the expectation that it will provide practical results.
Based on site investigation, the building displaced on its shallow foundation—no basement was present. It is observed that the present setting would reproduce the in situ condition specifically for the postfailure kinematic process. The study has assumed the building, with properties based on the Design Specifications for Concrete Structures Taiwan [50], standing on the ground without the basement or foundation. Additionally, based on the site investigation, the impact of the lateral dispersion of landslide debris at the toe was not considered critical as it was concentrated within a landslide width similar to that of the zone of rupture [57]; thus, 2D modeling was considered sufficient. By assessing the DEA results obtained in the study, the landslide progressed as a mudflow with an average velocity of 1.97 m/s, ranging from a class of very rapid to an extremely rapid [53]. The result was satisfactory as compared to field observations and case studies by Hungr et al. [58].
7. Conclusions
This study used a combination of FEA and DEA to investigate the process of Yanchao mudstone landslide including prefailure and postfailure regimes. The initiation time of landslide was investigated. The effect of water infiltration, the runout behaviors of mudstone landsliding, and general discussions regarding the numerical setting were also examined and discussed. The following conclusions can be drawn from the following results: (i)The initiation time of Yanchao mudstone landslide was revealed by the RCSD, and the likely failure timing was ascertained as 28 hr from the beginning of the typhoon rainfall(ii)The landslide hit the residential building within 16 s. During landslide runout, the primary increase in the sliding velocity was attributed to the steepness of the slope, suggesting that the runout behavior is dominated by the terrain(iii)The frictional rheology for simulating the kinematic process of the Yanchao landslide is reasonable, based on the DEA results. However, further investigations of the landslide kinematic features considering different rheology may also provide additional insights into the interplay of rainfall in the landslide(iv)The present study demonstrates that a combination of FEA and DEA can be an effective and useful approach to better understand the process of rainfall-induced landslides
The results of the study provided useful information about the role of water in the failure of mudstone slopes. This insight is valuable for formulating disaster prevention and mitigation strategies.
Data Availability
The data that support the findings of this study are available from the corresponding author on request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors gratefully acknowledge financial support from the Ministry of Science and Technology (MOST) in Taiwan: Grants 107-2636-E-006-003 and 106-2221-E-006-061. Special thanks go to the Young Scholar Fellowship Program by the Ministry of Science and Technology in Taiwan (The Pilot Directions for MOST Grant for the Columbus Program). The research, in part, received funding from the Headquarters of University Advancement at the National Cheng Kung University, which is sponsored by the Ministry of Education, Taiwan, ROC. The authors appreciate the colleagues in our laboratory for their constructive comments and assistance.