Study of landslides and soil-structure interaction problems using the implicit material point method

Mesh based methods such as the finite element method (FEM) are the most usually used techniques for analysing soil-structure interaction problems in geotechnical engineering. Nevertheless, standard FEM is unable to simulate large deformations and contact, hindering the realistic simulation of rotational, sliding, pull-out and overturning behaviours. Contemporary ‘particle ’ methods, such as the material point method (MPM), do not use a mesh to discretise the material, allowing large deformations to be simulated. In this paper, a recently developed technique to simulate contact using implicit MPM is tested by simulating soil-structure interaction problems and a landslide. First, the behaviour of a retaining structure is studied during the impact of a mass of soil for different foundation conditions. Then, a landslide triggered by construction procedures is analysed. This new formulation allows the development of deep and shallow complex failure mechanisms (a combination of passive and active soil failures) and therefore the means to assess the consequences of a slope failure.


Introduction
Landslides are natural hazards in which a large mass of debris, mud or rock moves down-slope at a range of velocities.The triggering causes are diverse; usually a landslide is activated because of (i) the rapid loss of material strength, as occurs when the pore pressure in the soil increases during extreme rainfall (Wang and Sassa, 2003;Collins and Znidarcic, 2004;Iverson et al., 2015), or (ii) the rapid application of external loads, as occurs during earthquakes (Rodríguez et al., 1999;Nakamura et al., 2014;Li et al., 2012).Nonetheless, the occurrence of slow landslides, in which movement can take place over several days or years, is also possible.These landslides are usually undergoing creep (Van Asch and Van Genuchten, 1990;Furuya et al., 1999), or the complex geometric and geologic characteristics of the site prevents fast sliding of the soil (Rico et al., 1976).In Fig. 1, typical slope failures are illustrated.Fig. 1a shows a retrogressive failure, in which an initial shear band (SB 1) develops, causing the material to slide (FV 1).Then, if the down-slope material cannot support the imbalance in loads of the backslope, a new SB develops, causing the slide of another block of material (FV 2).This process can repeat several times, causing multiple slope failures.In Fig. 1b, a translational failure is observed, in which a segment of material detaches and moves a short or large distance.In contrast to retrogressive failure, the mass of soil affected is larger and further sequential failures may not occur.Some of the most typical slope failure types, such as those depicted in Fig. 1, are described in Locat et al. (2011).
Regardless of the type of failure or triggering cause, damage to infrastructure and/or loss of human lives are a constant threat.As indicated by Kjekstad and Highland (2009), "there is an increased susceptibility of surface soil to instability as a result of more extensive human interaction of different kinds", which includes aspects such as uncontrolled land-use, increased forest clearance and climate change.Furthermore, based on data from the Centre for Research on the Epidemiology of Disasters (CRED), earth flows (i.e.landslides, avalanches, mudslides) caused up to 18,000 fatalities since the year 2000, with most of them occurring in America (North, Central and South) and Asia (CRED, 2009).In order to mitigate possible damage caused by landslides, it is important to develop tools capable of simulating real geotechnical scenarios in which slope failure or landslides can occur, including interaction with neighbouring structures.
The finite element method (FEM) is a numerical technique frequently used to study geotechnical problems, such as slope stability and the interaction with protection/retaining structures (e.g.Chen and Martin, 2002;Mishra et al., 2017).Nevertheless, in FEM, the connectivity between the mesh and the domain is essential and does not allow for the simulation of large deformations, thereby reducing the range of problems that can be studied.In contrast, with the development of particle-based methods (e.g.Harlow, 1957;Brackbill and Ruppel, 1986), it is possible to avoid the connectivity disadvantage of FEM and simulate large deformation problems.Some of the families of (continuum and non-continuum) methods used to simulate geotechnical large deformation problems are: the particle finite element method (PFEM) (Idelsohn et al., 2004), smooth particle hydrodynamics (SPH) (Pastor et al., 2014) and the discrete element method (DEM) (Tang et al., 2009).Nevertheless, the material point method (MPM) (Sulsky et al., 1994(Sulsky et al., , 1995)), has been demonstrated to be a better suited technique due to its simplicity and robustness (Tran et al., 2020;González Acosta, 2020).Nowadays, it is possible to find in the literature numerous studies simulating landslides similar to those shown in Fig. 1 (Wang et al., 2016b;Vardon et al., 2017).
MPM is formulated in a similar mechanical framework to FEM, keeping many of its attributes such as mass conservation.Moreover, via the incorporation of a contact algorithm (Bardenhagen et al., 2000), interaction between separate bodies is possible, allowing realistic simulation of soil-structure interaction problems (González Acosta et al., 2018;Müller and Vargas, 2019).To date, numerous papers have reported the use of MPM to simulate interaction between landslides and structures (Mast et al., 2014;Conte et al., 2019).Nevertheless, these attempts have generally failed to depict a realistic landslide simulation and the interaction with surrounding structures, due to the use of (i) unrealistic initial conditions (which leads to implausible failure triggers), (ii) simplified solution procedures (which ignore a realistic contact simulation), and (iii) poor stress recovery techniques (which reduces the accuracy of the simulations).
Recently, a procedure to simulate contact using an implicit scheme was developed (González Acosta et al., 2021), based on the explicit approach of Bardenhagen et al. (2000).This method was shown to be able to reduce (in most cases) the computational time for the same accuracy as an explicit scheme.Additionally, a method to reduce stress oscillation problems in MPM (González Acosta et al., 2020) was shown to increase the accuracy of simulations.This paper utilises these developments to validate the use of the implicit contact method for geotechnical problems.First, the theoretical background of implicit MPM is presented.Then, in Section 4, the equations and procedures to simulate contact using the implicit scheme are developed.Finally, in Section 5, two geotechnical applications are introduced to study the performance of the implicit MPM method in simulating soil-structure interaction problems.The first application is the simulation of a vertical cut which fails and impacts a retaining structure, and the second application is the simulation of a landslide that is triggered by a construction, which fails due to the inadequate design of the retaining structures, and subsequently interacts with multiple structures.

