Modelling the pre-and post-failure behaviour of faulted rock slopes based on the particle finite element method with a damage mechanics model

We develop a novel numerical model incorporating a damage mechanics formulation into the particle finite element method for analysing the pre-and post-failure behaviour of faulted rock slopes. In this computational framework, the stress-driven strain localisation and damage evolution are modelled based on an isotropic damage model; discontinuity structures like fault zones are represented as thin continuum layers with equivalent mechanical properties; the particle finite element method is used to solve and track the large deformation of rock masses. In the paper, we first present the mathematical formulation of the proposed model in the context of the Hellinger-Reissner variational principle. We conduct a thorough validation of our model for simulating the damage of brittle materials against well-documented experimental datasets of different failure scenarios. We then apply the model to simulate the deformation and failure phenomena of faulted rock slopes including both the pre-failure progressive damage and post-failure transient runout, demonstrating the strong capability of our model in physically capturing the initiation, evolution, and consequence of catastrophic rock slope failure across multiple temporal scales. In addition, our simulations can realistically reproduce the slope displacement time series with the interplay between rock damage and fault reactivation explored. The present work has important implications for understanding the physical mechanisms that drive the progressive destabilisation and cata-strophic failure phenomena of rock slopes in nature.


Introduction
Landslide phenomena usually involve the movements of a mixture of rock, earth, and/or debris downslope (Cruden, 1991), which develop over time experiencing from initiation to propagation, and to final deposition (Leroueil, 2001).During this course, a landslide may exhibit different mechanical responses, among which the most critical one might be the formation of a fully developed slip surface (Hungr et al., 2014), often referred to as failure.Before the failure, the slope deformation driven by stress changes and irreversible damages tend to show a progressive destabilisation behaviour (Eberhardt, 2008;Leroueil, 2001).During and after the failure, the geomaterials that have undergone a certain degree of strength loss may exhibit a significantly increased mobility, leading to rapid inertia-dominated movements (Hungr et al., 2014).To analyse the complex behaviour of landslide materials during these different stages, various modelling approaches have been developed (Chen, 2007;Denlinger, 2004;Duncan, 1996;Griffiths and Lane, 1999;Savage and Hutter, 1989;Sloan, 2013).For the pre-failure stage, the state of failure has been determined from static or dynamic perspective: (i) the static stability is mostly analysed using the limit equilibrium method (Duncan, 1996), limit analysis method (Chen, 2007;Sloan, 2013), and strength reduction method (Dawson et al., 1999;Griffiths and Lane, 1999); (ii) the dynamic stability is used to analyse the behaviour of geomaterials under complex conditions (Duncan, 1996), solved using different advanced numerical techniques.As for the postfailure stage involving large deformations, the moving mass is usually governed by the simplified Navier-Stokes equations due to their flowlike characteristics (Denlinger, 2004;Savage and Hutter, 1989).These two stages of slope movements (i.e.pre-and post-failure stages) have long been treated separately in landslide analyses due to the lack of a unified computational framework for modelling the entire landslide process that occurs across multiple temporal scales (Eberhardt, 2008;Soga et al., 2016).
Even though these methods can describe both the formation of a continuous shear surface during the failure and the post-failure movements of shallow sliding in soils (Cascini et al., 2010;Cuomo et al., 2021;Hungr et al., 2014), the numerical models with elasto-plastic relationships may not be able to describe the complex failure mechanisms of rock slopes.Indeed, the failure processes of rock slopes are often strongly affected by the state of in situ stresses and the distribution of pre-existing discontinuities (Eberhardt et al., 2004;Einstein et al., 1983;Hoek and Bray, 1981).For geological structures such as faults, the sliding movements along major structures can be simulated using the elasto-plastic relationships (Eberhardt et al., 2004;Stead and Wolter, 2015).While for the brittle fracturing developed in rock bridges in between discontinuity structures, the numerical model should include a mixture of deformation/failure modes (e.g.tensile, shear, or and mixedmode failure) generated during the brittle deformation prior to macroscopic failure in rocks (McBeck et al., 2020;Stanchits et al., 2006).Such a brittle behaviour is commonly implicitly represented in continuumbased approaches (Jirásek and Bauer, 2012), or explicitly simulated in discontinuum-based approaches (Eberhardt et al., 2004;Munjiza et al., 1995;Scholtès and Donzé, 2012).In discontinuum modelling (such as the discrete element method), interaction forces have to be resolved between each discrete element pair (Scholtès and Donzé, 2012), which is computationally very expensive and sometimes difficult to fully reflect the macroscopic behaviour of landslides (Soga et al., 2016).Furthermore, the time-dependent slope destabilisation behaviour during the pre-failure stage often cannot be modelled in real time (over a timescale of e.g.months to years) using a discontinuum-based model due to its relatively small time step (at a timescale of e.g.microseconds) required in the explicit time integration scheme (Stead et al., 2006;Sun et al., 2022).
To develop a unified computational framework for modelling the entire rockslide process, we implement a damage mechanics model (Jirásek and Bauer, 2012) into a PFEM framework (Wang et al., 2021b) for simulating the elasto-brittle behaviour of rock materials during the pre-failure stage as well as handling post-failure large deformation analysis.Our framework is established on the basis of the finite element method (FEM), where the PFEM is a mixture of the particle finite element technique and the Lagrangian finite element analysis (Oñate et al., 2004).The PFEM combines the robustness and accuracy of FEM for simulating complex problems and the advantages of particle-based approaches for large deformation analysis.So far, the PFEM has been applied to study various challenging landslide problems, such as landslide-induced waves (Cremonesi et al., 2011;Franci et al., 2020;Salazar et al., 2015), landslide propagation (Franci et al., 2020;Zhang et al., 2015), retrogressive failure in sensitive clay (Zhang et al., 2020), earthquake-induced landslides (Wang et al., 2021a), submarine landslides (Zhang et al., 2019), among others.However, the possibility of the PFEM for analysing the pre-and post-failure processes of rock slopes has not been explored yet.To demonstrate the proposed computational framework, we apply the model to simulate the progressive failure of a generic rock slope.We explore the performance of our model in covering both pre-and post-failure stages operating at different timescales (e.g. over weeks to months for the former stage and seconds for the latter stage).We investigate the kinematic characteristics of the initiation of slope failure (failure pattern and displacement time series) and the dynamical evolutions post failure (transition process and mass transport velocity), with the interplay between rock damage and fault reactivation captured.The rest of the paper is organised as follows: in section 2, we describe the methodology of our work; in section 3, we introduce the rock slope model setup; in section 4, we present and interpret the numerical results; in section 5, we discuss the implications of our numerical simulation results, followed by conclusions drawn in section 6.

