A MPM framework for large deformation seismic response analysis

: Landslides are often triggered by earthquakes and can cause immense damage due to large mass movements. To model such large-deformation events, the material point method (MPM) has become increasingly popular in recent years. A limitation of existing MPM implementations is the lack of appropriate boundary conditions to perform seismic response analysis of slopes. In this article, an extension to the basic MPM framework is proposed for simulating the seismic triggering and subsequent collapse of slopes within a single analysis step. Original implementations of a compliant base boundary and free- ﬁ eld boundary conditions in the MPM framework are presented, enabling the application of input ground motions while accounting for the absorption of outgoing waves and the free-ground movement at the lateral boundaries. An exam- ple slope is analysed to illustrate the proposed procedure and to benchmark it against the results obtained using an independent simulation technique, based on a three-step ﬁ nite element (FE) analysis. The comparison generally shows a good agreement of the results obtained from the two independent procedures and highlights advantages of the presented “ all-in-one ” MPM approach, in particular for long duration strong motions.


Introduction
Earthquake shaking can trigger the failure of slopes, which are stable under purely gravitational loading. The damage of such seismically triggered landslides can be large and, for many earthquake events, exceed the damage from all other seismic hazards combined (Kramer 1996). Numerous seismically triggered landslides highlight the severity of such catastrophic events (Havenith et al. 2016;Keefer 2002;Rodríguez et al. 1999). Seismic slope stability analysis is therefore a key component for earthquake hazard assessments. A number of techniques and methods have been developed for this purpose (see Section 1.1). A comprehensive risk assessment, however, includes not only an analysis of the slope stability during seismic shaking, but also of the further evolution of a potential slope collapse. The simulation of the movement of unstable soil masses often imposes a large deformation problem and hence requires different analysis techniques (see Section 1.2). Due to limitations in existing numerical approaches, these two stages are typically analysed separately, using different methodologies. In reality, however, the evolution of seismically triggered landslides is a continuous process, which cannot be strictly separated into separate steps. In this article, this issue is addressed and a methodology is presented, for the simulation of seismically triggered landslides as a continuous process using a single analysis step.

Modelling seismic response of slopes
Assessing slope stability under seismic loading requires an adequate modelling of the effect of the seismic action on the soil structure. Numerous techniques of varying sophistication have been developed in the past for such applications, ranging from pseudo-static limiting equilibrium approaches for stability analyses, to Newmark's sliding block analysis (Newmark 1965) for computing permanent displacements, to stress-deformation analyses for assessing the dynamic response of slopes (e.g., Hashash and Groholski 2010). The latter have become the state-of-the-art technique for deterministic seismic slope stability analysis. They can be divided into two main categoriesfrequency domain and time domain analyses (Hashash and Groholski 2010). For seismic slope stability analyses, where the shear stress reaches or even exceeds the shear resistance of the soil, the latter is often the preferred method (Kramer 1996) as it enables the use of nonlinear, elastoplastic constitutive models. This requires solving the dynamic equation of motion, for which the finite element method (FEM) and the material point method (MPM) are commonly used.
To simulate the propagation of seismic waves adequately, appropriate boundary conditions have to be applied to avoid the reflection and trapping of waves within the model domain. Special boundary conditions have been developed within the FEM in the past, for both the lateral boundaries and the base boundary of the model.

Lateral boundaries
At the lateral model boundaries, outgoing waves should be able to leave the model without reflection. At the same time, the freefield movement of the soil should not be restricted. One possible way to satisfy these requirements is by applying the free-field boundary method (Wolf 1989;Zienkiewicz et al. 1989). Free-field columns are included on both sides of the main model (see Fig. 1). The movement of the columns corresponds to the far-field response of the horizontally layered ground and is unaffected by the main model, but not the other way around. To absorb outgoing waves, the main model is connected to the free-field columns with dashpot elements (Nielsen 2006).

Base boundary
Typically, the input ground motion is applied at the base of the model domain to model the incoming waves (Fig. 1). The simplest way to do this is to prescribe the input acceleration or velocity directly. This, however, is equivalent to a rigid base, where downwards propagating waves are reflected. In many cases, this is not justified. This limitation can be overcome by applying a "compliant base boundary condition" (Zienkiewicz et al. 1989). It allows for the direct application of the input ground motion without prescribing the movement and causing a reflection of outgoing waves. First, viscous elements are connected to the base of the model; they are absorbing outgoing waves (Lysmer and Kuhlemeyer 1969). Second, instead of applying the acceleration time history, the input motion is transformed into a boundary traction. This traction time history is then applied at the base of the model and computed as ð1Þ where v su is the particle velocity of the upward propagating wave, i.e., half the "outcrop" motion (Mejia and Dawson 2006), r is the bulk density, and G is the shear modulus of the base material. The factor of two is required as half of the applied stress is absorbed by the viscous elements.

Solving large deformation problems in geotechnical engineering
With the methodology described above the seismic response of slopes can be adequately simulated. If the applied ground motion is sufficiently strong, it can initiate a collapse of the slope. In this case, the velocity of the unstable soil mass increases rapidly and the deformations in the model can become very large. Traditional geotechnical continuum analysis techniques are not well suited to analyse such large deformation problems. In recent years, however, a number of advanced techniques have been developed. A comprehensive review of different numerical approaches was provided by Soga et al. (2016). Two of the most frequently used approaches are the coupled Eulerian Lagrangian (CEL) FEM and the MPM.