MPM description
MPM is a numerical technique which uses a similar mechanical formulation to FEM.The main difference between these techniques is that, in MPM, the material in the problem is discretized using (material) points that are able to move through the background mesh.The background mesh is still used to solve the mechanical equations, but instead of modelling the shape of the body, it covers the entire computational domain.The material points perform the same role as FEM Gauss integration points, i.e. stress/strain recovery and numerical integration to form nodal equations.They also hold the state variables of the bodies throughout the simulation, e.g.kinematic variables, stresses, and material properties.Since the points carry the information of the bodies, it is possible to reset the mesh to its original position after each solution step while the points stay at their latest positions, thereby avoiding excessive mesh distortion and allowing the simulation of large deformations.
Fig. 2 shows that MPM consists of three steps, of which two are similar to FEM.Step one is the mapping (indicated with arrows) and integration of mechanical equations (using material points instead of Gauss points) and step two is the solution/upgrade phase (Fig. 2a and b, respectively).At the end of step two, the final deformed mesh is obtained, and the material point variables are updated based on the final positions of the mesh nodes.During the third step (i.e. the convection phase), the new global and local positions of the material points are computed, the background mesh is returned to its original position, and the elements containing material points and boundary conditions are activated (as shown in Fig. 2c).

Mathematical background
The equations of the implicit MPM are summarised for a single body and assuming a small strain formulation, although large displacements are allowed.For a detailed elaboration of these equations, the reader is directed to Bathe (1996), Belytschko et al. (2013), andGonzález Acosta et al. (2021).
To derive the equation for static equilibrium in MPM, the principle of  virtual work is used.This principle states that the difference between the internal and external work should be equal to zero (⊓ = W int − W ext = 0), and is written as where ε are the strains, D is the material elastic matrix, V is the volume of the body, u is the continuous displacement field, ρ is the mass density, b are the external forces, and τ s are the prescribed loads at the boundary Γ.To consider a dynamic scheme, the inertia term must be added to Eq.
(1).Using D'Alembert's principle, the inertia forces can be included as part of the body forces, and Eq. ( 1) is then rewritten as where a is the acceleration of the body.Following standard FEM discretisation, Eq. ( 2) can be expresed in elemental matrix form as where the element matrices are where nmp is the number of material points p in the domain, ρ p is the material point density, M is the element lumped mass matrix, g is the acceleration due to gravity, x p is the material point position, W p is the material point integration weight, J is the Jacobian matrix, N is the matrix of element shape functions (SFs), and B is the strain-displacement element matrix.The terms a and u are the vectors of nodal accelerations and incremental displacements, respectively.The stress-strain relationships are fulfilled since the material point stresses σ p are computed using the constitutive relationship from the strains ε (which have been calculated from the incremental displacements u).
To solve the global equations as a function of time, the elemental equations (Eq.( 3)) are assembled into a global equation and Newmark's time integration method is used.Using this scheme, the nodal velocities and displacements at time t + Δt are where v is the vector of nodal velocities, and α and γ are the timestepping parameters, which are here taken to be α = 0.25 and γ = 0.5, to represent a constant-average-acceleration approach.After isolating a t+Δt from Eq. ( 9) and rearranging, this leads to Substituting Eq. ( 10) into Eq.( 8) results in By substituting Eq. ( 10) into Eq.( 3) and incorporating the Newton-Raphson iteration procedure, the equilibrium equation is written as where and where m is the nodal mass, K is the modified stiffness matrix, F kin are the kinetic nodal forces, and the left superscript k is the iteration counter in the Newton-Raphson iteration procedure.It should be noted that Eqs. ( 15) and ( 16) use the material point mass (m p ), which is already incorporated in the main formulation in terms of the material point density, weight and element Jacobian (i.e.m p = ρ p W p |J|).Finally, the material point stresses are defined as where σ p t represents the previous material point stress state and Δσ p = DBu.