Methodology
In this section, we present the mathematical formulation, solution techniques and benchmark tests.

Governing equations
Consider a 2D domain V, delimited by a surface S. The set of governing equations for dynamic analysis is formulated as follows.
(a) The dynamic equilibrium equation: where σ is the stress, v is the velocity, b is the body force, ρ is the density, and v is the time derivative of v.
(b) The displacement-strain relation: where ε is the strain and u is the displacement.
(c) Boundary conditions: where N is the unit normal vector, t is the traction, u p is the prescribed displacement, and S t and S u are the surfaces where traction and displacement constraints are imposed, respectively.
(d) Elasto-plastic constitutive relationship: (5) where F is the yield function, ε e and ε p are the elastic and plastic strains, respectively, C is the elastic compliance matrix, λ is the plastic multiplier, and Ψ is the plastic potential.The associated flow rule assumes Ψ = F, and the so-called complementary condition for the incremental form of Eqs. ( 5) and ( 6) can be written as follows: (e) Elasto-damage constitutive relationship: where ω is a scalar damage variable varying from 0 (undamaged) to 1 (fully damaged).The incremental form of the formulation can be obtained based on the stress at the known (with a subscript 'n') and unknown (with a subscript 'n + 1') states: where α = 1/(1-ω) and we assume ω can only asymptotically approach 1 (the maximum value ω max < 1) to avoid singularity at ω = 1.
The difference between Eq. ( 10) and Eq. ( 9) leads to: (11) where Δσ = σ n+1 -σ n , Δα = α n+1 -α n and Δε = ε n+1 -ε n .The damage variable ω may be related to the so-called equivalent strain ε (Jirásek and Bauer, 2012).An internal variable κ is introduced to record the maximum level of this equivalent strain governed by the following loading-unloading conditions (Jirásek and Bauer, 2012): The equivalent strains for tensile (ε t ) and compressive (ε c ) conditions are respectively computed based on the smooth Rankine criterion (Jirásek and Bauer, 2012) as: where ||⋅|| is the norm operator, 〈⋅〉 are the Macaulay brackets representing the positive part, D e is the elastic stiffness tensor, and E is the Young's modulus.The compressive and tensile conditions are denoted as positive and negative, respectively, following the geomechanics sign convention.The evolution of the damage parameters (i.e.ω t and ω c for tensile and compressive scenarios, respectively) may be captured by an idealised elasto-brittle model (applicable for highly brittle materials like hard rock) (Lei et al., 2021;Tang et al., 2008): where ε t0 = − f t0 /E and ε c0 = f c0 /E are the limit elastic tensile and compressive strains, respectively; f t0 and f c0 are the tensile and compressive strengths, respectively; η is the residual strength ratio; κ t and κ c are the internal variables for tensile and compressive conditions, respectively.The value of damage (ω = 1-η) at the limit state (the first moment κ t < ε t0 or κ c > ε c0 ) is denoted as ω 0 , as shown in Fig. 1b.In addition, we have also implemented an elasto-brittle model with strain softening (often used for concretes) (Jirásek and Bauer, 2012): and for which the crack band approach (Jirásek and Bauer, 2012) may be used to constrain: where G I and G II are the energy release rates for mode I and mode II failures, respectively; h b is the size of localised crack band.The constitutive laws are illustrated in Fig. 1.

Space domain discretisation
By means of the standard finite element notation, we discretise the field variables using the mixed triangular element (Krabbenhoft et al., 2007) as shown in Fig. 2. The stress is interpolated using linear nodes, while the displacement and inertial force are interpolated using quadratic nodes, such that: (17) where the hat '^' represents variables at the mesh nodes, and N σ , N u and N r are the shape functions.The strain tensor ε can be expressed as: Even though the use of quadratic interpolation nodes may introduce additional degrees of freedom (compared to linear interpolation nodes), such a mixed triangular element has proven to be robust in overcoming the locking issue in elasto-plastic computations (Krabbenhoft et al., 2005;Zhang et al., 2019).In the future, the proposed damage mechanics model will be further implemented into a hydro-mechanical PFEM model (Wang et al., 2021c) for analysing the damage propagation  L. Wang and Q. Lei phenomenon in saturated rock slopes, where the same mixed triangular element performs well in simulating dynamic processes in saturated media.

