Analysis of the failure modes and development of landslides using the material point method

Landslides and mass movements are natural catastrophes that frequently occur in mountainous regions and urban areas, affecting the population and generating significant economic losses. In most cases, movements are caused when destabilizing forces exceed the resistance of the materials, by the alteration of the characteristics of the slope due to changes in the geometric factors (height, inclination), the materials conditioning factors or intrinsic factors (geology, hydrogeology and geotechnics) or by triggering factors (dynamic loads, variation in hydrogeological conditions, climatic factors, variations in geometry and reduction in resistance parameters). According to Skempton & Hutchinson (1969), mass movements develop in three stages: before the failure, the failure, and the landslides after it. Most research focuses on the study of the first two stages, i.e., the prevention of these movements by calculating the safety factor or the probability of failure. Some authors have focused their studies on the development of methodologies that quantitatively and qualitatively evaluate the susceptibility to landslides (Park et al., 2013; Malet et al., 2009; Simoni et al., 2008; Kanungo et al., 2006). Others have focused on the study of deflagrating factors and triggering conditions, including their description, classification, and time of occurrence (Aristizábal et al., 2016; Borfecchia et al., 2016; Aghda & Bagheri, 2015). Giupponi et al. (2015), Li et al. (2010), Zêzere et al. (2008), and Uzielli et al. (2008) have estimated the loss degree of an element set exposed to the occurrence of a landslide, assessing the vulnerability or propensity to such loss. All these works come together in the common objective of calculating the probability and severity of a mass movement, that is, in the calculation of the risk from the product between the probability and the consequences of an event. In order to improve the risk analysis and prevent mass movements, physical and mathematical models were developed to simulate the behavior of this type of phenomenon. Among the mathematical models, the Limit Equilibrium Method (LEM) is better known. It concentrates on the calculation of the safety factor using material resistance theories considering soil mass as a non-deformable rigid body. On the other hand, the Finite Element Method (FEM) has also been used as a more accurate tool for solving most problems in solid and soil mechanics (Nazem et al., 2006; Farias & Naylor, 1998). These methods concentrate only on the slope behavior before failure and at the beginning of instability, without considering the study of movement development, i.e., events after failure. However, to model mass movements after failure with good precision, a numerical method is necessary to accompany the movement of the simulated body and to describe its behavior, i.e., a method capable of simulating large distortions and deformations. Currently, the FEM is Abstract Mass movements are frequent natural phenomena and especially dangerous due to increased population and irregular settlements in mountainous areas. Carrying out studies using methods that allow the numerical evaluation of their behavior is essential to mitigate and prevent these events’ possible impacts. The traditional Limit Equilibrium and Finite Element methods fail to reproduce the problem from the beginning of the movement until the end of its development, so it is necessary to use new formulations. In this work the Material Point Method was selected to evaluate these movements in the stages of failure mode formation and development of landslide through an analysis of the soil shear strength and deformability parameters of a slope using the elastoplastic model with Mohr-Coulomb failure criterion. Maintaining a geometry and assuming soil mass behaves like a fluid, a single strength parameter and a single deformability parameter are analyzed to understand their influence before and after the failure.