MPM mapping, integration and stress recovery enhancement
Since standard MPM suffers from oscillation problems, the procedures elaborated in González Acosta et al. (2020) are used to mitigate integration inaccuracies and increase the accuracy of the stress recovery process.A double mapping technique is used, in which the (element) vectors of external and internal nodal forces are integrated using the generalized interpolation material point (GIMP) SFs (Bardenhagen and Kober, 2004), i.e.
where S is the element GIMP SF matrix, and ∇S is the straindisplacement matrix for GIMP SFs.The (element) stiffness matrix is computed, using double mapping (DM) procedures, as where ngauss is the number of Gauss integration points g in the element, smp is the number of material points with a support domain inside the element, S ip * is the local GIMP SF of node i (Charlton et al., 2017), and W FE is the weight associated with the Gauss point (as in FEM).Finally, the incremental material point stress is updated, using composite SFs (González Acosta et al., 2017), as and thereby where ∇N 2 is the strain-displacement matrix of CMPM SF gradients, and u ext and v ext are the vectors of nodal incremental displacements and velocities in the extended CMPM domain.It is important to remark that, if DM-GC is used, the rest of the state variables (i.e.mass, velocity, acceleration) are mapped using GIMP SFs.

Implicit contact
The procedures to simulate contact are based on the solution developed by Bardenhagen et al. (2000), where the velocity field of each body is computed and used to determine contact and frictional loads.In addition, the resulting velocities computed using Newmark's time integration scheme every Newton-Raphson iteration must be used to evaluate the new nodal velocity field, combined velocity and contact loads.The detailed formulation and steps followed to simulate implicit contact in MPM are elaborated in González Acosta et al. (2021).The first steps are the same as in Bardenhagen et al. (2000), in which the nodal velocities of each body and the combined velocities are used to detect contact as where where nb is the number of bodies in the computational domain, v i, bod and v p, bod are the nodal and material point velocities of each independent body, respectively, v i, C is the combined velocity accounting for the velocity of every body in the computational domain, and bod denotes the body.If contact is detected (i.e.Eq. ( 23) is satisfied), it is necessary to check if the bodies are approaching each other (compressing), i.e.
where n i, bod is the unit normal vector of the body, which is computed for two bodies as where the normal n i, bod considers the influence of both bodies, n is the individual normal direction to body bod at node i, m g is the accumulated mass at the Gauss position in the centre of an element, and x g are the coordinates of the Gauss point.Eqs. ( 27)-( 29) only consider contact caused by a maximum of two bodies; for multiple bodies at contact, readers are directed to Xiao-Fei et al. (2008).In Fig. 3, a sketch of the variables needed to detect contact is shown.Note that n m,A coincides with n k,B , which is not the case between n m,A and n k,B .This normal direction correction contributes to a better interaction between the bodies, mainly when contact at sharp corners occurs.If it is found that bodies are compressing (i.e.Eq. ( 27) is satisfied), the grid velocities in the normal direction are updated to ensure no interpenetration, i.e.
where *v i, bod t is the new nodal velocity corrected for normal contact.Finally, the normal contact force is computed as and the frictional force is computed as where F i, bod nc, t is the normal contact load, F i, bod stick, t is the test frictional load, fric, t is the real frictional load, and μ is the frictional factor that depends on the material properties.In Fig. 4, a sketch of the resulting contact forces is shown.Finally, after substituting Eqs. ( 31) and (33) into Eq.( 12), the final equation for equilibrium is written as