Time domain discretisation
In the framework of the θ-method (Wood, 1990), the stress and velocity evolution may be solved in a time marching fashion as: such that the equilibrium equation Eq. ( 1) can be discretised as: where Δt is the time step and 0 ≤ θ 1 ≤ 1 and 0 ≤ θ 2 ≤ 1 are two numerical parameters.The time integration scheme adopted in our model is unconditionally stable for θ 1 ≥ 0.5 and θ 2 ≥ 0.5.By introducing a new variable r, Eq. ( 21) is rearranged as with a traction boundary condition The update of velocity is based on Eq. ( 20), such that the new velocity is given as.

Mixed variational principle
A mixed variational principle has been previously established for dynamic elastoplastic problems (Zhang et al., 2013), which reads.
The Lagrangian of the above min-max optimisation problem (Zhang et al., 2013) is.
To also incorporate the damage model, we propose the following Lagrangian that: The solution of problem ( 28) is a saddle point of its functional and the associated Karush-Kuhn-Tucker (KKT) conditions can be given as.

Solution techniques
In this part, we present the solution strategy using the second-order cone programming (SOCP) and the particle finite element technique for large deformation analysis.

Reformulation and computation of the optimisation problem
The min-max problem (35) can be reformulated as a standard SOCP problem as.
where x = (x 1 , x 2 , ..., x n ) T contains the field variables, p represents the vector of coefficients involved in the objective function, L and y are known matrix and vector for linear equality, and ℵ is the tensor product of second order cones.Here, the second-order cones can be: } .Following the procedures developed in previous works (Krabbenhoft et al., 2007;Wang et al., 2021b;Zhang et al., 2013), the min-max problem ( 35) can be reformulated to a minimisation problem as: in which s 1 and s 2 are unknown positive scalars.By introducing , the above problem is further expressed as: where the underlined terms are rotated quadratic cones.The yield function F shown in (38) should also be formulated as cones.Here, we consider the Mohr-Coulomb criterion for plane-strain problems, given as: where c is the cohesion and φ is the friction angle.The equivalent standard cone for each stress node can be represented by introducing a new set of variables ρ = (ρ 1 , ρ 2 , ρ 3 ) T and prescribing ρ = Hσ + d: One can see that the yield function ( 39) is exactly the quadratic cone . Note that for non-associated flow rule, the plastic potential Ψ is different from the yield function F: in which ψ is the dilation angle.The associated computational scheme can also be applied using an approximation form of F (F ≈ F * →∂F * /∂σ n+1 = ∂Ψ/∂σ n+1 ): where c is treated as a constant and updated at each time step using stresses at the known state: This approximation has been used in typical geotechnical problems with very good performance obtained (Krabbenhoft et al., 2012;Wang et al., 2018Wang et al., , 2021c;;Zhang et al., 2013).Then, we have the following minimisation problem: where N G denotes the number of stress nodes (Gauss points).At each time step, the values of α n+1 and Δα are updated by the damage variable ω, calculated from Eqs. ( 13) -( 16).In the model, the material is either elasto-damage or elasto-plastic.We prescribe F(σ n+ 1 ) < 0 for elastodamage material and a constant α = 1 for elasto-plastic material.The problem ( 44) is of the same form shown in (36), and can be solved using the interior-point method (Andersen et al., 2003).In this study, we use the modern optimisation engine MOSEK (MOSEK ApS, 2019), and its detailed implementation can be found in (Wang et al., 2021b).

Particle finite element technique
The particle finite element method (PFEM) combines the particle approach and the Lagrangian FEM by treating mesh nodes as free particles and solving the governing equations based on the Lagrangian FEM framework (Idelsohn et al., 2004).To achieve this, the alpha-shape method (Edelsbrunner and Mücke, 1994) is used for identifying the boundaries of computational domains, which are further used for mesh generation and Lagrangian FEM analysis.The PFEM and its variants (Yuan et al., 2019;Zhang et al., 2018) have been applied to solve numerous large deformation problems in geotechnics, such as freesurface flows (Larese et al., 2008), fluid-structure interactions (Franci and Zhang, 2018;Idelsohn et al., 2008), granular matters (Dávalos et al., 2015;Zhang et al., 2013), penetration problems (Monforte et al., 2017), subaerial and submarine landslides (Cremonesi et al., 2011;Salazar et al., 2015;Wang et al., 2021a;Zhang et al., 2020Zhang et al., , 2019)), etc.For more details, the reader is referred to the recent review (Cremonesi et al., 2020).In this study, we adopt a PFEM version developed in (Wang et al., 2021b;Zhang et al., 2013).The procedures of PFEM analysis in each time increment include: (i) updating the mesh nodes according to the incremental displacement obtained in the last time step; (ii) identifying the boundary of the computational domain based on the updated mesh nodes; (iii) generating a new mesh based on the identified boundary; (iv) mapping the history field variables to the new mesh; (v) solving the governing equations based on the new mesh.The adopted PFEM technique is based on the infinitesimal strain assumption (see Eq. ( 2)) with an updated geometry, which is an approximation for large deformation analysis.Some errors may be introduced for special scenarios like rotations, while such a treatment has proven to be efficient in handling general large deformation problems.It was initially used for frame L. Wang and Q. Lei structure analysis with large displacements in the so-called sequential limit analysis (Yang, 1993) and further demonstrated in the plane-strain analysis of the von Mises model with non-linear isotropic hardening (Leu, 2005).The sequential limit analysis has also been applied to study pipe-soil interactions during large deformations, for which simulation results agree well with experimental data (Kong et al., 2018).Moreover, the infinitesimal strain assumption with an updated geometry is the core of RITSS (Remeshing and Interpolation Technique with Small Strain), which has been extensively examined and used to analyse geotechnical problems with large deformations (Hu and Randolph, 1998;Tian et al., 2014;Wang et al., 2015).This PFEM version has been applied to investigate challenging large deformation problems such as granular material collapse (Wang et al., 2021c;Zhang et al., 2013) and landslideinduced waves (Zhang et al., 2019), for which numerical results showed good agreements with experimental observations/measurements.

