Abstract
Although the mechanisms of slope failure caused by rising groundwater have been widely investigated, the kinematic behavior of landslides in the postfailure stage, which contains essential information for hazard mitigation and risk assessment, has not yet been fully studied. Thus, in this study, a series of numerical simulations using the material point method (MPM) were conducted to analyze the kinematic behavior and soil movement of shallow landslides (infinite slope problems). First, the proposed MPM formulation was validated in a full-scale landslide flume test. The simulated results of final slope profile, runout distance, deposit height, shear band development, slope displacement, and velocity accorded with the experimental results, suggesting that the MPM can quantitatively simulate large deformations. A parametric study of shallow slopes with various hydrological conditions and soil hydraulic and soil mechanical parameters was then performed to assess the influence of the aforementioned factors on landslide kinematics. The simulation results indicated that mechanical behavior at the slope toe is complex; the multiple plastic shear bands generated at the slope toe were due to a combination of shearing and compression. The deposition profile of the slopes was significantly influenced by all input parameters. Among the aforementioned parameters, soil cohesion, location of the groundwater table, and saturated soil permeability most greatly affected runout distance in the sensitivity assessment. Soil friction angle had a minor influence on the kinematic behavior of the slope.
1. Introduction
Landslides are major geotechnical disasters that occur worldwide. Forecasting landslide kinematics, such as the deposition profile, travel distance, and velocity of the unstable mass, is key in preventing and controlling the risks that landslides pose to infrastructure and human lives [1]. Many researchers have developed both physical and numerical models as well as conducted field observations to understand the kinematic characteristics of landslides [2]. However, most studies have focused on evaluating the failure mechanism of slopes rather than the postfailure behaviors of landslides [3]. Understanding the kinematic behavior of landslides remains a challenge that has yet to be completely surmounted.
Numerous experiments using small-scale [4–8] and large-scale [9] landslide experiments have been conducted to investigate debris flows behavior. Although existing methods for identifying the failure mechanisms of slopes yield undeniable advantages, most of these methods have mainly focused on the behavior of noncohesive soil [10]. Moreover, experimental approaches are time-intensive, and a scale model, unlike a field observation, cannot capture the same initial conditions.
With the development of increasingly powerful computers, numerical simulations have been increasingly used as powerful tools for slope stability analysis [11]. However, the numerical modeling of landslide kinematics (particularly deformation and runout) remains challenging for most geotechnical engineers [12]. Traditional numerical methods in stability analysis, such as the limit equilibrium method (LEM) [13] and finite element method (FEM) using the standard Lagrangian approach [14], remain useful only for analyzing slope response in the early stages of failure (i.e., the prefailure and failure stages), especially for identifying the location of failure surfaces and obtaining accurate safety factors, when the small strain assumption is considered reasonable [15]. These methods generally have limitations in problems associated with the development of large deformations after failure (i.e., the postfailure stage) because of excessive mesh distortions in the shear zone that lead to either inaccurate results or a calculation failure from a failure to converge [16]. In addition, advanced techniques that use the updated Lagrangian FEM formulation—in which first, the computational mesh is assumed to deform with the considered body and second, all static and kinematic variables that refer to the last calculated configuration in the solution require continual remeshing—lead to unstable and inaccurate results for large strain problems [17]. Troncone [18] simulated the failure mechanism of the Senise landslide at Southern Italy, and their simulated location of the failure surface, as identified in their finite element (FE) result, was consistent with the observed location; however, the predicted crest deformation was limited to a maximum of 0.53 m at failure onset. This predicted crest deformation was smaller than that in field observations. This small crest deformation was obtained even when a Lagrangian finite element formulation, updated for a large deformation, was used by Mohammadi and Taiebat [3] in a later study.
To overcome the drawbacks of mesh-based methods, a robust method for addressing large-deformation problems is warranted. Particle-based methods (i.e., meshless or mesh-free methods)—which include smoothed particle hydrodynamics (SPH) [19], the discrete element method [20] proposed by Cundall and Strack [21], and discontinuous deformation analysis—are robust approaches to solving large-deformation problems and to describing the postfailure mechanisms of landslides. These methods do not require data to be stored or the definition of a mesh to feature the connection of points. However, these methods incur a relatively high computational cost because neighboring particles must be intensively searched after each time step. This results in consistency loss [22]. In addition, SPH is limited by its inability to simulate multiphase interactions that involve failure evolution. Soga et al. [23] summarized the advantages and disadvantages of numerical methods for simulating landslide problems. By using both Eulerian (i.e., a fixed finite element grid) and Lagrangian (i.e., moving material points) formulations, Sulsky et al. [24] developed the material point method (MPM) to simulate the dynamic deformation of solid bodies. MPM has been successfully used to simulate geotechnical problems related to the large deformation of slopes and dams [1, 23, 25–31]. Although the MPM effectively models slope failure, few comprehensive overviews have been conducted on the influence of hydrological condition, soil hydraulic, and soil mechanical parameters on the postfailure behavior and kinematics of landslides.
The aforementioned discussion motivated the authors to conduct a series of MPM analyses to investigate the postfailure process, failure mechanism, and kinematic behavior of shallow landslides related to infinite slope problems. The sensitivity of the hydrological condition, soil hydraulic, and soil mechanical parameters in the kinematics of shallow slopes is evaluated and quantified. This remaining paper is organized as follows: first, an introduction of kinematic MPM formulations for a two-phase single-point MPM formulation used in the Anura3D MPM Research Community [32] is presented; subsequently, a validation of the MPM using the full-scale landslide flume test performed by Moriwaki et al. [33] is described; finally, a parametric study investigating the parameters influencing the postfailure kinematics of shallow slopes and a sensitivity assessment is detailed.
2. Methods and Validation
2.1. Basis of MPM Formulation
In this study, a two-phase single-point MPM formulation [34, 35] was used to represent the saturated porous medium, where the governing equations of the mass and momentum conservations of both the solid and liquid phases were written by Fern et al. [36] as follows: where is the material time derivative, is the material density, is time, is the velocity vector, and is the gradient operator.
The mass conservation for saturated soil is also known as the storage equation and can be expressed as where and are the velocities of the solid and liquid phases, respectively; is the porosity of the solid skeleton; and is the volumetric strain of the liquid phase.
The momentum conservation for the mixture and liquid phase (per unit of liquid volume) is expressed in Eqs. (3) and (4), respectively. where and are the solid and liquid densities, respectively; is stress tensor; is the volumetric concentration ratio of the liquid; is liquid pressure; is the viscosity of the liquid; is the intrinsic permeability of the liquid; is the gravitational acceleration vector; and is the drag force, which represents the interaction force between the solid and liquid phases. Based on the assumption that flow is considered laminar and stationary in the slow velocity regime (i.e., linear), drag force can be defined using Eq. (5) with Darcy’s law.
Figure 1 illustrates a computational cycle in the MPM. To solve the aforementioned momentum balance equation, the continuum body must first be spatially discretized by mapping information from the material points to the computational nodes of the mesh (referred to as the background computational mesh) (Figure 1(a)). At the beginning of each time increment, the momentum equations are solved on the predefined background mesh (i.e., the nodal accelerations) (Figure 1(b)). These nodal values are then used to update the acceleration, velocity, and position of the material points within a newly generated background mesh (Figure 1(c)). Because no permanent information is stored in the mesh, it can be freely redefined at the end of the time step. Finally, the assignment of material points to elements is updated after the background mesh is reset for the next time step (Figure 1(d)). The governing equation, time discretization, and solution procedure of MPM were detailed by Fern et al. [36].
2.2. Model Validation
To validate the proposed numerical model, the MPM result was compared against the results of a full-scale landslide flume test (Figure 2) conducted by Moriwaki et al. [33]. The real landslide case typically involves certain uncertainties associated with in situ conditions such as slope profile, hydrology, and subsurface ground conditions. Compared to real landslide cases, the reasons for selecting the experiment test reported by Moriwaki et al. [33] are because the selected experimental test was conducted in an indoor, well-controlled environment, resulting in better and more accurate test results for the use of the model validation. The model slope in the landslide flume test was 23 m long, 7.8 m high, 3 m wide, and 1.2 m deep, and the slope comprised three parts: an upper 30° slope segment, a lower 10° slope segment, and a horizontal segment at the toe of the model slope. The soil used in the test was loose Sakuragawa River sand, which is classified as a poorly graded sand (SP) according to the Unified Soil Classification System. A constant rainfall intensity of 100 mm/h was applied to the slope. A rapid landslide occurred on the 30° slope section 9267 s after the rainfall started. The collapsed soil mass slid downward and reached a final position within approximately 5 s from 9267 to 9272 s (postfailure stage). The test setup and experimental results were detailed by Moriwaki et al. [33].
Table 1 presents the mechanical properties of the input soil used in the MPM simulation. The soil properties of the Sakuragawa River sand were selected on the basis of those reported by Moriwaki et al. [33], Ghasemi et al. [37], and Yang et al. [38]. Sand behavior was analyzed using the elastic–perfectly plastic model. A 0.1 m thick steel flume was modeled as an impervious bedrock layer with a relatively high Young’s modulus to ensure that it behaved as a rigid material.
Moriwaki et al. [33] observed the critical water table immediately before slope failure, and its modeling is presented in Figure 2. Therefore, in the initial geometry in the MPM, the slope configuration comprised two layers under the assumptions that, first, the dry material (one phase, single point) is part of the upper layer above the observed phreatic surface at the time of failure and, second, the lower layer below the phreatic surface was a fully saturated material and was fully coupled (i.e., two-phase and single point). Slope deformation that occurred during the prefailure and failure stages was ignored. Figure 2 presents the initial geometry and discretization in the model validation. To obtain accurate results, the mesh was refined for sand material (i.e., element size of 0.1 m) because failure was expected to occur on the entire slope. The material point sizes were selected based on the relationship between the model height and the material point sizes for different cases of slope-related problems suggested by Llano-Serna et al. [1]. In total, 7567 elements and 3943 nodes with six material points per grid element were used in the model. Regarding boundary conditions, the bottom boundary was fully fixed, whereas roller boundaries were prescribed at the left and right sides of the model. The lateral and bottom models were impervious to the liquid phase. Because there is no information of interface friction coefficient reported by Moriwaki et al. [33], the frictional contact between the soil layer and bedrock was simply assumed to be full contact in the simulations.
The numerical simulation comprised two steps: (1) gravity loading to generate in situ stresses and (2) simulation of the 5 s propagation stage of the landslide. Because the Anura3D framework uses an explicit kinematic formulation to solve the governing equations [34], a time step of 0.1 s over 50 steps was used in the simulation to satisfy the requirements in the results. In the kinematic process, a local damping factor of 0.05 was used to simulate the slope failure, as suggested by Yerro et al. [30]. Furthermore, as reported by Abe et al. [39], during slope failure, an increase in pore pressure within the slope is synonymous with the change in soil permeability and stiffness. Therefore, a calibration parameter for both soil stiffness and permeability was used in this study to match the observed runout distance and final configuration in the experiment.
Figure 3 shows the landslide process in terms of the development of the deviatoric shear strain and deformation contours at different time points. In general, the observed and predicted final slope configurations were consistent. Both the landslide flume test and numerical simulation indicated that failure occurs in two stages. At the first stage of failure initiation (Figure 3(a)), the shear band is mainly formed along the soil–bedrock interface in the upper slope, and the sliding failure is a translational slide. In the second stage, failures occur within the moving mass when the shear zone is fully extended from the end of the upper slope to the end of the lower 10° of the slope (Figures 3(b)–3(c)). The shear band appearing in the lower part of the 10° slope was more complicated than that of the upper part because of the combination of compaction and shearing. As presented in Figure 3(h), the simulation results were consistent with the runout distance measured in the experimental test. At , the computed maximum depositional height was located near the toe of the 10° slope and was consistent with the experimental observations (Figure 4(a)). Figure 4 compares the simulated and observed shear band distributions in the final configuration. The kinematics of shear band development were similar to those observed in the MPM simulation.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(a)
(b)
(c)
Figure 5 shows the comparison between measured and predicted surface displacement and velocity for four monitored points on the slope surface. As illustrated in Figure 5, the predicted and measured surface displacement and velocity are close for the monitored points at the lower slope (i.e., D-7 and D-8), but the numerical results generally overly predicted both surface displacement and velocity at the upper slope (i.e., D-3 and D-6). The discrepancy is primarily attributable to the change of soil parameters (i.e., soil modulus and permeability) under different loading conditions during the landslide process. As the landslide progresses, the soil at the upper slope is mainly under extension, while the soil at the lower slope is subject to compression. In reality, soils under different loading conditions may have different void ratios, resulting in different soil mechanical and hydraulic parameters. However, this effect of change of soil parameters with soil void ratio (or loading conditions) is not modeled in this study due to the limitation of the current Anura3D MPM version. The above statement has been confirmed by conducting a trial simulation using different soil moduli and permeabilities for the soils at the upper and lower parts of the slope. The numerical results show the predicted and measured surface displacement and velocity match well for both the upper and lower parts of the slope. However, manually inputting two different sets of soil properties for soils at the upper and lower part of the slope is not a rigorous numerical procedure, and, thus, the results are not presented in this paper. For a more rigorous numerical procedure, advanced MPM functions should be implemented to resolve this limitation; for example, considering two phases, two points function to allow the soil properties to vary with the soil void ratio.
(a)
(b)
3. Postfailure Behavior of Shallow Landslides
3.1. Numerical Model and Inputs
A shallow slope with a 3 m thick residual soil layer on top of the impermeable rock layer, a length of 90 m, and an incline of 30° (Figure 6) was adopted in this study to evaluate the influence of soil parameters and groundwater level on the postfailure characterization and kinematic behavior of shallow landslides. The slope model and input soil properties were designed according to the statistical data reported in a companion paper by Yang et al. [38]. In this study, the soil properties, the thickness of the residual soil, and slope angle used in this study were statistically determined from a large soil database (i.e., 57 soil types) compiled from 35 landslide case histories of the literature (see Table 6 in Yang et al. [38]). The datasets cover a wide range of soil types, including residual weathered soils and transported colluvium deposits worldwide. For example, more than 80% of the landslide cases in the collected landslide database have soil thickness ranging from 1 to 5 m. Based on the above statistical results, a 3 m thick soil layer, the average value of 1-5 m, was selected in this study.
For the slope geometry and boundary conditions, the ratio of slope length to soil thickness in the slope model (L/H) was 30, which is larger than the suggested ratio of 25 to ensure no interference from boundaries in the calculation [40]. In addition, the length of the horizontal section at the toe of the model was extended by 100 m to ensure that the soil travel distance was unaffected. To save computation time, only a 1 m thick layer of bedrock was imposed in the model. Moreover, the water table was assumed to be located in the middle of the soil slope (Figure 6). The model comprised two parts. The upper layer of soil (i.e., that above the water table) was assumed to behave as a dry material, whereas the lower layer of soil was assumed to be fully saturated. The boundary conditions and numerical procedure were identical to those used in the validation model, which used a roller for the two side boundaries and a hinge for the bottom boundary.
The soil layer and bedrock were modeled under the assumption of the Mohr–Coulomb constitutive model with a nonassociated flow rule. During the analysis, the effective stress parameters under drained conditions were measured for the soil layer submerged in water by using the two-phase single-point MPM formulation. Figure 6 presents the geometric configuration and boundary conditions of the slope. The initial domain of the model was discretized using 8494 mixed triangular finite elements, 4453 mesh nodes, and 16902 material points in plane strain. The simulation lasted 25 s in total after landslide initiation.
The numerical analyses were performed in two series: a baseline case and a parametric study. For the baseline case, the input soil properties listed in Table 2 were obtained based on the mean soil property values from a previously compiled dataset of 35 landslide case histories in the literature [38]. A low soil cohesion value (i.e., ) was selected for the baseline case for a better observation of the postfailure behavior of landslide. For the parametric study, the effect of three groups of input parameters—namely, hydrological conditions (i.e., groundwater level), soil hydraulic parameters (i.e., saturated soil permeability), and soil mechanical parameters (i.e., soil friction angle, cohesion, and Young’s modulus)—on postfailure behavior such as travel distance and velocity was evaluated. In each case, only one parameter was manipulated, and all other parameters were constant in the baseline case. In total, 11 cases of numerical simulations were performed in this study (Table 3). The soil property values within one standard deviation (SD) of the mean () from the landslide database [38] were selected in the parametric study. If was negative or unreasonable, the minimum value of the parameter from the landslide database was used as the lower bound.
3.2. Baseline Case
Figure 7 presents the evolution of slope failure in the baseline case at , 5, 7.5, 10, and 25 s in terms of the deviatoric strain and displacement contours. At the start of the failure stage (Figure 7(a)), the shear strain developed evenly along the soil–bedrock interface on the slope section, and the sliding failure surface extended from the crest to the toe of the slope. The failure mode was categorized as a translational sliding failure. Over time, the landslide mass moved rapidly down the slope along the inclined surface to the toe and extended to the horizontal part of the slope as a complex failure surface. Displacement within a new sliding failure increased in each successive stage (Figure 7(e)–7(h)). Figure 7(d) presents the slope situation at , when the sliding failure on the slope was complete. Similar to the results observed in the model validation presented in Section 2.2, the distribution of shear strain in the horizontal part of the slope was affected by combined soil compression and shearing with the formation of multiple plastic shear bands. Figure 7(h) presents the final slope displacement. The runout distance, which was measured in terms of particle movement, was approximately , and the ratio of landslide height to travel distance (Hh/Ld) ratio was 0.4.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
To highlight the kinematics of the failure, Figure 8 displays a comparison of the displacement and velocity evolutions over time at material points T-1 and T-2 on the slope surface (Sections A and B in Figure 6(a)). In the numerical simulation, the soil mass on the 30° slope moved rapidly downward at a calculated maximum velocity of 10.2 m/s. Subsequently, (1) the kinetic energy of the slope decreased over time because shearing occurred between the flowing and deposited material at the horizontal section, and (2) the slope stabilized after (i.e., when slope failure was complete). The maximum thickness of the deposited material was 6.5 m in the final deposition profile (Figure 7(h)). Notably, the movement speed of material points on the 30° slope section was three times that of the material points in the horizontal section, and the magnitude of the runout and the maximum displacement of points on the slope surface were similar in this study.
(a)
(b)
3.3. Parametric Study
3.3.1. Effect of Water Table Location
To measure the effect of water table position, three positions on the phreatic surface (hereafter referred to as GWL) were adopted in the analysis, namely, GWL1, GWL2, and GWL3, which corresponded to the top, middle, and bottom of the soil layer, respectively (Figure 6(b)). The material model for the soil layer was assumed to be saturated-fully coupled and dry the corresponding water table was located at the surface slope (i.e., GWL1) and at the soil–bedrock interface (i.e., GWL3), respectively. Section A in the middle part of the slope and section B at 20 m away from the toe of the slope were selected to present the numerical results. In addition, the two quantities of runout distance (), which was defined as the distance between the end of the displaced material after failure and the toe of the initial slope, and the ratio of landslide height to travel distance (Hh/Ld) were used to evaluate the postfailure process in the present study.
Figure 9 depicts the failure evolution for both cases with contour plots of the deviatoric shear strain. The kinematic behavior of the slope after failure was significantly affected by water table location. In the GWL1 case, the deviatoric shear strain distribution in the horizontal slope section was significantly longer than that in the baseline case (i.e., GWL2). The final runout distance increased to approximately 72 m, and the corresponding Hh/Ld was 0.32. At , the maximum depositional height increased by more than 30% to 8.5 m relative to the baseline case. This likely occurred because (1) excess pore water pressure caused by the large distortions formed after slope failure and (2) the soil weight increased as the entire slope became saturated. Figure 9(b) illustrates the evolution of the shear strain contours for the GWL3 case. In contrast to the GWL1 case, the postfailure process in the slope did not occur when the phreatic surface level dropped to the base of the soil layer. This occurred because the initial slope conditions exhibited a minimum factor of safety (FS) larger than 1.0 (stable condition). Because the MPM used in this study cannot provide information related to slope stability, the initial FS value herein was calculated using the LEM. The obtained FS was 1.12 (Figure 9(d)).
(a)
(b)
(c)
(d)
A comparison of velocity at material points T-1 and T-2 between the GWL1 and GWL3 cases is presented in Figures 10(a) and 10(f). The GWL1 case exhibited a rapid development of velocity in the sliding particles on the slope. The velocity of point T-1 peaked at 12 m/s at 5 s (Figure 10(f)); this was much higher than the 0.15 m/s peak in the GWL3 case. Notably, unlike the baseline case, after reaching a peak, the velocity of point T-1 exhibited a slight discontinuity at when the material points moved to the toe of the slope. As illustrated in Figure 10(f), a new equilibrium of the slope was established after . As previously observed in the GWL3 case, because runout did not occur in this case, the displacement and velocity of point T-2 on the horizontal section of the slope was close to zero.
(a) Effect of GWL
(b) Effect of
(c) Effect of
(d) Effect of
(e) Effect of E
(f) Effect of GWL
(g) Effect of
(h) Effect of
(i) Effect of
(j) Effect of E
3.3.2. Effect of Soil Permeability
Regarding the effect of saturated hydraulic conductivity, as soil permeability increased, deformation decreased considerably (Figure 11). For the slope with high permeability (i.e., ), the failure plane in terms of deviatoric shear strain formed clearly on the 30° slope section (Figure 11(a)). In the final configuration (Figure 11(b)), the runout () was approximately 8 m. By contrast, the slope with low permeability exhibited substantially increased deformation with a maximum runout distance of approximately 47.5 m. This result occurred because, during kinematic postfailure, the generated pore water pressure increment in the slope with low permeability was larger than that in the slope with high permeability. The numerical results are consistent with the findings of Ghasemi et al. [37]. Notably, when , the final deposit profiles did not differ significantly between the simulated and baseline cases. In other words, when the permeability coefficient was less than a certain value, the final configuration was not affected by changes in . However, low soil permeability can significantly increase computation time, and this finding is consistent with that of Yerro et al. [29].
(a)
(b)
(c)
(d)
A comparison of the evolutions of deformation and velocity at different values are presented in Figures 10(b) and 10(g). Both velocity and deformation considerably decreased as increased. Compared with a slope with high permeability (i.e., ), the velocity at material point T-1 doubled when . In addition, as illustrated in Figure 10(g), the time to peak velocity for a slope with a high was 8 s, which was longer than the 5 s observed for the slope with a low . Overall, the numerical results indicated that permeability considerably influences the kinematic behavior and soil movement of slopes after the failure stage. Slopes with high values appeared more stable than those with low .
3.3.3. Effect of Soil Mechanical Parameters
(1) Effect of Soil Cohesion. Regarding the influence of soil shear strength, the kinematic behavior of slopes was strongly affected by soil cohesion and moderately affected by soil friction angle. Figure 12 presents the influence of soil cohesion values (i.e., and 7 kPa) on the final configurations of slope collapses. Generally, slopes with low effective soil cohesion had larger deformations. For the slope with noncohesive soil (), the behavior of the material on the slope section resembled that of a granular flow [41–45] with extremely large deformation. No distinct shear band was observed in the slope with noncohesive soil, and the time required to reach a new equilibrium condition () was longer than that of a slope with cohesive soil. As illustrated in Figure 12(b), the final runout distance for a slope with was 50 m, and the corresponding Hh/Ld was 0.37. Similar to the GWL3 case discussed earlier, for a slope with , the FS value of 1.05 obtained using the LEM confirmed that the initial state of the slope was stable, and the runout distance was not observed. The shear strain distribution (Figure 11(b)) was mainly observed on the 30° slope section along the soil–bedrock interface, and it did not extend to the horizontal part. Compared with the baseline case () and the previous case with , low cohesion resulted in a shallow slip failure surface, whereas higher cohesion resulted in a deeper sliding failure. Yerro et al. [30] and Shi et al. [22] also obtained similar findings.
(a)
(b)
(c)
(d)
The kinematics of the landslide were evaluated using the displacement and velocity evolution of the material points (i.e., T-1 and T-2), as presented in Figures 10(c) and 10(h). The velocity and displacement in the slope with were clearly more serious than those of other cases. For a slope with noncohesive soil, the speed of the mass movement on the slope was five times that of movement on the slope with high cohesion (Figure 10(h)). Notably, the time required to reach the peak value (i.e., ) was unaffected by the cohesion value. In addition, as observed in Figure 10(c), no displacement of material point T-2 occurred in the horizontal section of the slope when . In conclusion, soil cohesion substantially affected the postfailure landslide process.
(2) Effect of Soil Friction Angle. The effect of internal friction angle on runout postfailure characteristics was also investigated. Figure 13 compares several slope failure evolutions in terms of deviatoric strain for two soil friction angles (), 25° and 37°. As expected, the larger the internal friction angle, the weaker the movement of the sliding mass was. For both friction angles, when the failure stage initiated, the shear strain was concentrated mainly along the soil–bedrock interface (Figures 13(a) and 13(d)). However, the sliding of the landslide body toward the slope toe caused compressive forces in the soil mass in the horizontal portion. Consequently, a significant difference in the shear strain evolution was observed in the horizontal section at the final stage between the two friction angles. When , because of the smaller effect of the soil mass moving from the slope part and the higher soil resistance in the horizontal part, a shorter shear band at the soil–bedrock interface in the horizontal part of the slope was observed. By contrast, at a small internal friction angle, weak resistance in lateral support caused the sliding mass to accumulate in the horizontal part at . After 7.5 s, the shear strain continued to expand along the soil–bedrock interface in the horizontal section. The maximum calculated runout distances were 54 m and 34.5 m for and 37°, respectively.
(a)
(b)
(c)
(d)
(e)
(f)
The displacement and velocity evolutions of monitored points for various friction angles are presented in Figures 10(d) and 10(i). For different values, the evolutions of displacement and velocity were similar in the 30° slope section. The maximum velocity decreased as the friction angle increased, and the velocity stabilized at approximately 12 and 9 s for the friction angles of 25° and 37°, respectively. Notably, in contrast to soil cohesion, soil friction angle had ostensibly little influence on postfailure characteristics (Figure 10(i)). Additionally, the sliding failure surface and time required to reach peak velocity (i.e., 5 s) were similar for different friction angle values.
(3) Effect of Soil Modulus. Regarding the influence of soil modulus on postfailure behavior, the evolutions of slope failure at three times points in terms of deviatoric strain contours for different Young modulus values () are illustrated in Figure 14. Similar to the influence of soil mechanical parameters, the shape of the final deposition profile of the landslide body was greatly influenced by changes in soil stiffness. The runout distance and deposition height substantially increased as the soil modulus decreased. When , the maximum simulated depositional height was 7.5 m at ; by contrast, this value was only 5.8 m when . The difference was due to the change in the Young modulus that affected the effective stress distribution and generation of excess pore water pressure during slope failure. Figures 14(c) and 14(f) depict the final configuration of two simulations with different values. Similar to the previous cases, the slope with a low Young modulus (i.e., ) exhibited multiple shear zones in the horizontal part during failure. The final runout distance was approximately 64 m, and Hh/Ld ratio was 0.34; by contrast, the maximum horizontal runout distance for the slope with was approximately 19.5 m, and the Hh/Ld ratio was 0.48. These observations imply that the final configuration is highly sensitive to the value of the Young modulus.
(a)
(b)
(c)
(d)
(e)
(f)
The influence of on the velocity and displacement evolution of two monitored points over time is presented in Figures 10(e) and 10(j). For both values of , velocity peaked at after the start of the slope failure simulation (Figure 10(j)). As expected, the numerical simulations indicated that both the velocity and displacement of sliding failure considerably decreased as the soil modulus increased. For example, the computed maximum velocity of material point T-1 was 12.2 m/s when and was larger than that of the baseline case (10.2 m/s), whereas this value was only 7.3 m/s when .
4. Sensitivity Assessment
Identifying the most influential factor affecting the kinematic behavior of slopes after failure is crucial in landslide risk assessment. A sensitivity analysis of variables, such as soil hydrological and soil mechanical parameters, was performed to identify the influence of each parameter on the maximum velocity and runout distance. Figure 15 presents the sensitivity analysis results in terms of the percent change in maximum velocity and runout distance versus the percent change in input parameters. The results indicate that all input parameters influenced the postfailure behavior of the slope. Among the aforementioned parameters, and phreatic level location (i.e., GWL) had the greatest effect on runout distance, followed by and (Figure 15(a)). For instance, an increase in soil cohesion by up to 600% from its initial value resulted in a decrease in runout distance by up to 100%. Similarly, the runout distance decreased by 80% and 51% when and were increased by 10 times (i.e., 900%) and 2 times relative to the baseline case, respectively.
Figures 15(b) and (c) indicate that all input parameters, including , , , , and GWL, considerably affected the velocity of the tracked material points. The sensitivity assessment revealed that the velocity of material point T-1 on the 30° slope section was mainly affected by , , , and GWL (Figure 15(b)). The velocity of T-1 decreased by 29% and 77% when and were increased by 100% and 600%, respectively. Compared with material point T-1, the sensitivity of input parameters to the velocity of T-2 differed substantially with larger percent changes in velocity (Figure 15(c)). Notably, the sensitivity of the friction angle on the kinematic behavior of the slope was the lowest. A reduction or increase in friction angle by 20% changed the velocity of T-2 by ±33%, and runout distance increased by 35% and decreased by 14% relative to the baseline case, respectively. This was attributable to the negligible contribution of tan () on the shear strength of the soil when the change in the internal friction angle is small. In summary, the results of the sensitivity assessment suggested that the kinematic behavior and soil movement of slopes at the postfailure stage are highly sensitive to phreatic level location (GWL), hydraulic parameters (), soil cohesion (), and the Young modulus ().
The landslides simulated in this paper are hypothetical cases based on the landslide case histories compiled from the literature; therefore, the present results discussed in this study still cannot be directly applied or related to real landslides. However, this study identifies several important factors, such as phreatic level location, , , and , which can have a significant influence on the postfailure characterization and kinematic behavior of shallow landslides. These findings highlight the importance of obtaining the reliable values of these parameters and quantifying their uncertainties for a better prediction in the kinematic behavior of shallow landslides.
5. Conclusions
In this study, MPM formulations were used in an investigation of the failure mechanism, kinematic behavior, and postfailure processes of unstable shallow slopes. The influences of hydrological conditions, soil hydraulic parameters, and soil mechanical parameters on the kinematic behavior of a failure mass were evaluated in terms of deviatoric strain, displacement, velocity, and runout distance. The effects of each parameter on velocity and runout distance were quantitatively compared in a systematic sensitivity assessment. The following conclusions were drawn based on the numerical results. (1)The numerical model validation indicated that the MPM can effectively model situations that feature extremely large deformation. The runout distance and the shape of the final configuration in the simulation were consistent with the experimental data, demonstrating that the Mohr–Coulomb model used in the present study can describe material behavior(2)The numerical results indicated that the slip failure surface for shallow slopes constitutes a translational slide in the slope section. Additionally, the failure mechanism in the horizontal section is complex and generates multiple plastic shear bands(3)The parametric study and sensitivity assessment results revealed that all input parameters (i.e., hydrological conditions and soil hydraulic and soil mechanical parameters) considerably influence the kinematic postfailure processes. Among these parameters, had the largest effect on runout distance. The magnitude of runout distance decreases with increasing soil cohesion. The behavior of a slope with noncohesive soil (i.e., ) resembles a flow process with relatively high velocity in the sliding mass(4)Regarding the kinematic behavior of slopes in terms of velocity and displacement, the numerical results indicated that when the failure mechanism is shallow, the magnitudes of maximum surface displacement on the slope and runout distance are similar. For all slope cases, the velocity of the material points in the middle of the slope reached peak values at the same time (i.e., ), indicating that the time of peak velocity is not affected by changes in input parameters. Furthermore, the movement speed of material points in the slope section was approximately two to three times that of the material points in the horizontal section
In reality, soil shear strength parameters may decrease under large strains. This may result in differences in the failure mechanisms and postfailure characteristics of shallow slopes. In future studies, the strain-softening behavior (i.e., strain-dependent and shear strengths) of soil should be analyzed to quantify the influence of these parameters on the kinematic behavior of slopes after failure.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no competing interests.
Acknowledgments
This work was supported by the Ministry of Science and Technology, Taiwan under Grant No. MOST 107-2628-E-002-003-MY3. The authors would also like to thank the anonymous reviewers for their comments and suggestions that improved the quality of the manuscript.