Iterative contact procedures
The previous equations have been developed at a particular time t.Nevertheless, during the iterative procedures and after using Newmark's time integration scheme, new nodal velocities are obtained, and the contact conditions must be re-evaluated to account for the new time t + Δt.After computing the incremental displacements (Δu) using Eq. ( 34), new nodal accelerations and velocities are obtained (Eqs.( 10) and ( 11)) which must be used to estimate the new contact velocities and contact loads.Using the nodal velocities from Eq. ( 11), new combined velocities must be computed as After computing the new combined velocities, a loop over the previous contact nodes needs to be performed to test if contact persists (i.e. if Eq. ( 27) is still true at time t + Δt).If contact does persist, the nodal velocities must be corrected, and contact loads must be upgraded.Hence and the frictional forces are

Application of the implicit contact method
The proposed implicit contact method is applied to two geotechnical problems.The first problem consists of a vertical cutting, which fails under self-weight and collides against a rigid wall.The second problem consists of a large slope, which fails, due to the removal of material and inadequate construction of retaining structures, dragging with it all neighbouring structures in the failure.In both problems, the progressive failure of the soil and the response of the structures after contact are investigated.Note that the mesh size and time step selected are based on González Acosta et al. (2021), in which an extensive investigation was performed to provide guidance on the most adequate conditions to allow accurate soil-structure interaction.Finally, note that, besides the use of a small strain formulation, large deformation can be obtained through the accumulation of infinitesimal deformations and an update of the geometry.