Benchmarks
In this part, we verify our implemented damage model against available experimental datasets for mode-I and mixed mode failures in brittle materials.In this work, numerical simulations are all performed on a Dell laptop (Precision 3560, i7-1185G7 CPU, 32 GB RAM).The computational time (including mesh generation, computation, and data output) of each numerical example is also presented.

Model-I failure of a three-point bending beam
We apply the developed numerical model to simulate the three-point bending beam experiment reported in (Körmeling and Reinhardt, 1983).The geometry and boundary conditions are shown in Fig. 3a.The 2D domain of a plane stress condition is discretised by a grid of 1074 elements (see Fig. 3b).The material properties include (Jirásek and Bauer, 2012): the Young's modulus E = 20 GPa, the Poisson's ratio ν = 0.2, the tensile strength f t0 = 2.4 MPa, and the mode-I energy release rate G I = 90 J/m 2 .An incremental displacement Δu = 5 × 10 -4 mm is applied at the centre of the top surface and the quasi-static simulation runs for 1000 steps.The crack band size is assumed to be 1.316 and Bauer, 2012), where A e is the area of the local element.Fig. 3c shows that the numerical simulation well captures the distribution and evolution of damage driven by mode-I failure, consistent with the previous simulation results in (Jirásek and Bauer, 2012).As we can see in Fig. 3d, the displacement-loading response obtained from our numerical simulation also agrees well with the experimental dataset reported in (Körmeling and Reinhardt, 1983).The average computational time for one time step simulation is 0.612 s.

Mixed-mode failure of an L-shaped panel
We examine our model for simulating the failure of an L-shaped panel under quasi-static plane stress condition.The geometry and boundary conditions are shown in Fig. 4a.The material properties of the original experiment (Winkler, 2001) are: the Young's modulus E = 25.85GPa, the Poisson's ratio ν = 0.18, the tensile strength f t0 = 2.7 MPa, the compressive strength f c0 = 27 MPa, and the mode-I energy release rate G I = 95 J/m 2 .Here, we follow the work in (Nguyen et al., 2018) and use the calibrated values for some parameters, i.e. the calibrated Young's modulus E = 21 GPa and mode-I energy release rate G I = 115 J/m 2 , while the mode-II energy release rate G II is assumed to be the same as G I .The loading is applied at the lower surface of the L-shaped panel via an incremental displacement Δu = 8 × 10 -4 mm (see Fig. 4a), and 1000 steps of simulation is performed based on a computational mesh of 2897 elements (see Fig. 4b).The crack band size assumes to be 1.316 As we can see in Fig. 4c, the damage distribution is similar to the experimentally observed 'crack' pattern as illustrated in Fig. 4a (shaded area) and also consistent with the simulation results of (Nguyen et al., 2018).Furthermore, the displacement-loading response captured by our model agrees well with the experimental data (Winkler, 2001).The average computational time for one time step simulation is 1.64 s.

Model Setup
We further apply our model to simulate rock slope failure with the model setup shown in Fig. 5, where a roller boundary condition is applied at the left and bottom of the domain.The horizontal stress σ h at the right side of the domain is defined as σ h = α h ρgh, where α h is the lateral stress coefficient (ratio of the horizontal stress to the vertical one), ρ is the density, g is the gravitational acceleration, and h is the depth calculated from the ground surface.To capture the progressive evolution of the rock slope from a stable to an unstable state, we gradually increase α h by a prescribed formula: α h = min(2 + t/t 0 , 10) where the time t has a unit of days, and t 0 is a reference time set as 15 d, which may resemble the effect of an increased tectonic forcing or a possible surcharge load in the region right to the study domain.The hypothetical slope contains a fault zone (with a thickness of ~ 4 m) bounded by a polygon P 1 P 2 P 3 P 4 .The fault zone is simulated as an elasto-plastic layer governed by the non-associated flow rule (Eq.( 42)) with the friction angle φ = 30 • and the dilation angle ψ = 0 • .We model the damage evolution in the region P 1 P 5 P 6 P 7 , aiming to capture the interplay between the rock damage propagation and the fault zone reactivation within this failure-prone region (i.e.below the slope surface and above the basal fault zone).Other parts are assumed to be elastic for simplicity.
Material properties are assumed as follows: the density ρ = 2500 kg/m 3 , the Young's modulus E = 20 GPa, the Poisson's ratio ν = 0.2, the tensile strength f t0 = 2.4 MPa, and the compressive strength f c0 = 42.3MPa.We use the elasto-brittle constitutive law (Eq.( 14)) for modelling damage evolution in the rock slope with the residual strength ratio set as η = 0.01.The model is discretised using 22,281 uniform triangular elements.An adaptive time stepping scheme is used with the incremental displacement constrained to be smaller than the element size and the maximum time step set as 1 d.The average computational time for one time step simulation is 109.4s.
We consider two scenarios where the fault zone dips or anti-dips with respect to the slope, respectively, with the purpose of analysing the effect of fault geometry on the slope deformation and failure.To capture the behaviour of geomaterials changing from a solid-like to a fluid-like state during the failure and post-failure stages, fully damaged geomaterials (with the damage reaching its maximum value) are described as non-associated frictional material with a residual friction angle φ r = 5 • and a dilation angle ψ = 0 • .Such a frictional rheological model is commonly used for analysing the post-failure processes of flow-like  L. Wang and Q. Lei landslides (Pastor et al., 2009;Savage and Hutter, 1989) and low frictional resistance may explain the high mobility of real landslides (Iverson, 1997;Lucas et al., 2014).