Material point method (MPM)
The MPM is a hybrid method, where the material is represented as Lagrangian particles. The equations of motion are solved on a Eulerian grid. The Eulerian solution procedure allows for the material to undergo large deformations, while the Lagrangian material representation presents a convenient way of tracking material properties and constitutive state variables. MPM is a generalization of the particle in cell (PIC) and the fluid implicit particle (FLIP) method for solid mechanics, first proposed by Sulsky et al. (1994). Nowadays, numerous variations of the original algorithm have been introduced, tailored for a range of different applications (Bardenhagen and Kober 2004;Jiang et al. 2015;Steffen et al. 2008). In recent years the MPM has been successfully applied to model landslides (Andersen and Andersen 2010;Soga et al. 2016) and significant advances have been made, such as the development of hydromechanically coupled approaches (Abe et al. 2014;Bandara and Soga 2015;Zabala and Alonso 2011). Also seismic slope failure (Bhandari et al. 2016;Ering and Sivakumar Babu 2020) and landslides triggered by earthquakes (Alsardi and Yerro 2021;He et al. 2019) have been modelled with MPM. In these studies, the ground motion was prescribed as a velocity boundary condition and the lateral boundaries were modelled as either rigid or free moving, leading to a reflection of outgoing waves back into the model domain. Whereas this correctly represents the boundary conditions of the simulated shaking table experiments, it is a source of inaccuracy when modelling the seismic response of real slopes.