Vertical cut
A 2D elasto-plastic vertical cut has been simulated using a von Mises constitutive model incorporating post-peak softening (Wang et al., 2016c), which, after failure, collides against a rigid wall.As depicted in Fig. 5, the model is linear-elastic before yield and then softens proportionally to the plastic shear strain invariant ε pr (at a rate defined by the softening modulus H s ).Fig. 6 shows the generic problem domain (including dimensions in meters) and boundary conditions.This problem has been analysed twice, by considering two different geometries.In the first simulation, the wall is founded on the ground surface (i.e.s = 0),   and the foundation soil layer is shallow (h 2 = 0.25 m).In the second simulation, the foundation soil layer is deeper (h 2 = 1.0 m), and the rigid wall is founded at a depth of s = 0.5 m.The background mesh element size is Δx = Δy = 0.05 m and each element contains initially four equally distributed material points (i.e. at the centre of each element quadrant).
The unit weights of the soil and the wall are γ s = 18 kN/m 3 and γ w = 20 kN/m 3 , respectively.The elastic parameters for the soil are Young's modulus, E = 1.0 × 10 3 kPa, and Poisson's ratio, υ = 0.35, and the peak shear strength parameters for the soil cut and foundation soil are s pv = 10 kPa and s pf = 30 kPa, respectively.The residual shear strength and the softening modulus are the same for both soils, and are equal to s r = 3 kPa and H s = − 30 kPa, respectively.The elastic parameters for the rigid  wall are E = 1.0 × 10 4 kPa and υ = 0.38.
At the left and right boundaries of the domain, the nodes are on rollers to avoid displacement in the horizontal direction, whereas the nodes are fully fixed at the bottom boundary.The initial stresses in the domain were generated by linearly increasing the gravity load.During these steps, the kinematics were not considered and the material points stayed in their original positions (i.e.no movement was allowed).After the gravity load had reached its maximum value, the kinematics were included, and the material points were released.Then, due to the development of large shear stresses at the cutting toe and the low strength of the material, the failure was triggered.The failure process was modelled using a time step of Δt = 1.0 × 10 − 4 s for the simulation.Fig. 7 shows the interaction between the vertical cut and the rigid wall, by way of the plastic deviatoric strain and deviatoric stress distributions, for the simulation with no wall foundation.Note that the deviatoric stress range is fixed to a maximum of 25 kPa, allowing a Fig. 8. Slope initial state (a, b), collision instant t = 0.8 s (c, d), interaction at time t = 1 s (e, f), and final configuration (g, h), showing contours of plastic deviatoric strain and deviatoric stress.
better visualisation of the stresses in the vertical cut.In Fig. 7a and b, the initial condition of the simulation is shown (at the instant the material points are released).It is seen that the material points have not yet moved, although there is a stress concentration at the base of the cut.In Fig. 7c and d, the simulation at the instant of the contact (t = 0.7 s) is shown.At this stage of the analysis, the band already traverses the height of the vertical cut and a large block of soil is moving towards the retaining structure.In Fig. 7e and f, the simulation has reached time t = 1 s.In this case, the contact has already occurred and, due to the kinematic forces applied by the soil, the wall is being pushed and slides to the right.It is seen that a large portion of soil close to the contact zone undergoes plasticity due to the large contact forces developed.Finally, Fig. 7g and h show the final step of the simulation corresponding to time t = 3 s.By this stage of the analysis the rigid wall has fallen over after being pushed a considerable distance.Furthermore, at the base of the wall, the soil and the wall are not in contact and the soil has formed a circular scarp.This occurred during the fall of the wall, due to the rotational movement of the wall pushing the soil away and, due to the plasticity of the material, the soil not moving back into the gap.
Fig. 8 shows the results of the simulation with the wall foundation and a deeper soil foundation.In Fig. 8a and b, the initial plastic deviatoric strains and deviatoric stresses are shown, while Fig. 8c and d show the stresses and strains at the initiation of contact (t = 0.8 s).At this point, the results are similar to the previous simulation (as expected), but the contact occurs 0.1 s later.This time difference is attributed to the inclusion of the wall foundation, which has a small influence on the growing velocity of the failure mechanism.In Fig. 8e and f, the results during the contact at a time of t = 1 s are quite different from the previous simulation.It is seen that the rigid wall, rather than sliding, rotates because of the support given by the foundation soil.A small zone in the foundation soil at the back of the wall develops into a passive failure wedge, and the wall separates from the soil in areas where tension would occur, highlighting the advantage of using the contact algorithm.Finally, at the end of the simulation (Fig. 8g and h) the wall is not far from its original position and is able to prevent further movement of the failed soil mass.
In Figs. 9 and 10, the magnitude of the contact loads developed at the wall surface during the collision and the displacement of the wall, respectively, are shown for the second simulation.Fig. 9a shows the contact loads before the collision.In this case, the loads are caused only by the self weight of the wall and the soil, since the collision has not yet occurred.Fig. 9b shows the contact loads at the instant of the collision (t = 0.8 s); since the area of the contact is small and the velocity of the soil is large, the contact loads are high.Fig. 9c shows the evolving situation during the collision at time t = 1 s.In this case, since the soil begins to accumulate at the wall, the contact loads are distributed and the magnitudes reduce.It is also seen that, at the base of the wall on the left side, there are large forces due to the rotation of the wall.Finally, Fig. 9d shows the final position of the wall (t = 1.5 s).In this case, the contact  J.L. González Acosta et al. loads have reduced since the velocity of the material is equal to zero (i.e. the kinetic forces are equal to zero).It should be noted that the magnitudes of the vectors representing the contact loads have different ranges in each sub-figure to enhance visualisation.
In Fig. 10, the average horizontal displacement of four material points at the top of the wall (labeled mps) is plotted for both simulations.The selected material points are located at the centre of the wall.As observed, the displacement of the surface wall is at least twice as big as the displacement of the embeded wall.At the end of the simulations, the maximum displacements are d max S = 1.59 m and d max E = 0.67 m, for the surface and embedded walls, respectively.The simulations were run using an Intel Xeon E5-1620 processor and were completed in 9 h and 4 min for the problem considering the wall at the surface (using 11,040 material points), and in 12 h and 55 min for the problem considering the embedded wall (using 18,400 material points).
This first example has demonstrated that the contact formulation can be used to accurately assess the adequacy of a retaining wall structure, especially at ultimate failure conditions.This is a problem studied since the very conception of geotechnical engineering (Rankine, 1857) due to the need to stabilise slopes and earth fills.Nonetheless, methods which are commonly used have substantial limitations due to their deformation limits and their lack of an adequate contact formulation.