Pre-failure stage
We characterise the pre-failure stage based on the distributions of damage (in the brittle rock) and equivalent plastic strain (in the fault zone), as shown in Fig. 6 and Fig. 7 for the dip and anti-dip cases, respectively.The elapsed time during the pre-failure stage is denoted as t pre .
As we can see in Fig. 6a, the equivalent plastic strain first accumulates along and within the fault zone at t pre = 1 d.Then, the damage initiates from the toe of the slope (see t pre = 4 d), while the damaged rock and reactivated fault are still not connected yet.Later, the damage further propagates and intersects with the fault zone at t pre = 17.7 d, such that the sliding body becomes mobilised for downslope movement with the mean velocity v exceeding 0.05 m/s at t pre = 59.2 d.It can be seen that at this moment, the damaged region in rock is well connected to the fault zone that has been reactivated from the toe up to the crest.The corresponding velocity distributions at these four instants, i.e. t pre = 1, 4, 17.7, and 59.2 d are shown in Fig. 6b, where the displayed velocity is normalised to the mean velocity magnitude v.At t pre = 1 d, a uniform distribution of velocity in the slope is observed with a mean value of 3.5 × 10 -7 m/s.Then, at t pre = 4 d, the velocity field becomes slightly heterogeneous with the mean velocity decreased to 2.5 × 10 -8 m/s.This may be attributed to the resistance of the slope toe region against the upslope sliding body initiated along the fault.With the further propagation of damage in rock and accumulation of plastic strain in the fault zone, the mean velocity increases to 4 × 10 -7 m/s at t pre = 17.7 d, for which the upper body moves faster than the lower part.This is also reflected by the accumulation of plastic strain within the upper part of the fault zone in Fig. 6a.Further we can see a rapid movement (0.05 m/ s) of the sliding body at t pre = 59.2 d.The slope toe region starts to move faster than the upslope region, indicating a transition of the slope kinematic from along-fault sliding to lateral movement along the basal damaged zone.
For the anti-dip case, we find that the damage propagates faster than the accumulation of plastic strain, as shown at t pre = 19 d in Fig. 7a.Then the damaged rock connects with the fault zone at t pre = 25 d (Fig. 7a), leading to the mobilisation of the slope characterised by an increased mean velocity from 1.7 × 10 -8 m/s at t pre = 19 d to 5 × 10 -7 m/s at t pre = 25 d (Fig. 7b).Even though the damage further propagates from the toe up to the crest and the plastic strain continues to accumulate within the fault zone (see t pre = 100 d and t pre = 236 d in Fig. 7a), the trapezoidshaped mass is still able to remain stable.As shown in Fig. 7b, at t pre = 236 d, the velocity distributions are quite uniform with the mean velocity finally decreased to a very small value of 4.6 × 10 -14 m/s, indicating the stabilisation of the slope.

Post-failure stage
We present the post-failure evolution of the rockslide (note: only the case with a dipping fault is shown since the anti-dipping one does not fail) after the failure at t pre = 59.2 d (Fig. 6a), for which a velocity threshold (0.05 m/s) for the mean velocity is used to distinguish the preand post-failure stages.For the post-failure stage, the time is denoted as t post such that the starting time of the post-failure stage at t post = 0 s corresponds to the ending time of the pre-failure stage at t pre = 59.2 d.To elucidate the transition of the slide from a solid-like to a fluid-like state: we first present the damage and velocity distributions within the sliding body at selected instants in Fig. 8; we then display the subsequent runout processes of the fully mobilised mass in Fig. 9.
As one can see in Fig. 8, the damage further propagates during the early stage after failure (see t post = 0.15 s), which is characterised by the large velocities distributed along the slip surface (see Fig. 8b).In a few seconds, we observe the generation of damage within the sliding body (see t post = 7 s in Fig. 8a), and the velocity increases from 0.12 m/s (at  L. Wang and Q. Lei in Fig. 8a and b).Further, the mass becomes totally damaged at t post = 9.56 s.
We present the simulation results of the runout process of the failed mass in Fig. 9, the normalised velocity distribution of the flowlike materials at four instants are shown.The flow has a front/centreconcentrated shape with a thin tail and the mean velocity of the mass increases to 8.32 m/s at t post = 10.5 s (compared to the mean velocity of 7.01 m/s at t post = 9.56 s in Fig. 8b).Then, the mass accelerates to 8.68 m/s (at t post = 11.5 s) due to the gravity-driven movement.Afterwards, it spreads in an elongated shape and deaccelerates as a result of the energy dissipation from plasticity, as shown at t post = 13.5 s.The mass flow further impacts the left boundary (where a roller condition is imposed, representing a possible cliff at the left of the valley) with the mean velocity magnitude being 5.86 m/s at t post = 14.86 s.The run-out distance is not analysed here due to this boundary effect.

Discussion
In this section, we discuss our modelling results from geological and geotechnical perspectives.We first present the capability of the model in capturing realistic displacement time history and discuss its implications.We then demonstrate that our model can provide useful insights into the rockslide transition process (from a solid-like state to a fluid-like state), which may not be captured by the conventional post-failure analyses.Finally, we discuss the capability of our PFEM model for landslide modelling.

