Modeling The Post-Failure And Run-Off For a Translational Landslide In Taiwan By Material Point Method

The material point method (MPM) is an extended nite element method that can be used to simulate large deformation scenarios. A massive translational landslide in Taiwan was adopted to validate the numerical technique as thorough investigations, including the digital terrain models (DTMs), laboratory experiments, and numerical analyses, were available in a forensic report. The MPM code Anura3D was used to mimic the landslide’s kinematics, post-failure, and run-off process. An unstable sandstone/shale interlayer was found to lead the slope sliding; therefore, the before-and-after DTMs from the report mentioned above were used to examine the run-off distance and deposition to determine the best t of reduced material properties for this layer. The sliding paths, displacements, velocities of the sliding can be evaluated by dividing the material points into several groups to differentiate the kinematic among them. Meanwhile, the simulations were compared with different numerical methods. The landslide duration and possible maximum safety distance were also assessed. This study has demonstrated that the MPM can analyze the large deformation, post-failure, and run-off distance of landslides. The critical timing of a slope failure is possible to be an essential index on national spatial planning for future disaster reduction.


Introduction
Global warming and extreme climate have profoundly in uenced our daily life in recent decades. They lead to many natural disasters, where landslide is one of the most severe ones. If a possible failure process of a landslide can be predicted in advance, an incident of such a disaster can be prevented accordingly. Understanding the post-failure behavior of a landslide becomes critical in slope stability assessment for disaster reduction. A successful prediction of a post-failure behavior provides useful information, including both kinematics and run-out distance. Warning given to the potentially vulnerable area may be issued before the incident to protect people's lives and properties. (e.g., Guo  Limit equilibrium is a traditional method used for slope stability analyses. It is a fast and straightforward method with an assumption that materials are perfectly plastic. Force and/or moment equilibrium can be satis ed depending upon the method selected for each sliding surface. It re ects the slope's limiting condition where the slope is about to fail. It also provides the so-called factor of safety associated with this slope failure. Researchers have commonly used the limit equilibrium to carry out slope stability analyses or back analysis to estimate soil parameters at the triggering of a landslide (e.g., Fredlund and Krahn, 1977;Duncan, 1996;Krahn, 2003).
Finite element method has been a widely used for decades, where the balance of forces and deformation continuity can be satis ed. Besides, constitutive models can be applied to mimic different material behaviors. Due to the anisotropic and inhomogeneous natures in geotechnical engineering, the nite element method is capable of simulating the process of a boundary value problem for a speci c stress-strain relation. Duncan (1996) compared the limit equilibrium and the nite element methods for analyzing dams, embankments, and slopes. He showed the advantages and limitations of using the two methods in practical engineering problems. Troncone (2005) performed a landslide case using the limit equilibrium method (LEA) and nite element method (FEM). The analyzed result using the LEA showed that the average resistance along the sliding surface was between the peak resistance and residual resistance from laboratory tests. Using the LEA was found challenging to analyze the deformation and soil softening behavior of the translational type of landslides. In order to simulate a progressive failure process, an elasto-viscoplastic constitutive model implemented in the FEM was used.
Landslides are natural disasters with multiple processes in which the uid and solid mechanisms should be taken into consideration simultaneously; therefore, Eulerian and Lagrangian approaches were introduced (Durst et al., 1984;Mostafa and Mongia, 1987;Tinti et al., 1999). The Eulerian approach is a xed observation-only concerned with the uid properties at a speci c space and time point. It is capable of handling extreme deformation and avoiding mesh distortion. In the Lagrangian approach, the mesh is xed to the material geometry and moves with the deformed material so large deformation may lead to mesh distortion. However, large deformation problems simulated by the nite element method normally require longer computational time and lead to convergence problems. Uni ed methods were then proposed to avoid the numerical di culties on boundary setting and convective terms, and simulate large deformation problems without the constraint of mesh distortion. Cuomo  Discontinuous approach, so-called the discrete element method and some meshless methods (e.g., smoothed particle hydrodynamics, the particle nite element method) were proposed for large deformation problems. These methods are applicable to track particle movements and analyze It is suitable for any large deformation problems, especially landslides, because it can simulate kinematics movement during sliding. In addition, excess pore pressure and effective stress can be determined and monitored during the entire failure process. Mirada Larroca (2015) used the MPM to evaluate the slip surface and the post-failure behavior of soil excavation. This method was compared with different numerical methods, including discrete element method and limit equilibrium method, through simulation of granular ows and moving landslides (Ceccato et Fig. 1. The computation process is demonstrated by a set of Lagrangian material points passing through the Eulerian mesh and delivering the information of their properties, mass, velocity, and acceleration to the nodes around them to have new updates during deformation. Furthermore, with assigned constitutive models in material points, stress-strain relations, stress paths, and excess pore pressure changes can be obtained during computation. The aim of this study is to reproduce the possible failure scenario of a translational landslide in Taiwan using the MPM code Anura3D through back analysis. This case has been analyzed and simulated using The geological setting along the sliding direction was classi ed into the following strata, from the slope surface to the deepest layer, sandstone (SS), sandstone/shale interlayer (SS/SH), and the shale layer (SH), as shown in Fig. 2. It was believed that the sliding zone was somewhere between the sandstone/shale interlayer (SS/SH) and the shale layer (SH). A dip slope with an average dip angle between 14° to 16°.
The forensic report based on the in-situ investigation indicated that the sandstone in the studied area were weathered with joints and vertical cracks. These joints and cracks provided a natural shortcut for the rainfall to pass through and reach the relatively impermeable alternation of SS and SH layer. There was no rain before the landslide; however, the in ltrated water had secretly pushed the slope to edge. The in ltration raised the groundwater table and gradually reduced the strength in SS/SH layer. This softened layer lay on the extremely low permeable SH layer; thus, the interface between SS/SH and SH layers was judged to be the landslide sliding zone.
Moreover, the investigation also reported that most of the ground anchors used to stabilize the slope was rusted because of the groundwater seepage. As part of the anchors was overloaded and about to fail, eventually, almost 90% of the anchors were found ruptured. This was believed to be the trigger of the landslide. Because of the strength reduction in the sliding zone SS/SH layer, the slope nally became unstable. The rusting on the anchorage system could not stabilize the slope and led to 75 m displacements, and over 210,000 m³ sliding mass collapsed on the freeway just in a short period of time. Several vehicles failed to escape in time.
In this paper, section AA' (shown in Fig. 3) displays the topography before and after the landslide was selected for the MPM simulation. The digital terrain model (DTM) was obtained through the aerial photograph and orthoimage. The DTMs were obtained at the time before the freeway was constructed, after slope was excavated, after freeway viaducts were constructed, and after the landslide took place in 2010. The complete DTM data help this study adequately to understand the initial state of the landslide and its post-failure topography as shown in  Table 1.  The stratum of the studied slope was determined by extending the lower half geological pro le in section AA', so the thickness of SS/SH layer was set as 5 m. Strength reduction was applied in the simulation to mimic the water accumulation above the impermeable SH layer in SS/SH layer, which means the dry condition (or one-phase) was set up in MPM. This method accounted for both effects of reducing the material strength properties and raising the groundwater level.
The simulation aims to evaluate and discuss when all ground anchors were destroyed in the studied section; therefore, no anchoring system was considered in the model. For the boundary conditions, rollers were assigned on the right and left sides of the model to constrain the horizontal displacements, while hinges were set along the bottom boundary to constrain both vertical and horizontal displacements. Rollers were also placed along the top boundary to constrain vertical displacement. Figure 5 shows the mesh with triangular elements, and one material point was considered within each triangular element. In this case, the model veri ed by convergence analyses has 25,844 elements, 13,185 nodes, and 3,646 material points in total, which were assigned in the computational domain. Figure 6 illustrates the initial vertical stresses of the studied area. DTM data were collected before and after the landslide. The post-failure topography constructed by DTM was used to compare with the simulations, as shown in Fig. 7. The calculation procedures and numerical parameters are described as follows. In the rst phase, gravity loading stage, a relatively high damping coe cient, 0.75, was applied linearly from 0 to 1 to generate the quasi-static initial stress. In the second phase, a relatively low damping coe cient, 0.01, was selected so that the landslide failure can be identi ed when a highly dynamic behavior was observed. The entire simulation began as an initially unstable situation due to reduced material properties in the SS/SH layer.

Simulation Procedures
Each time step was set as 0.5 sec. A total of 80 calculation steps were needed for a 40-sec event. Finally, a strain smoothing algorithm was applied during the simulation to avoid kinematic locking. The strain smoothing algorithm controlled the increase of ctitious stiffness and smoothening the volumetric strains in the elements nearby and led to more accurate displacement results.
According to the summary from Table 2, the strength reduction in the SS/SH layer in this study started from c = 0 kPa, φ = 12 o to c = 0 kPa, φ = 6 o . Seven simulation cases were performed. After comparing the simulation results with the DTM of the post-failure stage, the simulation with c = 0 kPa and φ = 10 o best ts the corresponding run-off distance and deposit of the sliding mass shown in Fig. 7. It is noted that all these parametric studies have the same settings described in this study, except the SS/SH material properties of the studied slope. Figure 8 illustrates the simulation results of the freeway No. 3 landslide by MPM, presenting the displacement at 3.5 s, 6.0 s, 9.5 s, and the nal step, respectively. Table 2 Final maximum displacement of each simulation.

Results And Discussions
Six groups of material points were selected to assess their sliding paths during the landslide. Each group consists of six material points. As shown in Fig. 9, Group 1 is located near the toe and ground surface of the slope, while Group 2 to Group 6 have materials point inside the slope. In each group, the average displacements and velocities at various time instances were calculated.

Displacements and Velocities
For the simulation case that best ts the eld observation reported in the forensic report (c = 0 kPa, φ = 10 o ), the displacement and velocity correspond to time histories are shown in Fig. 10. In the beginning of the landslide event, say within 10 sec, Group 1 had more accumulated displacement than other groups. This also re ects on its behavior of velocity time history, where Group 1 had larger velocity in the initial movement of the slope. The front of the landslide (Group 1) reached its maximum sliding velocity around 13 m/s, 6 sec after the landslide started. It can also be seen from Fig. 10, the displacements of Group 1 to Group 5 show a similar pattern. However, due to the limitation of the two-dimensional model in this study, the material points in Group 6 were stuck by the deposition of the sliding mass in front. In a threedimensional simulation, the material points in Group 6 may keep sliding along the bedding plane to the north after being stuck by material points in front. The soil mass velocity plays an important role in the landslide. The dip slope failure usually causes severe damage, where the translational slide takes place suddenly, with the sliding velocity that can easily reach 100 km/hr or 28 m/s (Hung, 2010). The horizontal displacement, vertical displacement, maximum velocity, and average velocity of each group are listed in Table 3. The average velocity of the landslide is around 5 m/s, and the maximum velocity can be more than 12 m/s. The maximum velocity can reach 12.34 m/s in Group 1 and have an average of 9.00 m/s for Group 2 to Group 4. Table 3 The displacement and velocity data as c = 0 kPa and φ = 10 o for the SS/SH layer.  Figure 11 presents the sliding path of the selected material points in each group. Group 6 had much less displacement than the other groups. As indicated earlier, it may be due to the limitation of plane strain modeling. Despite this limitation, the post-failure simulation by the material point method still gives a good agreement with the photographs taken after the landslide.

Comparison with Other Studies
Back analyses were carried out to clarify what strength reduction in SS/SH layer led to the best run-off distance and deposition quantity. Previous studies done by FLAC3D and PFC3D (Chang, Table 4 summarizes all the results simulated in two approaches, including the strength reduction in the SS/SH layer and friction coe cient in the sliding surface.

Duration Analysis
In the discussion of the disaster duration, Group1 to Group 3 took around 18 ~ 20 sec to reach a new equilibrium in the simulation, which is quite close to 17 sec simulated by Chang (2012) in PFC3D. Table 5 lists the landslide duration in three studies carried out by using PFC3D and MPM. Examining the location of the southbound and northbound lanes in Freeway No. 3 (Fig. 12), the aerial view before the landslide and the side view after the landslide shows that the freeway has roughly around 50 m wide at section AA'. The simulation shows that at time 3.5 sec the landslide reached the edge of the Freeway No. 3 southbound lane (shown in Fig. 8). The sliding mass took only 2.5 sec and 3.5 sec to cover the southbound and northbound lanes.
Since the speed limit of the Freeway No. 3 Cidu section is 90 km per hour (25 m/s), combining the critical timings to cover the southbound and northbound lanes with the vehicle speed limit, a possible safety distance can be estimated. As a result, the possible maximum safety distance for vehicles passing through the southbound lane and not getting involved in the landslide is only 62.5 m, while it would be 87.5 m for the northbound lane. Through the MPM analysis, the run-off distance with the time history of the landslide was reproduced, and the critical timings were identi ed. Hence, knowing the landslide process and duration helps disaster reduction and future national spatial planning

Conclusions
The MPM was used to mimic the process of the Freeway No. 3 Landslide, where run-off distance/velocity time histories, the in uence area, failure deposition, and sliding paths were determined. In this study, the rock slope was simulated as the material points. The sliding zone was de ned by the full immersion in all SS/SH layers instead of assigning a friction coe cient on the sliding surface between the SS/SH and SH layers. The following conclusions can be drawn from the ndings of this study.
1. The simulated result shows that as the strength reduction in all SS/SH layer reaches 10 o in friction angle with 0 kPa cohesion, it has the best correspondence in deposition and run-off distance compared with the DTM measured at the post-failure stage.
2. Through the simulation of the landslide process using MPM, the critical timings and the duration of the landslide can be analyzed. It can serve as a vital index on slope failure and national spatial planning for disaster reduction.
3. According to this study and other literatures indicated, when dealing with slope stability design along the freeway or any road, its normal section should not be the only design section. The orientation of the stratum, including strike and dip directions, should also be considered.   Geometry and geological strati cation in MPM model.  The landslide process simulated by MPM at the displacement at 3.5 s, 6.0 s, 9.5 s, and the nal step, respectively.

Figure 9
Selected MPs in each group to represent the sliding mass: (a) initial stage before applying gravity loading; (b) post-failure stage as c = 0 kPa and f = 10o for SS/SH layer.

Figure 10
Calculated displacement and velocity of group 1 to group 6 as c = 0 kPa and f = 10o for SS/SH layer: (a) displacement versus time; (b) velocity versus time.