Coupled Eulerian Lagrangian (CEL) method
With the CEL finite element (FE) analysis technique a traditional Lagrangian phase, in which elements deform with the material, is followed by a Eulerian phase, during which elements with significant deformation are remeshed and material flow between neighbouring elements is computed (Dassault Systèmes Simulia 2012). The material is tracked as it moves through the Eulerian mesh by computing its volume fraction, which represents the proportion of an element filled with a specific material. This allows for multiple materials to move within the model domain. The methodology has been successfully applied in recent years to simulate the progressing collapse of both subaerial and subaqueous slopes (Dey et al. 2015(Dey et al. , 2016Stoecklin et al. 2021). Similar to the MPM, existing CEL codes do not include appropriate seismic boundary conditions, and hence cannot be applied directly for seismic response analyses.

Goal and objectives
As outlined above, different techniques have been developed in recent years to simulate and analyse different processes related to earthquake-induced landsliding. However, to simulate seismically triggered slope failures as one continuous process, a unifying methodology is required. In recent attempts to develop such a framework (e.g., Stoecklin and Puzrin 2020), the co-seismic and post-seismic slope collapse analysis were performed in separate analysis steps. The effects of this simplification, however, could not be assessed. In this article, an attempt is made to develop a unified approach for performing co-and post-seismic landslide analysis, based on the MPM method. The main challenge lies in the implementation of appropriate seismic boundary conditions within the MPM framework. The resulting procedure is applied to simulate the seismic triggering and subsequent collapsing of an example slope and the results are compared to the Stoecklin and Puzrin (2020) type analysis. It is shown that for individual stages, the two independent approaches provide very similar results, thus validating the developed methodology. Furthermore, the proposed MPM framework enables the application of ground motion loading throughout the post-failure stage, and for the first time an assessment of the effects of seismic shaking on the co-seismic landslide evolution.

General MPM framework
Rather than using existing codes, an in-house MPM framework has been implemented in C++. This provides the required flexibility to implement the seismic boundary conditions and the necessary adaptions for performing seismic response analyses. The general implementation closely follows Stomakhin et al. (2013) and Jiang et al. (2015), using explicit time integration. Subsequently, the basic methodology and main equations are outlined, to provide the necessary background information. This serves as a summary of the essential equations rather than a stepwise derivation starting from the weak form of the balance equations. For simplicity, the subsequent equations are presented for the plane strain case where the deformations are restricted to two dimensions. The formulation, however, can be readily extended to the three-dimensional (3D) case.

General procedure
The concept of MPM is to use particles (so-called material points) to discretize the material and track mass, momentum, deformation, and constitutive state variables. To solve Newton's law of motion, a Eulerian background grid is introduced that allows for the computation of the derivatives needed for the stress-based force evaluation. Initially, each material point is assigned a position x p = (x p , y p ), volume V p , and mass m p . All the other quantities (i.e., velocity v p , deformation gradient F p , and Cauchy stress tensor r p ) are set to zero. Following the initialization, the motion of the material points is computed, based on the basic four step MPM algorithm (see Fig. 2).

Interpolation functions
Interpolation functions are defined over the nodes of the background Eulerian grid. The interpolation function at grid node i is denoted with N i (x) and is evaluated at the material point position x p = (x p , y p ). Often the more compact notation N i (x p ) = w ip is used, referring to the weight of material point p associated with the grid node i. As proposed by Steffen et al. (2008), a dyadic product of one-dimensional interpolation functions was used where h is the grid spacing and x i = (x i , y i ) is the grid node position.
To avoid the common problem of so-called grid-crossing errors, quadratic or cubic B-splines should be used. The application of Bsplines offers a straightforward extension of standard MPM and greatly improves the accuracy of the solution (Steffen et al. 2008;Tielen et al. 2017). Moreover, it allows for an efficient parallel computing implementation. Due to the longer span of the B-splines, special considerations have to be taken regarding boundary conditions. However, in the presented framework only rectangular truncated models are considered and therefore the concept of mirror particles (Schulz and Sutmann 2019) can be easily applied. Despite being computationally more expensive, cubic B-splines are applied in this work, as they are numerically more stable. Especially in cases where strain localization is expected to occur, this is of particular importance. The cubic basis function is defined as NðxÞ ¼ In addition, for computing the internal force vector and the velocity gradient, the gradient of the interpolation functions evaluated at the material point position rN i ðx p Þ ¼ rw ip is needed where N 0 (x) is the derivative of the basis function N(x). In contrast to the FEM, the interpolations functions and its gradients have to be recalculated in each time step, as the material points are moving relative to the grid. Therefore, in the following derivations the index n referring to time step n is added to both the weight (w n ip ) and the gradient (rw n ip ).

Particle to grid
In a first step, mass and linear momentum are transferred from the material points to the grid nodes using the weight w n ip at time step n ð5Þ where m p is the mass and v n p is the velocity vector of material point p. In the case where the affine particle in cell (APIC) method is used, the transfer of linear moment involves an additional term to preserve affine velocity fields and therefore conserve the angular momentum (Jiang et al. 2015) ð7Þ where B n p is a matrix stored at each material point and is defined in Section 2.5. D n p ¼ 1=3h 2 I for cubic interpolation functions. In geotechnical problems, large rotations are normally associated with plastic yielding and failure. These phenomena introduce considerable dissipation and therefore the loss of angular momentum is of minor importance. To be fully consistent, the APIC transfer should be adopted.
Furthermore, the internal force vector at the grid nodes can be computed as where V n p is the volume and r n p the Cauchy stress tensor at material point p. External forces (body forces and tractions) can be introduced at the material points by an external force vector f p . The external force vector is transferred from the material point to the grid as where it must be distinguished between body forces b p and traction forces t p . Body forces can easily be defined as a function of the volume or mass of the material point, on the one hand. For the definition of the tractions, on the other hand, the corresponding surface area is needed for distributed forces; it has to be recomputed in each time step in case it changes. It should be noted that in general cases, the direction of the surface traction can also change in each time increment. Tractions can also be prescribed at the grid and are considered separately in the following section as they do not have to be transferred.

Grid velocity update
Following the transfer described above from each material point inside the domain, calculations on the grid are performed. The resulting force at the grid nodes is computed as where f n bc;i denotes possible traction boundary conditions prescribed directly on the grid. For quasi-static analyses (e.g., computation of the initial stress field) so-called local damping can be used (Cundall 2002). Therefore, an additional damping force f n d;i proportional to the resulting grid force (or out-of-balance force) f n i and in the opposite direction of the velocity is introduced as where each component is evaluated separately and b is a dimensionless damping factor. For the seismic response analysis, where a compliant base and far-field boundary conditions are used, the application of damping is crucial for the computation of the initial stress conditions under static conditions. This ensures that velocities remain under a certain numerical threshold value, preventing the model from drifting away during the seismic analysis after replacing the kinematic boundary conditions with tractions. In the case of an implicit solver, the equilibrium equation could be solved instead of the equation of motion and no numerical damping would be required. Applying explicit time integration, the linear momentum can be updated according to At this point, kinematic boundary conditions on the grid can be applied. The corresponding nodal values of a kinematic boundary cannot simply be overwritten by its desired value, because the B-spline basis functions range outside of the domain and hence do not maintain unity inside the domain (Steffen et al. 2008). To overcome this problem, the method of mirrored particles has been applied in this work (Schulz and Sutmann 2019). Material points are reflected across the boundary in different ways depending on the type of boundary condition (slip or no-slip). These additional material points are not included in the calculation and serve purely for the formulation of the boundary conditions. Further details can be found in Ding et al. (2020).
Finally, the grid velocities can be computed as 2.5. Grid to particle and particle update After updating the grid velocities, they are transferred back to the material points. For a PIC or APIC transfer scheme the material point velocity is given by whereas for a FLIP scheme only the velocity increment is transferred In the case where APIC is used, the B-matrix has to be updated to Not only velocities, but also the velocity gradient L nþ1 p and deformation gradient F nþ1 p have to be updated is called the incremental deformation gradient. The updated deformation gradient also serves to update the volume of each material point The material point position is updated as which can be simplified for a PIC or APIC scheme as Finally, the stress tensor and potential internal variables are updated. Constitutive laws are implemented analogous to the FEM, where a relation has to be defined between the deformation gradient F n p (or the velocity gradient L n p in the case of viscosity) and the Cauchy stress tensor r n p . Depending on the constitutive model, this might include some additional internal variables such as plastic strain and other state variables. However, it must be emphasized that MPM is a large strain solution method and hence the choice of stress-and strain measures should be assessed carefully in terms of objectivity. For geotechnical engineering applications of the MPM, small strain constitutive models are often extended by using the Jaumann rate of Cauchy stress where _ r is the rate of Cauchy stress and the spin tensor W = skew[L] is the skew symmetric part of the velocity gradient. Applying an explicit numerical integration, the stress update can be written in incremental form as The stress updated includes a rotational component due to rigid body motion and a change in stress due to straining. The latter is described by the constitutive law in form of a relationship between the Jaumann rate of Cauchy stress and the rate of deformation tensor D = sym[L] Often, the stress update is implemented in a slightly different way. The current stress is rotated according to r n rot ¼ r n 1 Dt W n11 r n 2 r n W n11 ð Þ and subsequently updated using well-established stress integration algorithms (e.g., elastoplastic return mapping) It must be pointed out that eq. 22 lacks incremental objectivity (Hughes and Winget 1980) and is only justified if displacement increments are sufficiently small, which is true for the explicit time integration in the MPM algorithm. In this work the incremental objective stress integration after Hughes and Winget (1980) was used. Not only will this guarantee incremental objectivity, but it will also allow for a better comparison with results of the FE analyses, as the same algorithm is implemented in FE code ABAQUS (Dassault Systèmes Simulia 2012). Equation 22 is therefore replaced by where the rotation tensor is given by The application of small strain plasticity models in a finite strain framework comes with certain drawbacks. It has been shown that this can lead to different inconsistencies, such as spurious stress oscillations or improper energy dissipation (Bažant et al. 2012). However, it should be emphasized that these inconsistencies remain negligible for most applications in geotechnical engineering where elastic strains remain rather small (Simo and Pister 1984). The aim of this work is to present a basic methodology to perform seismic response analyses with the MPM, rather than discussing the appropriate choice of constitutive law and its implementation. Hence, this topic is not discussed any further here. For the elastoplastic seismic slope analysis performed in this study, a von Mises constitutive model with isotropic softening was used. The model was implemented using an implicit elastoplastic return mapping algorithm, which is consistent with the implementation in the FE code ABAQUS (Dassault Systèmes Simulia 2012). Further information regarding the implementation can be found in Simo and Hughes (1998).

Extension of MPM framework for seismic analysis
In this section, the extensions of the general MPM framework to perform the seismic response analysis are outlined and the implementation of seismic boundary conditions are described.

General modelling procedure
The analysis is commenced in a static step to compute the initial static stress field, followed by a dynamic analysis step. Although the term static is used here, the equation of motion is solved in this step (eq. 12) rather than any equilibrium equation. Static conditions are ensured by smoothly ramping up gravity forces over a sufficiently large time period and using local damping according to eq. 11. For the static analysis kinematic boundary conditions are applied. At the bottom boundary a no-slip condition is used whereas at the lateral boundaries the material points are allowed to move freely in the vertical direction using a slip condition (see Fig. 3a). Both the static and the dynamic analysis steps are performed on the same model. Only the boundary conditions are changed in between steps. The material points can be carried forward with all the stored information of their constitutive state variables. To ensure a correct treatment of the boundary conditions, it is crucial to align the external model boundaries with the background grid. For the dynamic analysis, two free-field columns are included in the model to simulate the far-field movement (see Fig. 3b). The base is assumed to behave elastically, which is a prerequisite for the ground motion to be prescribed as stress-time history according to eq. 1. Hence, the model truncation boundary has to be chosen sufficiently far below the surface where plastic deformations remain negligible. For the soil above the base, any suitable constitutive model can be used to represent the soil behaviour.

Base boundary
As pointed out in Section 2.4, the MPM requires a different treatment of kinematic boundaries compared to FE analyses. In this work, the method of mirrored particles was applied (Schulz and Sutmann 2019). As shown in Fig. 4a, the grid is extended with additional nodes and particles are mirrored across the boundary to mimic no-slip boundary conditions. However, the mirrored particles are not explicitly included in the model and only the boundary nodes need a special treatment.
For the seismic model the concept of wave-absorbing viscous elements (Lysmer and Kuhlemeyer 1969) was adapted to the MPM framework. Therefore, the kinematic boundary condition used for the static analysis is replaced by a traction boundary condition and absorbing elements are connected to the model boundary. In previous studies, the latter were directly connected to the boundary grid nodes (Al-Kafaji 2013; Shen and Zhen 2005). Due to the use of B-spline rather than bilinear shape functions in this study, the same procedure cannot be applied, as the unity condition is not fulfilled at the grid nodes. Only rectangular boundaries are used, which are aligned with the material points. For this layout, a novel approach is proposed here. Dashpots are directly applied to the bottom row of material points in the form of a traction force t p (see Fig. 4b). Similarly, the input motion is applied to the same row of material points as a traction force time history. Hence, the traction for a bottom material point can be expressed as the sum of three components: (i) a viscous surface traction representing the dashpot element, (ii) a reaction force derived from the static analysis, and (iii) a traction time history representing the applied input ground motion where d p is the surface area related to the material point ¼ h= ffiffiffiffiffiffiffiffi n MP p À Á (see Fig. 4b), which can be easily determined for a regular distribution of material points along the bottom boundary (h: grid spacing; n MP : number of material points inside a grid cell). The material point velocities in the xand y-directions are v px and v py , respectively, whereas the corresponding viscosities of the dashpots h s = r c s and h p = r c p are given by the shear-and pressure-wave speed of the base material, c s ¼ ffiffiffiffiffiffiffiffiffi G=r p and c p ¼ ffiffiffiffiffiffiffiffiffiffi M=r p (M is the elastic P-wave modulus). The input ground motion is represented by the particle velocity of the upwards propagating shear and pressure wave (v su and v pu ). Finally, to ensure the static equilibrium of the model, the static shear stress t p,s and normal stress  s p,s of the corresponding material point are applied as an initial condition and kept constant through the analysis. For the sake of simplicity, the index of the time step n is omitted.
It should be pointed out that switching from a kinematic boundary condition (using mirrored particles) to a traction boundary, where forces are applied to material points, does not usually represent a perfect transfer as this traction is mapped to several grid nodes in the particle to grid phase. This could lead to a slight lack of balance static equilibrium at the lower boundary, causing a slow drift of the model. Nevertheless, the influence of these redistributions is usually negligible and the problem of drifting can be solved by replacing the dashpots by Kelvin-Voigt elements (Al-Kafaji 2013). In this work, however, a stricter approach is proposed, where the static reaction forces are directly applied to the nodes as the equation of motion is solved on the grid. After the stress initialization in the static analysis, the particle to grid transfer is performed to get internal and external nodal forces, but this time without the mirrored material points. Hence, at the boundary nodes, equilibrium of force is not strictly fulfilled anymore. The resulting out-of-balance forces are then applied as a traction boundary f i = (f ix , f iy ) on the grid for the dynamic analysis with the opposite sign (see Fig. 4b).

Lateral boundaries
Similar to the base boundary, the lateral boundaries are also modelled by the mirrored particle approach for the static analysis. In contrast to the base, a slip boundary is used to allow for vertical deformations. Therefore, only the horizontal components of the nodal forces at the boundary nodes (see Fig. 5) need a special treatment (Schulz and Sutmann 2019). For the free-field columns (Wolf 1989;Zienkiewicz et al. 1989) separate models are created because a different grid is needed. Their movement is not influenced by the main model and hence can either be precomputed or computed in parallel to the main calculation. The free-field columns represent infinitely flat ground, which can be model by using so-called periodic boundary conditions. This is achieved by giving the node at the right boundary the same number as the corresponding node on the left boundary and hence both nodes refer to the same memory. The column is modelled with a width of four grid cells, which is exactly the span of the cubic shape function. In theory, a width of only one cell would be sufficient, but this would lead to a rather cumbersome implementation. Following Nielsen (2006), the freefield columns are connected to the main model using dashpot elements (see Fig. 5). The connection is imposed directly on the boundary material points in the form of a traction force similarly to the bottom boundary. However, the viscous component is defined as a function of the relative motion between the corresponding material points of the free-field (index fp) and the main model (index p). The traction follows as where d p is the surface area related to the material point, v fpx and v fpy are the velocity components of the material point in the freefield column, and v px and v py the corresponding material point velocity in the main model (remark: the time step index n is omitted for simplicity). The viscosities of the dashpots h s = r c s and h p = r c p are given by the shear-and pressure-wave speed of the material, which might be depth dependent. Similar to the base boundary, static equilibrium is ensured by the static normal stress component s p,s . Moreover, the dynamic traction due to the seismic excitation of the free-field column has to be applied on the main model (Nielsen 2006) and is retrieved from the dynamic stress tensor of the corresponding free-field material point where r fp is the actual stress tensor and r fp,s is the stress tensor at the end of the static analysis of the free-field column. The components of the dynamic surface tractions s fp,d and t fp,d follow by multiplying the dynamic stress tensor with the outer normal vector of the lateral surface.
Analogously to the bottom boundary, the static support is implemented in a more stringent approach by a traction force f ix , which is directly applied to the boundary nodes (see Fig. 5). It should be emphasized that this approach is limited to models where the horizontal deflections of the lateral material points are small compared to the grid size. Otherwise, material points might fall out of the influence zone of the boundary nodes defined by the shape function span and static support is not ensured anymore. In this case, static tractions should be applied on the boundary material points instead.
In contrast to the bottom boundary, the computation of the surface area related to the material point d p in eq. 30 is not straightforward and has to be performed for each boundary material point Fig. 5. Illustration of lateral boundary condition including free-field column (left). (v fp : velocity vector of corresponding material point in free-field; f ix : out-of-balance force as traction on grid; t px , f py : traction on material point due to dashpot elements and dynamic surface traction).
separately. As the deformation gradient is stored for each material point, d p can be computed as where l is the stretch of a normalized line element in the initial configuration and d 0 p is the initial surface related to the material point before any deformation occurred (which is before the static analysis). The stretch is given by the normalized line element s = (0 1 0) T and the Right-Green deformation tensor C = F T F. In case large deformations at the boundary material points are expected during the seismic analysis, d p could be updated in each time increment. However, for most applications this is rather unlikely. Although in the current work this would not be necessary, the stretch is updated in each time increment for consistency.
Material points and grid nodes at the corners (i.e., where the lateral and the base boundaries intersect) do not need to be treated specially. The contributions of both boundaries can simply be superimposed. Special consideration is only required for the mirrored particle approach in the case of a parallel computing implementation where possible race conditions must be carefully checked and avoided.

FE-Eulerian methodology
To benchmark the results obtained using the MPM procedure described above, the slope failure process is simulated with an alternative methodology in a sequence of three steps: (i) a static analysis step to compute the pre-failure static stress field within the slope, (ii) a dynamic analysis step to simulate earthquake events followed by (iii) a post-failure analysis step to compute the motion of the collapsing slope. Each analysis step is based on a different approach because the seismic boundaries are not available within the applied CEL framework. They are connected by prescribing the results from the preceding step as initial conditions for the subsequent step (see Fig. 6). The framework was developed within the ABAQUS computing environment (Dassault Systèmes Simulia 2012). It has been successfully applied in previous studies to analyse the behaviour of earthquake-induced subaqueous landslides (Stoecklin and Puzrin 2020;Stoecklin et al. 2021). In this section, the three main analysis steps are described briefly.

Static analysis
In this first step, the pre-failure stress field within the slope is computed under static conditions, using a total stress-based, implicit, plane-strain Lagrangian FE analysis (Fig. 6a). The same kinematic boundary conditions are applied as in the MPM simulation: no-slip condition at the bottom boundary, where the displacements in the xand y-directions are fixed, and slip conditions at the lateral boundaries, where the displacements are only fixed in the y-direction.
The resulting stress values are then transferred for each integration point to the subsequent seismic analysis step.

Seismic analysis
In this second analysis step, the impact of a seismic event is simulated by subjecting the slope to an earthquake ground motion (Fig. 6b). A dynamic, implicit, total stress-based, plane-strain FE Lagrangian analysis procedure was employed. At the lateral boundaries, the free-field boundary method was applied to avoid reflection of stress waves without restricting the free field movement (see Section 1.1). The free-field columns were a large out-of-plane thickness to ensure that they move nearly independently from the main model, which was connected to the lateral free-field columns with dashpot elements. At the base of the model, viscous elements were connected to the model to achieve the effect of a compliant boundary condition. To ensure stress equilibrium under static conditions, the boundary tractions obtained from the static analysis are prescribed as initial stress conditions at the lower and lateral model boundaries. The ground motion was applied as a traction time history at the nodes where the viscous elements are connected to the main model (Mejia and Dawson 2006). The modelling procedure was described in more detail by (Stoecklin and Puzrin 2018).
If the analysed slope becomes unstable as a result of the applied ground motion and the degradation of the shear resistance, the displacements in the unstable soil mass accumulate at an increasing rate. Hence, the solution loses accuracy due to excessive mesh distortion and the analysis has to be terminated. To simulate the motion of the collapsing slope past this point, the results are transferred to the post-failure analysis.

Post-failure analysis
In this final step, the motion of the collapsing slope is simulated until it reaches a static equilibrium once again. The velocity, strength, and stress fields are mapped to the post-failure model, serving as starting conditions (Fig. 6b). The analysis is based on the CEL FE approach, allowing for the materials in the model to undergo extreme straining without suffering from excessive mesh distortion (see Section 1.2). A void space was included above the soil layer, facilitating the free movement of the material within the specified Eulerian domain. The procedure was described in more detail by Stoecklin et al. (2021).
While the method enables the computation of the large displacements within the collapsing soil mass, the absorbing boundary conditions cannot be implemented easily. To avoid trapping seismic stress waves in the model, the input ground motion was not applied anymore during this final step. At this stage, however, the movement of the unstable soil mass is self-driven and, in many cases, only marginally affected by the ground motion. In some Fig. 6. Illustration of analysis procedure in three main steps: (a) computation of pre-failure stress-field within slope followed by (b) simulation of earthquake event and (c) post-failure evolution of collapsing slope.
cases, however, this simplification can lead to an underestimation of the predicted displacements (see Section 5.4).

Analysis of an example case and comparison between different procedures
To compare the results obtained with the MPM and the CEL methodologies, an example case is analysed in this section. The comparison of the results serves as a benchmark for the developed MPM extension and illustrates how the procedure can be applied. Following a description of the example case, the MPM versus CEL comparison is presented, providing validation for both approaches.

Description of example case
The example case comprises a 10 m thick soil layer on top of an elastic base (Fig. 7a). The slope curvature was chosen to follow a Gaussian shape (Adams and Schlager 2000), with a maximum inclination of about 22° (Fig. 7b). The applied soil parameters are listed in Table 1. The slope was subjected to two input ground motions: (i) a motion recorded from the Loma Prieta event in 1989 (RSN 769, H2 direction) and (ii) a motion recorded from the Imperial Valley event in 1979 (RSN 165, H2 direction). The time histories were retrieved from the Pacific Earthquake Engineering Research (PEER) strong motion database (Ancheta et al. 2014). The recordings were chosen arbitrarily for the purpose of demonstration of the procedure.
For the MPM model a square grid with element size h = 0.25 m and a regular distribution of 3 Â 3 material points has been chosen to reasonably represent the slope curvature and to accurately model higher frequencies in the seismic analysis (see Fig. 8a). For more complex slope geometries, a finer grid may be chosen or the concept of isoparametric elements could be introduced (Tjung et al. 2020). A structural mesh with quadrilateral four node FEs (Q4) that follows the slope curvature is used for the CEL model (see Fig. 8b). The mesh size h el is chosen identical to the MPM.
Undrained soil behaviour was assumed, represented by a von Mises constitutive model with isotropic softening in combination with an isotropic, linear elastic model. The stressÀstrain response is linear and elastic up to the peak shear strength, followed by linear softening until the value of the remoulded shear strength is reached and the strength remains constant (Fig. 9).
It should be mentioned that the presented example case is a theoretical example and does not represent a particular slope.
The slope shape and curvature as well as the used soil parameters were chosen to represent typical values (see Table 1).
The residual shear strength is defined by the peak shear strength q p and the sensitivity of the material S.
As the specified strain-softening constitutive behaviour can lead to mesh-dependent results, a smeared crack approach was employed as a regularization technique (Rots et al. 1985). This approach is often employed for MPM simulations featuring strainsoftening materials (e.g., Soga et al. 2016;Yerro 2015) and provides the advantage that the numerical shear band thickness does not have to match the in situ zone of intensive shearing. Therefore, a much coarser mesh can be used. It should, however, be assured that the mesh is fine enough to allow the shear bands to develop and propagate. A shear displacement in simple shear loading was used rather than a shear strain to define the post-peak softening curve. Assuming simple shear conditions, the plastic shear strain at which the material is fully softened was scaled as a function of the shear band thickness  where d r is the displacement at which the material is fully softened and h shear the shear band thickness in the computational model (Fig. 9b). A typical value of d r = 0.2 m was chosen for the calibration (Skempton 1985). For the FE analysis it was assumed that the shear band develops within a single element. Hence, the element size is equal to the shear band thickness (h shear = h el ). Due to the longer span of the B-Spline in the MPM analysis, the scaling parameter h shear was calibrated prior to running the analysis, based on a numerical simple shear element test. The calibration showed that this leads to a shear band thickness of about twice the grid cell size (h shear ffi 2h), which corresponds to exactly half the span of the cubic B-spline. However, this topic is not discussed in detail here and will be the subject of future investigations.

Results of seismic response analysis
In a first step, the implementation of the seismic boundary conditions was benchmarked. For this purpose, the resulting motion within the slope was compared for the MPM and the FE methodologies. The comparison shows whether (i) the input ground motion is prescribed correctly, (ii) the wave propagation is simulated correctly, and (iii) whether outgoing waves are absorbed at the model boundaries. To exclude effects of strain localization and scaling, both the soil layer and the base were assigned a purely isotropic, linear elastic behaviour.
A comparison of the computed motion is shown in Fig. 10 for two arbitrarily chosen points within the model domaina point at the ground surface and a point at the boundary between the soil layer and the base (see Fig. 7a). The comparison shows that the solutions obtained with the MPM and the FEM analysis approach are nearly identical. The same holds for any other point within the model and other output quantities, such as acceleration or deformations. The results therefore validate the implementation of the seismic boundary conditions in the MPM framework.

Results of slope failure analysis
Having validated the implementation of the seismic boundary conditions, the process of failure initiation and propagation by seismic ground motion was analysed next. This comparison provides an indication whether (i) the degradation of the shear resistance by cyclic loading, (ii) the initiation of a failure mechanism, and (iii) its propagation along localized shear bands are captured adequately in the simulations. The soil material was assigned a strain-softening constitutive behaviour (see Section 5.1), whereas the base was assumed to behave elastically.
The resulting displacements and distribution of plastic strain within the soil layer are shown in Fig. 11. It can be observed that the slope failure is initiated nearly at the same time in both simulations. Following failure initiation, a very similar failure mechanism develops in both simulations. As the slope collapse evolves further, some differences emerge between the two solutions in terms of the development of secondary shear bands. A more distinct difference can be identified at the newly formed toe. In MPM the shear band localizes less and a steeper toe is formed, whereas in the CEL simulation the toe collapses and results in a milder gradient. These differences are most likely a result of the often reported sticky behaviour of MPM (Huang et al. 2011;Soga et al. 2016) as material points remain numerically in contact when they have a node in common. The long span of cubic B-splines might lead to an even more pronounced phenomenon. However, this does not affect the main results of this analysis. Modelling strain localization problems using MPM in general has room for improvement and is the subject of ongoing research studies.   9. Illustration of (a) assumed constitutive behaviour and (b) adopted simplified scaling method.
Nonetheless, the geometries of the final deposit are remarkably similar. The comparison thus improves the confidence in the obtained solution and validates the two independent methodologies against each other. For consistency, the ground motion was applied only until the failure became self-driven for both simulations. This simplification can be avoided with the MPM methodology, which is discussed in the next section.

Advantages of MPM procedure
As shown in the previous sections, the developed MPM framework provides very similar results to the FE analysis procedure. As the MPM analysis is computationally more expensive, the FE analysis procedure remains, in the authors' opinion, the preferred method for seismic response analysis of problems where deformations remain limited. However, for applications where earthquake loading leads to the initiation of instabilities and large deformations within the ground, the MPM approach offers a significant advantage. It offers an "all-in-one" approach, allowing for the simulation of seismic ground excitation and the evolution of large deformations within the ground at the same time. For some applications it can therefore provide more accurate results.
To demonstrate this effect, the results of two simulations were comparedone where the ground motion was applied throughout the entire duration of the analysis and one where the ground motion was applied until the slope failure becomes self-driven. When applying the relatively short Loma Prieta (1989) ground motion the results are nearly indistinguishable, showing that the ground motion does not influence the movement of the slide significantly upon failure initiation. However, as shown in Fig. 12, applying the Imperial Valley (1979) ground motion, which has a much longer strong motion, leads to significantly different results. It can be observed that, early on, a similar failure mechanism is triggered. However, with a continuing ground excitation, the failure mechanism propagates further down-and uphill into stable parts of the slope, resulting in a considerably larger landslide. Therefore, applying the unified MPM procedure can provide more reliable risk assessment results for certain cases.

Limitations
The presented approach has a number of limitations, which should be refined in future studies. For the validation of the procedure, both geometry and the employed constitutive models were kept as simple as possible and are therefore not necessarily appropriate to Fig. 10. Comparison of resulting horizontal velocity v x for point at ground surface (top) and boundary between soil layer and base (middle) using MPM and FEM analysis procedure. The Loma Prieta (RSN 769) recording was applied as an input ground motion (bottom). [Colour online.] represent real-life slopes. Nevertheless, the methodology can readily be applied to include more complex geometries and more advanced constitutive models.
Furthermore, other effects, which have been shown to influence the failure behaviour of slopes, such as the effect of heat and excess pore-water pressure generation (e.g., Abe et al. 2014;Bandara and Soga 2015;Zabala and Alonso 2011) were out of the scope of this study. The presented framework could be applied in combination with such techniques in future investigation to investigate for these effects.
Another potential shortcoming is that the ground motion was applied as a stress history distributed uniformly along the entire base of the model, implying that the earthquake signal arrives simultaneously at each point of the base. A way to overcome this shortcoming in future investigations would be to apply a nonuniform input signal from a larger seismological model (Bielak 2003).
It should be mentioned that using the compliant base boundary and free-field columns lead to an unsupported model during the seismic analysis, regardless of the numerical method. This can cause a drifting of the model. This problem might be more pronounced in MPM and must be carefully considered when evaluating permanent displacements. It is crucial to apply local damping for the stress initialization to keep the initial material point velocities small. In combination with prescribing the static support directly on grid nodes during the seismic analysis, this often helps to reduce the model drift to an acceptable level. For the presented example analysis, the associated stress redistributions at the boundaries and the effect of the model drift are negligible. In cases where the simulation is continued after the seismic event for a longer time span, kinematic boundary conditions could be reapplied or the dashpots exchanged with KelvinÀVoigt elements.

Conclusions
Dynamic response analysis in the time domain is considered the state-of-the-art technique for assessing the stability of slopes under seismic loading. For a comprehensive risk assessment, however, not only the stability of a slope during earthquake loading is of interest, but also the subsequent mass movement. Traditional approaches to seismic response analysis, such as the standard finite element method (FEM), solve the equation of motion using Lagrangian framework and fail to provide accurate solutions in case large deformations occur, due to excessive mesh distortion.
The material point method (MPM), on the contrary, has become a successful tool to model such large deformation problems. However, so far, it has not been widely applied for seismic response analyses, partly due to the unavailability of appropriate boundary conditions. In this article, this gap is addressed and a procedure is Fig. 11. Comparison of evolution of plastic strain within example slope at different time points during analysis for both (a) MPM and (b) CEL simulation. The slope was subjected to the Loma Prieta (1989) ground motion until the failure became self-driven.
presented to model the entire process from initiation to post-failure evolution of seismically triggered landslides within a single analysis step.
For this purpose, an extension to the basic MPM framework is presented. Special considerations for the computation of initial stress conditions and the switch to dynamic boundary conditions are outlined and a formulation and implementation of suitable dynamic boundary conditions are presented. For the latter, boundary conditions that are established for dynamic FE analyses were adapted and implemented in the MPM framework, namely the compliant base boundary and the free-field boundary method.
A comparison between the computed seismic response of an example slope using the MPM and a traditional FEM procedure shows that the obtained results are nearly indistinguishable. Furthermore, the results of the co-and post-seismic analysis obtained from the two procedures compare remarkably well in terms of the time of triggering, the resulting failure mechanism, and the post-failure geomorphology for the analysed case. An important advantage of the proposed unified MPM approach is that the seismic shaking can be applied continuously throughout the analysis, even after the slope collapsed and large deformations have occurred. For the first time it is possible to quantify this phenomenon for a case of a long duration strong ground motion and to demonstrate that it can have a significant effect on the landslide evolution.