Pre-failure stage
To further quantitatively analyse the pre-failure evolution of the slope, in Fig. 10, we plot the displacement time series of the monitoring point P 5 at the top surface (see Fig. 5 for the location of this point).For the case with a dipping fault, we observe a progressive slope deformation with a displacement rate of ~ 10 mm/d starting from t = 20 d and then a rapid acceleration of slope movement exceeding 1 m/d at t = 59 d when it is close to failure.For the case with an anti-dipping fault, the slope episodically accelerates but globally self-stabilises.Such progressive slope deformation/failure phenomena captured in our model resemble many of the features of natural rockslides as observed in the field (Crosta et al., 2017;Crosta and Agliardi, 2003;Loew et al., 2017;Petley et al., 2005).Indeed, the use of continuum methods with damage mechanics models have long been regarded as an efficient way to simulate the progressive development of rock slopes (Brückl and Parotidis, 2005;Lacroix and Amitrano, 2013;Riva et al., 2018).In this work, we generically study the rockslide failure triggered by an increased lateral stress, which may resemble the triggering situations due to tectonic loading, building construction, landfill or other engineering activities.Note that the prescribed stress change rate may affect the slope response (Spreafico et al., 2021), which will be explored in our future work.To better mimic the long-term strength degradation of large rock slopes in mountainous regions, time-dependent damage models (Lacroix and Amitrano, 2013;Riva et al., 2018;Spreafico et al., 2021) can be integrated into our computational framework in the future for simulating the creep behaviour of rock slopes.

Post-failure stage
For the post-failure stage, a large mass of rock may experience significant strength degradation and the rock debris can move very long distances at extremely high speeds.The long runout mechanisms may be attributed to friction degradation, liquefaction, and entrainment processes (Davies and McSaveney, 1999;Hungr and Evans, 2004).With the incorporation of different rheological models that can describe the strength degradation, various numerical models have been developed to predict the mobility of rockslides (Denlinger, 2004;Hungr and Evans, 2004;Lucas et al., 2014).For such post-failure analyses, the slip surfaces are usually pre-assumed and the frictional materials are released from a static state (Denlinger, 2004;Mergili et al., 2017;Pastor et al., 2009), which means that the failure-to-runout transition process is not physically captured.To demonstrate the differences in modelling concepts/ procedures between our simulation of the full rockslide process (including pre-and post-failure stages) and the conventional post-failure analysis, we conduct a comparative analysis for the scenario of rock slope with a dipping fault.To describe the friction degradation behaviour during motion, we adopt a strain-weakening model that updates the friction angle with the deviatoric plastic strain as (Potts et al., 1990;Troncone, 2005): where φ m , φ p and φ r are the mobilised, peak, and residual friction angles, respectively.k shear is the equivalent deviatoric plastic strain that is evaluated by the time integration of its rate k shear = ∫ kshear dt.The rate  .The rate of deviatoric plastic strain tensor ėp ij is defined as ėp ij = εp ij − 1 εp kk δ ij , in which δ ij is Kronecker's delta, and εp ij is the rate of the plastic strain tensor.k p shear and k r shear are two thresholds.In the simulation, strength parameters are set as the same in the full process simulation: φ p and φ r are 30 • and 5 • , respectively, while the dilation angle ψ = 0 • ; we choose k p shear = 0.1 and k r shear = 1.The developed PFEM model with elasto-plastic relationships has been validated against other numerical simulations and observation data, for run-out analysis of landslides (Wang et al., 2021a(Wang et al., , 2021b(Wang et al., , 2021c)).
We display the equivalent deviatoric plastic strain and the mean velocity distribution of the masses at four instants in Fig. 11.As can be seen in Fig. 11a, k shear immediately accumulates along the toe and the fault (see t post = 0.32 s) and rapidly connects with each other, forming a shear surface at t post = 0.72 s.Macroscopically, this process is reflected as a sliding body that rapidly moves on a stable zone in Fig. 11b, which agrees well with the observations of granular collapses in laboratory experiments (Balmforth and Kerswell, 2005;Lube et al., 2004).In a very short time (i.e.t post = 0.32 s), the mean velocity increases to 1.6 m/s, and then goes to 3.29 m/s at t post = 0.72 s.Along this shear surface (see t post = 0.72 s in Fig. 11a), the strength of material decreases as a result of strain-weakening effects and the lower and upper parts evolve into the residual state (characterised as dark red scatters).Further, several internal shear zones are observed and the rapid sliding movement causes a distinct change in the slope geometry (t post = 1.63 s in Fig. 11a).The sliding body moves downward along the slip surface with an extremely rapid mean velocity of 6.55 m/s (see t post = 1.63 s in Fig. 11b).During the rapid motion of the rockslide, the equivalent deviatoric plastic strain quickly spreads within the sliding body (see t post = 2.56 s in Fig. 11a) and the mean velocity increases to 8.58 m/s (see t post = 2.56 s in Fig. 11b).
Compared to the transition process captured in the full process model (the material degradation within a sliding body along a well-developed slip surface) shown in Fig. 8, in the post-failure only simulation, the rock masses are rapidly accelerated (mean velocity increases from 0 to 1.60 m/s in 0.32 s) before the accumulation of equivalent deviatoric plastic strain (see t post = 0.32 s in Fig. 11), due to the gravitational sliding.The rapid movements of materials significantly change the slope geometry in ~ 2 s (see t post = 1.63 s and t post = 2.56 s in Fig. 11), generating shear bands within the sliding body.As shown in Fig. 11, the movement of the rockslide is rotational, and the strength degradation does not affect the failure pattern.While for the full process simulation shown in Fig. 8, we observe the entire sliding block evolves into a flow resulting from the damage propagation, and this transition occurs over a certain time window (~8 s).Such a transition from a solid state to a granular state has been highlighted in the progressive evolutionary changes and different dynamic regimes of fault zones, and damage rheology provides a bridge between these two states (Ben-Zion, 2008).In this work, our model mainly phenomenologically represents such a transition, since we do not calculate the residual friction angle for damaged elements according to the local stress state.More precise treatments (e.g.stressdependent residual friction) will be adopted in our future works.
In Fig. 12, we present the simulation results of the failed rock mass (at the residual state), showing the normalised velocity distribution of the sliding body at four instants.The flow has a centre-concentrated shape with the mean velocity of the mass increases to 9.01 m/s at t post = 3.29 s (compared to the mean velocity of 8.58 m/s at t post = 2.56 s in Fig. 11b).Then the mass spreads and deaccelerates to 8.21 m/s at t post = 4.89 s, at which a thin front and tail are observed.The sliding body rapidly moves with the spread of the concentrated mass, shown as a thin layer of materials (see t post = 6.41 s).Afterwards, the flow reaches the left boundary of the model with the mean velocity being 6.47 m/s, as can be seen at t post = 7.19 s.
From the translational movement of the sliding body shown in Fig. 9 and Fig. 12, we observe that the rockslide profiles (t post = 10.5 s in Fig. 9 and t post = 3.29 s in Fig. 12) are different at the early stage after gravitational sliding.Once the rock masses are fully mobilised (residual frictional angle is applied), the simulated rockslides both rapidly spreads and reach the left boundary in ~ 4 s.In contrast to the performed post-failure only analysis showing the mass sliding along a circular slip surface, the full process simulation can capture the transition of the rockslide from a solid-like to a fluid-like state, resulting from the damage propagation.In rockslide problems, this may related to the mechanical fluidisation and/or dynamic fragmentation of debris that are commonly used to explain the high mobility of rock avalanches (Davies and McSaveney, 1999;Davies, 1982;Taboada and Estrada, 2009).After the transition (the rock masses are at residual state), the flow-like rockslide evolves into a translational mode, while two different shapes of the flow are observed at the beginning of translational movement, according to the full process simulation and the post-failure simulation.Such a difference in the longitudinal variation of material distribution can be further investigated to better understand the material flow surges by considering solid-fluid interactions (Hungr, 2000;Iverson, 1997).