Landslide
The analysis of an initially stable slope, which undergoes several construction phases before failure, has been performed.The soil follows the same type of constitutive behaviour as in the previous example.Also, a depth dependent strength is included, in which the peak shear strength increases linearly with depth, from s p = 15 kPa at the soil surface to s p = 80 kPa at the base of the soil layer.The residual cohesion and the softening modulus are s r = 5 kPa and H s = − 30 kPa, respectively.The elastic parameters for the soil are Young's modulus, E = 1.5 × 10 3 kPa, and Poisson's ratio, υ = 0.40.The elastic parameters for the structural elements are E = 5.0 × 10 4 kPa, and υ = 0.35.
In Fig. 11, the main features of the slope are shown, including the construction stages and dimensions (in meters).Fig. 11a shows the initial slope, comprising a single homogeneous soil layer of constant depth H 1 .The element size is Δx = Δy = 0.5 m, and each element initially contains four equally distributed material points.The entire layer is supported by a rigid material, simulated by the fixed bottom boundary.The two vertical boundaries allow only vertical displacement.Fig. 11b indicates the excavation process, which involves three steps: (i) excavation 1 (exc -1) at a distance L 4 from the top of the slope, (ii) excavation 2 (exc -2) at a distance L 5 from the first excavation, and (iii) excavation 3 (exc -3) at a distance L 6 from the second excavation.Finally, Fig. 11c shows the wall foundation depths and the geometry of the building.The thickness of the columns and floors of the building are 0.75 m and 1.0 m, respectively.
The initial stresses in the slope are assigned in a single step using a static approach (i.e. by removing the kinematics), i.e. σ y = γ soil h y and σ x = σ y K 0 , in which γ soil is the soil unit weight, equal to 17.5 kN/m 3 , h y is the vertical distance between the soil surface and the material point, and K 0 is the coefficient of earth pressure at rest, equal to 0.7.Meanwhile, the material points are kept in their original positions (i.e.no movement allowed) in this step.After the first step, the material points are free to move, including during the construction steps.
The excavation procedure and the installation of the structural elements was simulated as follows.Each excavation was executed by removing columns of material points every 50 steps, in which each column had a width of one element.Fig. 12 shows an example of the excavation procedure.Note that this example is only illustrative, and that the dimensions are not the same as in the simulation.The excavations were executed sequentially, i.e. after excavation exc -1 was completed, exc -2 began.Each retaining wall is installed in a single step by changing the properties of ground material points to structural material points.The insertion of each wall is performed immediately after the removal of the last portion of soil in each excavation.The thickness of each retaining wall is equal to 1.5 m.To reduce the large loads caused by the inclusion of the wall in a single step, it is considered that the weight of the wall is initially defined by a reduced gravity of g ini = 0.3 g, which is increased in increments of Δg = 1.0 × 10 − 4 g each time step until reaching 1 g.The building is placed in a single step by adding material points in the domain.Similar to the retaining walls, the initial structure gravity is g ini = 0.3 g and this is increased in increments of Δg = 1.0 × 10 − 4 g each time step until reaching 1 g.The time step used for the simulation is Δt = 1.0 × 10 − 4 s.
Figs. 13 and 14 show the plastic deviatoric strains and the deviatoric stresses, respectively, during the landslide.The results are obtained from the instant after concluding the construction process (i.e.excavations, and installation of retaining walls and building) until the end of the landslide.Fig. 13a shows the instant immediately after adding the building, in which the plastic strains are not yet visible.Nonetheless, Fig. 14a shows a large concentration of shear stresses below the walls.Fig. 13b shows the simulation after 0.7 s.It is observed that the wall W-3 at the base of the slope begins to rotate and separate from the soil it retains.This is a consequence of both the weight of the building, which increases the horizontal pressure applied to the wall, and the lack of support on the downslope side.Moreover, the soil pushing the wall W-3 moves, thereby causing the movement of wall W-2, which, at the same time, causes plastic strains to develop in the soil beneath its base.Additionally, Fig. 14b shows that, at the soil-boundary interface at the base of the soil layer, large shear stresses develop.Note that, since it was the inclusion of the building that caused the wall W-3 to fall over, and thereby the subsequent soil displacements and the development of large shear stresses at the soil-boundary interface, it can be concluded that the building triggered the landslide.
Fig. 13c shows the simulation after 1.1 s.It is seen that the wall W-3 has rotated a considerable amount, and that the soil behind the wall begins to slide downwards.The movement of the soil triggers shear bands to develop at the base of each wall and at the soil-boundary interface.It is also seen that the building is tilting and that the walls W-1 and W-2 are separating from the retained soil.Fig. 13c shows that the plastic strains at the soil-boundary interface extend almost the whole length of the inclined boundary, and that the plastic strains at the base of the walls have grown significantly, reaching almost the soil surface.This is more clearly seen through the softened material zones in Fig. 14c.
Fig. 13d shows the simulation after 1.75 s.It is seen that the wall W-3 has fallen over, allowing the soil behind to slide more freely.The shear bands have extended to the soil surface, as well as between the middle and bottom walls where a shear band connects the bases of both walls, and numerous shallow failures (which exhibit typical rotational mechanisms) have also formed.Fig. 14d shows that the deviatoric stresses along the length of the sloping soil-boundary interface have dropped to residual values.However, it is seen that the deviatoric stresses are still large at the leading edge of the basal shear band, which indicates that the shear band is still advancing.
In Fig. 13e and 14e, it is seen that the shear band at the base of the soil layer now exceeds the length of the slope, to form a translational failure mechanism which encloses the multiple shallow failures.The walls W-1 and W-2, and the building, are being dragged by the soil, while the wall W-3 is being pushed away in front of the landslide.In addition, due to the large kinematic forces, another shear band begins to grow at the base of the lower horizontal section of the soil layer.Finally, Figs.13f and 14f show the final step, in which the soil has nearly reached static conditions.It is seen that most of the soil is in the plastic condition, and that it exhibits numerous failure surfaces.The shear band at the base of the horizontal section of the slope has grown enough to reach the soil surface at the end of the domain.Note that the right vertical boundary causes the vertical growth of the plastic shear band at this location, which would not have occurred had a larger horizontal domain been used.Finally, it is seen that all the structures have ended up buried in the ground.
Fig. 15 shows the displacements and velocities of the structures, which where obtained from an average of the results from four material points located at the centre of each structure (denoted with a solid circle in Fig. 15a).In Fig. 15a, the displacement of each structure is shown, with the labels A, B, C, D, E and F indicating the times corresponding to the results shown in parts (a) to (f) of Figs. 13 and 14.It is observed that, after 10 s, the structures have moved nearly 50 m.The vertical wall nearest the bottom of the slope (S3) has moved the furthest because it has fewest obstacles in front of it.Similar results are observed in Fig. 15b, in which the vertical wall S3 reaches the maximum velocity (close to 10 m/s), whereas the rest of the structures reach similar velocities (around 6 m/s).The type of information illustrated in Fig. 15 is necessary to evaluate the consequences of a landslide.This simulation was completed in 215 h with a total of 25,766 material points, and used the same processor as in the vertical cutting problem.
This second example has attempted to replicate the scenario of a landslide triggered by human activity rather than geological or external natural conditions.It has been demonstrated that, by adding the contact formulation to the MPM, an exhaustive study of such landslides can be performed.In this case, the initiation, development and final position of the landslide was observed and analysed.Based on the results, it is possible to evaluate the risk to nearby communities of severe damage, or to estimate the economic losses in such scenarios.