Introduction
Landslides and mass movements are natural catastrophes that frequently occur in mountainous regions and urban areas, affecting the population and generating significant economic losses. In most cases, movements are caused when destabilizing forces exceed the resistance of the materials, by the alteration of the characteristics of the slope due to changes in the geometric factors (height, inclination), the materials conditioning factors or intrinsic factors (geology, hydrogeology and geotechnics) or by triggering factors (dynamic loads, variation in hydrogeological conditions, climatic factors, variations in geometry and reduction in resistance parameters). According to Skempton & Hutchinson (1969), mass movements develop in three stages: before the failure, the failure, and the landslides after it. Most research focuses on the study of the first two stages, i.e., the prevention of these movements by calculating the safety factor or the probability of failure.
Some authors have focused their studies on the development of methodologies that quantitatively and qualitatively evaluate the susceptibility to landslides (Park et al., 2013;Malet et al., 2009;Simoni et al., 2008;Kanungo et al., 2006). Others have focused on the study of deflagrating factors and triggering conditions, including their description, classification, and time of occurrence (Aristizábal et al., 2016;Borfecchia et al., 2016;Aghda & Bagheri, 2015). Giupponi et al. (2015), Li et al. (2010), Zêzere et al. (2008), and Uzielli et al. (2008) have estimated the loss degree of an element set exposed to the occurrence of a landslide, assessing the vulnerability or propensity to such loss. All these works come together in the common objective of calculating the probability and severity of a mass movement, that is, in the calculation of the risk from the product between the probability and the consequences of an event.
In order to improve the risk analysis and prevent mass movements, physical and mathematical models were developed to simulate the behavior of this type of phenomenon. Among the mathematical models, the Limit Equilibrium Method (LEM) is better known. It concentrates on the calculation of the safety factor using material resistance theories considering soil mass as a non-deformable rigid body. On the other hand, the Finite Element Method (FEM) has also been used as a more accurate tool for solving most problems in solid and soil mechanics (Nazem et al., 2006;Farias & Naylor, 1998). These methods concentrate only on the slope behavior before failure and at the beginning of instability, without considering the study of movement development, i.e., events after failure. However, to model mass movements after failure with good precision, a numerical method is necessary to accompany the movement of the simulated body and to describe its behavior, i.e., a method capable of simulating large distortions and deformations. Currently, the FEM is the most used tool for the calculation of deformations in geotechnical problems; however, it is not possible to obtain reliable results on large deformation problems using this method with its traditional formulation due to excessive mesh distortions (Nazem & Sheng, 2005). To overcome this difficulty, other methods were proposed, such as the Discrete Element Method (DEM), Galerkin Element Free Method (EFMG), Smoothed Particle Hydrodynamics (SPH), Arbitrary Lagrangian-Eulerian Method (ALE), and the Material Point Method (MPM).
The MPM is a tool capable of solving problems of large deformations related to geotechnical engineering, where the material point is seen as a representative volume element. The MPM uses constitutive models based on the continuum mechanics, such as the elastoplastic models Mohr-Coulomb, Modified Cam-Clay, among others. In addition, the use of a background mesh allows the implementation of boundary conditions similar FEM, and compared to other mesh-free methods, MPM has computationally efficient features (Abe et al., 2013).
The MPM method has been used to model different geotechnical problems such as foundations (Lorenzo et al., 2013), anchors (Coetzee et al., 2005), problems of granular flows (Więckowski, 2003), fault deformations (Johansson & Konagai, 2007), analysis of soil flows propagation induced by earthquakes (Konagai et al., 2004) and geomembrane response to settlements (Zhou et al., 1999). Other authors have worked in mass movements, such as González Acosta et al. (2018), who studied the effect of a landslide colliding with a rigid wall, considering multiple initial conditions to identify the most critical case. Vardon et al. (2017) conducted several slope failure simulations that included the initiation and development of the movement to better quantify the risk and consequences of these events. Llano Serna et al. (2016) demonstrated the capabilities of the MPM to evaluate landslides and their behavior after failure, numerically validating the failure of the Tokai-Hokuriku highway in Japan and the Vajont landslide in Italy. Bhandari et al. (2016) adopted the MPM to simulate the progressive failure of a slope caused by an earthquake. Gabrieli & Ceccato (2016) simulated the impact of a dry granular flow on a rigid wall and compared their results with DEM, obtaining good approximations with both methods. Llano Serna et al. (2015) simulated the landslide of an urban slope to calculate velocity and energy variables to quantify the vulnerability of structures and people subjected to impact, and Mast et al. (2014) used the MPM to simulate gravity-driven slope failure to assess the interaction of these events with the built environment.
As in the mentioned studies, this work uses the MPM to evaluate mass movements, specifically in the stages of formation of the failure mode and development of the movement through an analysis of the soil shear strength and deformability parameters of a slope using the elastoplastic model with Mohr-Coulomb failure criterion. The results obtained in the failure are compared with the results calculated by LEM. Finally, the estimates of velocity and displacement are evaluated and compared with different values of Young's modulus and cohesion to understand, at a general level, their influence on the formation and development of a mass movement.