PFEM model for landslide modelling
In the past, various methods and tools have been developed and applied to analyse the complex mass movements at several spatial and temporal scales (van Westen et al., 2006).To analyse and forecast landslide hazards, mathematical models that can quantify the relationships between the inputs (e.g.triggers) and outputs (e.g.slope movements) can be very useful (van Asch et al., 2007).Even though it is almost impossible to completely represent the characteristics of slope movements subject to multiple factors (e.g.geometrical, geotechnical, structural, and hydrological properties), the mathematical models solved by numerical techniques have been found useful for landslide hazard assessment (Eberhardt, 2008;Guzzetti et al., 2002;Mergili et al., 2017;Pastor et al., 2009;Thiebes et al., 2014;van Asch et al., 2007), where the pre-and post-failure processes are, however, separately modelled in most previous studies.Consequently, this lack of a unified mathematical model leads to difficulties and uncertainties in landslide hazard assessments (Cascini et al., 2013), because the destabilised rock slope may not evolve into a rapid avalanche after the generation of a slip surface (Cascini et al., 2013;Cuomo et al., 2013;Soga et al., 2016;Wang et al., 2021b).This indicates that the conventional runout analysis (assuming the identified landslide volume is mobilised along the preassumed slip surface) may not be applicable to some real landslide scenarios.For example, if we perform a post-failure simulation (e.g. using PFEM or depth-averaged models) on the slip surface identified by FEM analyses (e.g. the slope with an anti-dipping fault shown in Fig. 7), we will get similar numerical results for the post-failure stage as shown in Fig. 11 and Fig. 12, which however do not reflect the self-stabilised post-failure behaviour of this slope.Indeed, the modern numerical methods can serve as a basis for developing more advanced models for landslide hazard assessment, since it provides a computational framework that unifies the simulation of the pre-and post-failure stages (Cuomo et al., 2021;Dey et al., 2016;Islam et al., 2018;Soga et al., 2016;Wang et al., 2021a).If the variation of material strength can be constrained by the real-time monitoring data (Crosta et al., 2017), such a unified computational framework may be integrated into geoscientific models that are used for describing or predicting the entire landslide process in practice.Nevertheless, numerous efforts have been paid to geomechanical analyses of landslides with most in the run-out analysis of landslides (Cremonesi et al., 2011;Franci et al., 2020;Huang et al., 2012;Llano-Serna et al., 2016;Pastor et al., 2009;Yerro et al., 2016;Zhang et al., 2020) and few in the entire process simulation (including both pre-and post-failure stages) (Bandara et al., 2016;Cuomo et al., 2021;Islam et al., 2018;Wang et al., 2021a).To overcome this issue, the PFEM might be a good choice, since the standard FEM code can be readily extended to include the particle finite element technique (Cremonesi et al., 2010;Zhang et al., 2013).
The application of PFEM to landslide problems starts from its capability of simulating free-surface flow problems, where the post-failure behaviour of landslides can be described using the Newtonian or non-Newtonian rheological laws (Cremonesi et al., 2011;Franci et al., 2020;Idelsohn et al., 2008Idelsohn et al., , 2004;;Larese et al., 2008).By adopting the efficient remeshing procedure (Cremonesi et al., 2010) and a practical remapping method (Hu and Randolph, 1998), a geomechanics-based PFEM has been developed for large deformation problems (Zhang et al., 2013), in which the widely used elasto-plastic models for solids are used (e.g.Mohr-Coulomb criterion).Additionally, a smoothing strain technique is implemented into such a geomechanics-based PFEM to avoid remapping and reduce remeshing operations (Yuan et al., 2019;Zhang et al., 2018Zhang et al., , 2022)).For the pre-failure simulation, the geomechanics-based PFEM naturally inherits the flexibility and robustness of FEM for the stress-strain analysis of slope problems, since the numerical solutions are obtained from the Lagrangian FEM analysis.As for the post-failure simulation, the geomechanics-based PFEM has been applied to the runout analysis of real landslide cases (Jin et al., 2020;Wang et al., 2021b;Zhang et al., 2020Zhang et al., , 2015)), where numerical results agree well with other independent numerical solutions and observation data.Further, the shallow and deep failure mechanisms and subsequent runout processes in layered clayey slopes under seismic loadings have been studied (Wang et al., 2021a).Among these simulations, the prefailure long-term destabilisation has been rarely studied (partly due to the small time step required in some explicit time integration schemes), which therefore has not fully demonstrated the capability of PFEM for landslide modelling.Moreover, some modules are still under development, and the developed hydro-mechanical modules (Monforte et al., 2019;Wang et al., 2021c) have not yet been applied to real landslide cases.
In this work, we extend the previously established PFEM model to simulate the pre-and post-failure behaviour of rock slopes.Here, the computational framework is implemented in a fully implicit time integration scheme based on a generalised Hellinger-Reissner variational principle, which permits a large time step with guaranteed numerical L. Wang and Q. Lei stability.Consequently, the model can capture the landslide evolution during the pre-and post-failure stages operating across different timescales.The incorporation of damage mechanics models enhances the capability of the PFEM model in realistically simulating progressive failure processes in rock slopes and Parotidis, 2005;Lacroix and Amitrano, 2013;Riva et al., 2018).Moreover, such a unified computational framework opens the possibility of capturing the complex transition processes from a destabilised rock slope to a catastrophic runout.In the future, this computational framework will be applied to systematically study the progressive rock slope failure phenomena in natural systems based on available site characterisation and in-situ monitoring data.