Discussion and limitations
The methods and analyses presented above provide the means to analyse local (e.g.vertical cut) and global (e.g.landslide) geotechnical problems.Nevertheless, it is pertinent to mention that the simulations herein could be further developed by considering the following: • Reliability assessment: a comprehensive analysis considering heterogeneous soil properties would provide information on the range of possible outcomes.This is particularly relevant in geotechnical engineering, in which the soil properties are found in nature to be randomly distributed.The reader is directed to Wang et al. (2016a) and Liu et al. (2019) for a further study of the implementation of spatially variable property fields (random fields) using the material point method.• 3D simulations: 2D simulations have been used here due to their adequate representation of a slope with a reasonably uniform crosssection, which translates to lower computational times and economical use of memory.Nevertheless, 2D simulations can only consider a cross-section of the full domain.This means that 3D features, due to geometry or heterogenous material properties, cannot be simulated.In the literature, geotechnical research using 3D MPM is available (Mast et al., 2014;Remmerswaal et al., 2019;Li et al., 2020), but it has remained scarce and at a relatively simple level due to computational limitations.• Validation: it has been demonstrated through several geotechnical benchmarks (e.g.bearing capacity, consolidation, water seepage) that MPM accuracy is satisfactory, encouraging the study of more complicated problems.Currently, MPM has been used to simulate landslides and slope stability problems.The results obtained seem to follow geotechnical engineering expectations but, due to the absence of methods to verify such problems (considering large deformations), or suitable field observations for validation purposes, the method  Acosta, 2020), and (iii) laboratory tests (Li et al., 2016).Nonetheless, these attempts cannot be considered full validation due to (i) the lack of scientific rigour of visual comparisons, especially during the failure process, (ii) the impossibility of performing large deformation simulations using FEM, and (iii) the difficulty of ensuring adequate initial conditions and (kinematic and stress) measurements.• Mesh dependency: it is well known that MPM simulations (as in FEM) are mesh-dependent when the material is undergoing strain- softening behaviour.Furthermore, contact simulation also experiences mesh-dependency effects since contact detection depends on the mesh-size.The consequences of having such dependencies are that (i) the initiation and growth of shear bands can be different for the same problem when different meshes are used, and (ii) the interaction between structures may be unrealistic when large mesh discretisations are used.To diminish these effects, it is recommended to consider (i) regularisation techniques to control plastic deformations (e.g.Galavi and Schweiger, 2010;Burghardt et al., 2012), and (ii) proximity contact techniques (González Acosta et al., 2021) to ensure that contact is activated once a certain proximity between the bodies has been reached.