Materials and methods
As previously stated, the MPM developed by Sulsky et al. (1994Sulsky et al. ( , 1995 was used to perform the simulations. The MPM brings together ideas and procedures of the Particle in Cell Method (PIC) and FEM. Material bodies are discretized in a collection of particles not connected between them that transport a mass whose value is kept fixed to guarantee the mass conservation. Other parameters necessary to define the body's state, such as stress, density, and the history variables, are also assigned to the material points (Zabala & Alonso, 2011). The interaction between the particles is performed in the nodes of a stationary Eulerian computational mesh like those used in the FEM, which remains constant for the entire calculation eliminating the problem of distortion ( Figure 1) This mesh is used to determine the governing equations' incremental solution by means of an Eulerian description (Al-Kafaji, 2013).
In this method, the equations of motion are solved in the background mesh that covers the entire domain of the problem. In each analysis step, the quantities transported by the material points are interpolated to the mesh nodes, using the functions associated with them as in the FEM. The boundary conditions are imposed on the mesh nodes, and the equations are solved incrementally in it. Then, the magnitudes of the variables in the material points are updated using the weighting of the nodes' results, using the same functions of form again. In the MPM, the mesh information is not required in the next steps of analysis, so it can be discarded (Zabala, 2010).
To evaluate the influence of the strength and deformability parameters on the formation of the failure mechanisms and on the development of a mass movement using the MPM, a slope was designed with the geometry presented in Figure 2 and  Table 1 (previously used in studies by Alelvan et al., 2020). As this method uses constitutive models based on the continuum mechanics, the elastoplastic model with Mohr-Coulomb failure criterion in the slope and the linear-elastic model in the foundation were used.
All analyses were performed with the ANURA3D® software, a 3D implementation of MPM, used to simulate the phenomenology involved in soil-water-structure interaction and large deformation problems. This tool was developed by Anura 3D  Redaelli et al. (2017), and .
Initially, it is necessary to evaluate the influence of the number of elements of the bottom mesh, and the number of material points in each cell, i.e. a sensitivity study of the discretization in the MPM should be performed. In this case, a series of analyses were carried out in which the slope instability due to its own weight is simulated, this being discretized with five types of meshes of different cell sizes, each one evaluated with one, four, and eight material points. The geometry used was the same as presented in Figure 2. The affix of each of the simulations (expressed by the MXMPY code where "M" is the mesh, "X" the reference to the size of the elements, "MP" the material points, and "Y" the quantity of them), the discretization used and the analysis time is presented in Table 2.
When the number of elements increases by decreasing the cell size, the analysis time rises significantly. n the other hand, the visual result is also influenced because as the cell size increases, the distance traveled by the unstable mass on the foundation decreases regardless of the number of material points. When the cell size is smaller, the distance reached by the mass is much greater, behaving practically the same with the three amounts of material points analyzed. Figure 3 shows the slope's discretization and the final distance reached by the mass when the type of mesh changes.
Based on these results, the M5MP4 mesh was selected to perform the following simulations. Although this simulation's analysis time was very long, the movement's development was much more detailed. It should be noted that the mesh has a strong influence on the results, which is more important than the number of material points assigned.

Estimation of an initial failure surface
To estimate an initial failure surface, the slope of Figure 2 was analyzed with the LEM proposed by Morgenstern & Price (1965) using the GEOSLOPE software. The LEM is one of the most used to analyze the stability of a natural or artificial slope by calculating the safety factor. This method is based on statics principles, i.e., the static equilibrium of forces and moments without considering the soil mass's displacements (Fredlund & Rahardo, 1993). With this method and the properties of Table 3, the calculated safety factor was 0.66, meaning the slope in the initial conditions was unstable. The cohesion was then increased to 1.8 kPa, the minimum value to guarantee the slope's stability with a safety factor equal to 1 (Figure 4).
When running the MPM simulations, it was observed that the minimum cohesion value to keep the slope stable was 1.25 kPa. These results being different from those obtained with LEM due to the estimates of stresses and deformations in the material that this method does not consider (van Asch et al., 2007, quoting Bromhead, 1996.    Different situations were raised from these divergences to be analyzed with MPM through two stages: formation of the failure mode and development of the mass movement. The final properties of the materials for each case are presented in Table 4.
Each simulation was carried out in two stages: in the first one, the initial stress state was generated, and in the second, the failure was simulated due only to the effect of gravity.

Results
In all cases studied were selected three material points on the foot, center, and crown of the geometry as shown in Figure 2. For the ease of the reader, the results of the points located in the crown are presented, which exemplify the general behavior of each simulation.

Analysis of the failure mode
In this stage the first instants of the simulations are studied, where it is determined if the movement is progressive,  Analysis of the failure modes and development of landslides using the material point method that is, if a slide is going to occur or not. For this purpose, two types of analysis were performed: one varying the cohesion values and the other increasing Young's modulus values.
In Figure 5, it is possible to observe the small displacements of the material points for each simulation of the cases with different cohesion values. When cohesion decreases, the displacements' magnitude is increasingly greater and much more concentrated in the region close to the slope face, defining more precisely what will be the failure surface and the mass that will move after that failure. In this case, the boundary effect was considered as not influence the results. Figure 6 shows that with an increase in cohesion, displacements are smaller. As the slope approaches the failure, the displacements, and velocity increase, although they stabilize again after approximately 0.7 s.
The displacements and velocity recorded with the cohesion of 1.25 kPa to analyze the Young's modulus variation are shown in Figure 7. The magnitude of the displacements decreases when the module increases, i.e., when the material has a more rigid behavior, the particles' internal accommodation is reduced, bearing in mind that the slope is in the failure's imminence. Likewise, velocities are higher with low modulus.

Failure
To analyze the behavior of the slope at failure, accelerations before and after instability were calculated. Figure 8 illustrates the selected points' behavior on the sloped crown for different Young's modulus values and cohesion of 1.25 kPa (imminent failure). In this figure, it is possible to observe two successive increases of acceleration at the beginning of the analysis, followed by braking generated by the mass's accommodation that prevents other displacements. The movement then slows down and quickly returns to 0 m/s 2 , just when the slope is stable again. It is also possible to observe that with the lowest modulus of 30 MPa the material point reaches the highest acceleration and takes longer to stop its movement compared to the other points; otherwise, it happens with the material point that represents the highest module of 300 MPa, which reaches a much lower acceleration and therefore stabilizes before the others. Figure 9 shows the accelerations of the material points analyzed with a cohesion of 1.0 kPa. The initial behavior of these points resembles the cohesive action of 1.25 kPa. In this case, the processing time was longer. Only after 3s of analysis, it is possible to observe changes in the acceleration, which remains constant after reaching its first peak at the     beginning of the movement. The most critical decelerations occur in the time interval 3.0s-3.65s, being -1.31 m/s 2 the most significant magnitude reached by the material points of module 30 MPa and 60 MPa. After 4.25s, all points remain stable at an acceleration of 0 m/s 2 .
The comparison of behavior with the two cohesions and the four modules simultaneously (0 s-1.2 s) is presented in Figure 10. In all these graphs it is possible to observe the formation of the failure mechanism for slopes with cohesions equal to 1.0 kPa and the delay of this formation for slopes with the cohesion of 1.25 kPa, just when the acceleration is reduced. The increase or decrease in cohesion determines the number of points that are being plastified, this number being lower when cohesion is greater.
The difference in each case lies in the redistribution of stresses through strains that generate an acceleration reduction. For slopes with 1.25 kPa of cohesion, this distribution (where the points do not being plastified) is enough to produce a braking process that stops the sliding, unlike slopes with 1.0 kPa, where the points yielding and, therefore, the reduction of acceleration is not enough (still with positive magnitudes) to prevent movement.  In Figure 10, it is also possible to observe how the modulus is responsible for the increase or decrease of the acceleration in the formation of the failure mode: the greater the modulus, the lower the acceleration, both in the maximum peak reached and in its decrease until reaching the failure.

Development of the landslide
The Young's modulus does not influence the distance covered; otherwise, cohesion strongly influences it and is present in Figure 11. This mass reaches a greater distance when cohesion decreases and is barely displaced when cohesion is equal to 1.0 kPa, with sensitivity to cohesion being one of the essential characteristics of the cases studied. Figure 12 compares the points on the sloped crown for cases with variation in cohesion after failure. In these analyses, there is no significant difference in the displacement and velocity results when Young's modulus increases, as the figure show for an example with E=30 MPa and E=300 MPa. In all cases, small increases in cohesion generate large decreases in the displacement magnitude. Likewise, this behavior is observed in the velocity of the three points, the greater the cohesion, the lower the velocities' magnitudes. Figure 13 shows that Young's modulus does not affect the magnitude of displacement either velocity. On the other hand, cohesion strongly influences the magnitude of displacement, velocity, and mass stabilization.

Conclusions
The material point method used in this investigation was efficient to analyze the three stages of a mass movement: before the failure, at failure, and the development of landslide. With this method, it was possible to evaluate the influence of strength and deformability parameters.
Despite its simplicity, the elastoplastic model with the Mohr-Coulomb failure criterion, allowed to identify the formation of the failure mode and the development of landslides by analyzing the influence of strength and deformability parameters. With this constitutive model, the behavior of the slope was evaluated using MPM and LEM methods. The results obtained with MPM are not congruent with the results obtained with LEM due to estimates of stress and strains in the material that this method does not consider. According to the MPM results, the slope would remain stable for a more extended time. Likewise, the unstable mass that would slide is overestimated by the LEM, which calculates a failure surface at the slope's foot and not in the body as approximated by the MPM and shown in Figure 14.  The analysis of the strength and deformability parameters showed that Young's modulus has no influence on the development of the movement and that only the strength parameters (cohesion) intervene in the landslide's behavior. After the failure, the soil mass flows and the shear strength properties are the only ones that influence the movements, being like the viscosity in a fluid.
The formation of the failure mode and the development of a mass movement is susceptible to variations in cohesion. Small increases in this property's values define the stability or instability of a slope (Figure 15). It is possible to identify with the acceleration analysis if a property has, or not, importance in the formation of the failure mode and in the distance reached by the movement once this is developed. Therefore, the acceleration analysis is one of the most important in these results because it is the one that allows identifying entirely if the slope is stable or not. For example, the acceleration and deceleration mechanism is susceptible to small variations in cohesion since this property determines whether the points are being yielded (Figure 15).
Deformations on the slope before failure decrease as Young's modulus increases, as is expected in materials with high stiffness. On the other hand, once the slope slides, i.e., for Analysis of the failure modes and development of landslides using the material point method slopes with safety factors less than 1.0 the distance traveled and velocity are slightly influenced by Young's modulus, due to in a perfectly plastic elastic model, once the material is broken, the rigidity of the system tends to zero regardless of the modulus adopted. For this reason, after the break, the displacement and velocity results do not differ significantly from the module variations of up to 10 times the initial value.
The results obtained in this study are also being evaluated in other geometries with other types of materials (including those with different friction angles and higher cohesion values). They, therefore, will be part of future publications.