Conclusions
In this work, we established a novel computational framework for simulating the pre-and post-failure behaviour of rock slopes by implementing a continuum damage mechanics model (Jirásek and Bauer, 2012) into the particle finite element method (Wang et al., 2021b).The computational framework is compatible with the accurate and efficient implicit SOCP-based finite element methods (Krabbenhoft et al., 2012(Krabbenhoft et al., , 2007;;Wang et al., 2021b;Zhang et al., 2013), for which the mathematical formulation is rigorously derived based on a generalised Hellinger-Reissner variational principle.The developed model is validated against experimental data of brittle material failure, and further demonstrated for rock slope problems based on a hypothetical model considering rock damage and fault reactivation.We show that the model can simulate the entire process of rock slope failures, with realistic prefailure displacement time history and post-failure runout evolution captured.The present work opens a door for the development of a unified computational framework for analysing the initiation, propagation, and deposition of rockslides.L. Wang and Q. Lei

Fig. 1 .
Fig. 1.Schematics illustrating the constitutive laws of damage in brittle materials: (a) stress-strain and (b) damage-strain relationships in the damage model with/ without a post-peak softening.

Fig. 2 .
Fig. 2. The mixed triangular element adopted in the numerical model.

Fig. 3 .
Fig. 3. Mode-I failure of a three-point bending beam: (a) experimental set-up; (b) computational mesh; (c) computed damage pattern; (d) comparison of the force-displacement response between the experiment data (Körmeling and Reinhardt, 1983) and our numerical solution.

Fig. 4 .
Fig. 4. Mixed-mode failure of an L-shaped panel: (a) experimental set-up; (b) computational mesh; (c) computed damage pattern; (d) comparison of the force--displacement response between the experiment data(Winkler, 2001) and our numerical solution.

Fig. 5 .
Fig. 5. Geometry and boundary conditions of the hypothetical rock slope.

Fig. 6 .
Fig. 6.Pre-failure evolution of the slope with a dipping fault: (a) damage pattern and equivalent plastic strain distribution at four instants; (b) normalised velocity distribution of the sliding body at four instants.

Fig. 8 .
Fig. 8. Post-failure evolution of the slope with a dipping fault: (a) damage distribution at four instants (only the damaged areas with ω > 0.1 is plotted); (b) normalised velocity distribution of the sliding body at four instants.

Fig. 9 .
Fig. 9. Normalised velocity distribution of the flowing mass during the post-failure stage at four instants.

Fig. 11 .
Fig. 11.Post-failure only simulation of the rockslide with a dipping fault at four instants: (a) equivalent deviatoric plastic strain (only areas with k shear > 0.1 is plotted); (b) normalised velocity distribution of the sliding body.

Fig. 12 .
Fig. 12. Normalised velocity distribution of the flowing mass at four instants in the post-failure simulation.

Fig. A2 .
Fig. A2.The stress distribution around a single fracture in a body under normal stresses.