Conclusion
An implicit MPM, that incorporates a recently developed stress oscillation reduction technique and an implicit contact formulation, has been utilized to simulate complex geotechnical problems.It was observed that the method was able to capture realistic interaction behaviour, including features such as the sliding and rotation of structures pushed by the soil, the development of passive and active soil failure mechanisms, and the combination of rotational and translational slope failures.Considering that the results exhibit the expected behaviour, MPM can be used for the assessment of structures that offer protection against the consequences of landslides.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 2 .
Fig. 2. Sequence of MPM algorithm steps occurring in each time step: (a) mapping and integration, (b) solution and material point location update, and (c) mesh reset.

Fig. 3 .
Fig. 3. Schematic of two bodies represented by material points and the variables used to detect contact calculated on the background grid, i.e. the body velocities, combined velocity and nodal normals.

Fig. 9 .
Fig. 9. a) Initial contact loads, b) contact loads at instant of collision (t = 0.8 s), c) distributed contact loads (t = 1 s), and d) contact loads at the end of the simulation.Note the different scaling in each sub-figure.

Fig. 11 .
Fig. 11.a) Slope dimensions and boundary conditions, b) construction stages,and c) building and wall foundations.Note the unequal scales to enable better visualisation; dimensions in meters.

Fig. 12 .
Fig. 12. a) Initial slope, b) after two excavation steps, c) after four excavation steps, and d) excavation finalized, including installation of retaining wall (depicted by black material points).