Mathematical Optimization Problems for Particle Finite Element Analysis Applied to 2D Landslide Modeling

Notwithstanding its complexity in terms of numerical implementation and limitations in coping with problems involving extreme deformation, the finite element method (FEM) offers the advantage of solving complicated mathematical problems with diverse boundary conditions. Recently, a version of the particle finite element method (PFEM) was proposed for analyzing large-deformation problems. In this version of the PFEM, the finite element formulation, which was recast as a standard optimization problem and resolved efficiently using advanced optimization engines, was adopted for incremental analysis whilst the idea of particle approaches was employed to tackle mesh issues resulting from the large deformations. In this paper, the numerical implementation of this version of PFEM is detailed, revealing some key numerical aspects that are distinct from the conventional FEM, such as the solution strategy, imposition of displacement boundary conditions, and treatment of contacts. Additionally, the correctness and robustness of this version of PFEM in conducting failure and post-failure analyses of landslides are demonstrated via a stability analysis of a typical slope and a case study on the 2008 Tangjiashan landslide, China. Comparative studies between the results of the PFEM simulations and available data are performed qualitatively as well as quantitatively.


Introduction
As a widely found geophysical phenomenon, the term "landslide" refers to various mass movements on slopes, including rockfalls, topples, debris flows, etc. (Varnes 1978). The entire process of a landslide event contains three stages, namely initiation, motion, and deposition. It is common to study the initiation stage by means of various forms of stability analysis, while the motion and deposition stages are considered by means of approaches capable of addressing extreme material deformation.
Based on either solid or fluid mechanics, numerous efforts have been devoted to the investigation of complicated landslide behavior at different stages. Typical numerical approaches for analyzing the initiation of landslides include, but are not limited to, the limit equilibrium method (Fredlund and Krahn 1977), limit analysis method (Chen 1975), and finite element method (Dawson et al. 1999). Among these, use of the FEM, with its flexibility in handling complex geometries, sophisticated constitutive models, various loading conditions, and provision of the time evolution information of the slope body, is prevalent for slope stability analysis. For the motion and deposition stages, the sliding mass is commonly represented by Lagrangian block models (Hungr 1995;Tinti et al. 1997) or flow models (Iverson 1997;Savage and Hutter 1989) based on depth-averaged variables, which are solved numerically using finitedifference schemes (Tai et al. 2002), or by smoothed particle hydrodynamics (SPH) approaches (Pastor et al. 2009). Other approaches capable of modeling the post-failure processes include the discrete element method (Staron 2008), material point method (Llano-Serna et al. 2016), coupled Eulerian-Lagrangian method (Crosta et al. 2003), etc. Notably, although the traditional Lagrangian FEM performs well for slope stability analysis, it cannot capture the motion and deposition stages, since severe mesh distortion is encountered when the sliding mass experiences large changes in geometry. To tackle this issue, the particle finite element method (PFEM), which combines standard finite element analysis and a particle-based technique, was proposed by Idelsohn et al. (2004). Although originally developed to consider fluid-solid interaction problems with free surfaces in fluid mechanics, the PFEM has demonstrated its capability to model the motion and deposition stages of landslides due to the works by Cremonesi et al. (2011), Oñate et al. (2011), andZhang et al. (2015). Moreover, since the FEM is utilized in the PFEM for each incremental analysis, slope stability analysis can also be performed in the framework of the PFEM, implying an ability to model the entire landslide process, from its initiation through the sliding process to the deposition, in a single simulation.
According to the solution scheme, the versions of the PFEM developed for landslide modeling can be classified into two categories. The PFEM developed by Oñate et al. (2004) and by Cremonesi et al. (2010) belongs to the first category that solves the resulting nonlinear finite element equations using a nested scheme based on Newton's scheme or a variant thereof. In this solution scheme, iterations are carried out between the level of global structures (where the unbalanced force is minimized) and the level of material points such as the stress integration points (where the stress-strain relationship is fulfilled). This is also the most widely used scheme for the implementation of the FEM in both in-house and commercial numerical packages. In the second category (optimization-based PFEM), the FEM formulation is resolved in mathemat-ical programming. Specifically, the boundary-value problem for landslide analyses is converted into an equivalent optimization problem, for instance a second-order cone programming (SOCP) problem, which is then solved by using the interior point method (IPM). This version of the PFEM inherits some unique merits of the SOCP-FEM, such as natural treatment of the singularities in the Mohr-Coulomb and Drucker-Prager yield criteria and robust extension to handling of multisurface plasticity. One more noticeable advantage of this solution scheme, which is vital to landslide modeling, is its admirable convergence properties. Due to its variational nature, the existence, uniqueness, sensitivity, and stability of the solution can be analyzed mathematically. Indeed, efforts have been devoted to the analysis of the convergence properties of the utilized IPM for solving nonlinear optimization problems; For example, the global and local convergence properties of the IPM for nonlinear programming were demonstrated by Tits et al. (2003), while its stability and convergence rates were analyzed mathematically by Alizadeh et al. (1998). Recently, the possibility of using the optimization-based PFEM for simulating challenging landslide problems has been explored through modeling a flow-like landslide (Zhang et al. 2015), progressive failure of landslides in sensitive clays (Zhang et al. 2017), as well as submarine landslides and their consequences such as induced tsunamis and the impact on ocean pipelines .
Notwithstanding its powerful capability, few contributions have been published regarding the implementation of the optimization-based PFEM. Indeed, the theory of the optimization-based PFEM has been well documented and verified by Zhang et al. (2013). However, its numerical implementation, which is widely different from the version of PFEM based on Newton's iteration schemes, has not been introduced in detail, which therefore handicaps its further applications. We remark that the solution scheme for the finite element formulation in the optimization-based PFEM differs considerably from that in the version developed by Oñate et al. (2004) and Cremonesi et al. (2010). Additionally, we note that, although the capability of the optimizationbased PFEM for landslide modeling has already been explored, there are only a few efforts devoted to studies comparing the optimization-based PFEM versus the common approaches for landslide modeling in geosciences (e.g., versus depth-averaged approaches). This work aims to contribute to fill these gaps with a twofold objective: (1) to detail the numerical implementation of the optimization-based PFEM to assist researchers from geoscience in developing their own version, and (2) to provide a quantitative comparison between the optimization-based PFEM and the techniques commonly adopted in geosciences to treat landslides, including both stability analysis and post-failure analysis.

Governing Equations
In this section, the governing equations for quasistatic and dynamic analyses of rateindependent elastoplastic problems under plane-strain conditions are summarized. The equations consist of the equilibrium equation, strain-displacement relation, constitutive equation, and boundary conditions. Additionally, the contact between the deformable solid (e.g., sliding mass) and rigid surface (e.g., basal surface) is also con-sidered in the formulation. For the dynamic analysis, the θ time-integration method (Bathe and Wilson 1973) is used to discretize the governing equations in the time domain.

Static Analysis
For a two-dimensional (2D) domain V , delimited by a boundary S, the set of governing equations for a static analysis is as follows: (a) The equilibrium equation where the strain is ε ε x x , ε yy , 2ε xy T and the displacement is u u x , u y T .
(c) The boundary conditions on S where N n x , 0, n y ; 0, n y , n x , t is the traction force, and u p is the prescribed displacement. (d) The constitutive relationship in which 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 G is the plastic potential. As shown in Eq. (2.4b), the additive decomposition of the strain is used, whose incremental form is For an associated flow rule, we have G F. When the material undergoes purely elastic deformation, the plastic multiplier increment λ 0 and F (σ ) < 0, whereas when the material yields, we have λ > 0 and F (σ ) 0. This constraint is the so-called complementary condition that λF (σ ) 0, λ ≥ 0, (2.6) which implies that the incremental form of the constitutive relationship is as follows

Dynamic Analysis
For dynamic analyses, the inertial force should be included in the equilibrium Eq. (2.1), that is (2.8b) By means of the θ -method (Bathe and Wilson 1973), the stresses and velocities can be rewritten as σ (2.9b) Hereafter, the subscripts n and n + 1 denote the known and unknown states, and t is the time increment. By introducing a new intermediate variable, i.e., the inertial force γ , whose definition is shown in Eq. (2.12), and substituting Eq. (2.9b) into Eq. (2.9a), the latter can be rearranged as (2.10) and the according traction boundary condition Eq. (2.3a) becomes (2.12) Equation (2.9b) can also be rearranged as which is used to update the velocity v n+1 after the displacement increment u is obtained. It should be mentioned that the time integration scheme is unconditionally stable when θ 1 ≥ 1 2 and θ 2 ≥ 1 2 .

Contact Analysis
The rigid no-penetration contact scheme is adopted to consider the interaction between the sliding mass and the basal surface, implying that, for any points of the sliding mass, the incremental displacement of the point cannot exceed the gap between the point and the boundary, as illustrated in Fig. 1. Therefore, the corresponding contact conditions are as follows where n is the outward unit vector of the rigid surface, g 0 is the gap between the boundary point of the deformable material and the rigid surface at t t n , p is the normal force, q is the tangential force, and μ is the friction coefficient.

Equivalent Min-Max Program
The aforementioned governing equations can be reformulated as equivalent min-max optimization problems according to the works of Simo et al. (1989), Krabbenhoft et al. (2007), and Zhang et al. (2013), where detailed derivations are documented. In this section, a brief summary is provided for those equivalent min-max optimization problems.
The Hellinger-Reissner variational principle (Reissner 1950) is adopted to reformulate the quasistatic elastoplastic problem. According to Krabbenhoft et al. (2007), the equivalent optimization problem is in the form of As shown, both the displacement and stress are independent master fields for the above optimization problem. The optimal solution of problem (3.1) is a saddle point that renders neither the maximum nor the minimum of the functional. The incremental form of (3.1) is which is the one used for the incremental finite element analysis of quasistatic elastoplastic problems. Efforts have also been devoted to reformulate the governing equations for dynamic elastoplastic problems. It was demonstrated by Zhang et al. (2013) that the min-max optimization problem equivalent to the governing equations discretized by the θmethod for dynamic analysis is where the inertial force γ is included as an independent master field in the maximum part of the optimization problem. The contact constraints, i.e., (2.14), will be taken into account later.

Mixed Triangular Element
Using standard finite element notations, the master fields of the min-max program are approximated as follows whereσ ,û andΥ are the stress, displacement, and inertial force at mesh nodes and N σ , N u , and N γ are the matrices consisting of the interpolation shape functions. Specifically, the mixed isoparametric triangular element shown in Fig. 2 is selected in our simulation, where a quadratic shape function is used for the fields of the displacement and inertial force (N u N γ ) whereas a linear shape function is used for the stress field.

Finite Element Discretization
Using the finite element discretization, the discrete form of the min-max problem (3.3) for the dynamic analysis of elastoplastic problems is In the optimization problem (4.4), the yield function (4.4b) is imposed at the stress interpolation points, with n σ andσ j n+1 being respectively the total number of stress interpolation points and the stress states at the jth stress integration point at t n+1 . The minimization part of the min-max problem (4.4a) can be resolved analytically, meaning that the min-max problem (4.4) is equivalent to the following maximization problem which is the one for dynamic elastoplastic analysis. The optimization problem (4.6) can be degraded to the one for quasistatic elastoplastic analysis by removing the terms relevant to dynamics and setting θ 1 θ 2 1. Specifically, the discretized optimization problem for quasistatic elastoplastic analysis is Zhang et al. (2013), further consideration of contact constraints, i.e. (2.14), transforms the optimization problem (4.6) into in which ρ (ρ x , ρ y ) T are the nodal contact forces in the global coordinate system relating to the normal and tangential contact forces p and q via n andn. n is the unit vector normal to the rigid boundary, andn (−n 2 , n 1 ) T . Assuming that the inclination angle of the slope is θ s (Fig. 1), the corresponding components are n 1 sin θ s and n 2 cos θ s . The logical index set of contact nodes is denoted by E c , and n c is the total number of potential contact nodes. The potential contact nodes are set as the nodes on the surface S of the computational domain V . The present maximization problem is a type of convex optimization problem that can be recast as a standard second-order cone programming (SOCP) problem and then resolved using the interior point method. For simplicity, the present scheme for the finite element solution is abbreviated as SOCP-FEM.

Numerical Implementation
Although the theories of the SOCP-FEM analysis for quasistatic and dynamic elastoplastic problems have been well demonstrated, there are only few efforts devoted to its numerical implementation. This section focuses primarily on the numerical implementation of the SOCP-FEM, with emphasis on some key aspects. In particular, we focus on (1) the implementation of the boundary conditions in the SOCP-FEM, (2) the transformation of the resulting maximization problem for quasistatic/dynamic elastoplastic analyses into a standard SOCP problem, and (3) the solution of the resulting SOCP problem using the IPM available in MOSEK (MOSEK ApS 2019).

Implementation of Boundary Conditions
The traction and displacement boundary conditions have to be imposed so that the boundary-value problem can be resolved. The traction boundary condition, for instance Eq. (2.3a), is handled by integrating tractions along the boundary surface, resulting in equivalent nodal forces [e.g. Eq. (4.5e)]. In other words, the implementation of the traction boundary condition in the SOCP-FEM is exactly the same as that in the traditional displacement-based FEM (Bathe 2006).
Nevertheless, the imposition of the displacement boundary condition (2.3b) in the SOCP-FEM differs from that in the displacement-based FEM. In the traditional displacement-based FEM, the displacement boundary condition is implemented either by the penalty method or by modifying the global stiffness matrix. In contrast, the SOCP-FEM requires the introduction of a new field variable, i.e., the nodal reaction forcer n+1 , for this purpose. Specifically, for quasistatic problems, the discretized optimization problem with displacement boundary conditions being fulfilled is in the form of where E u , consisting of 0 and 1, is the index matrix to indicate the prescribed displacementû p . The underlined terms in the objective function and the equality constraint in problem (5.1) are induced by the displacement boundary conditions. To check its validity, the Lagrangian associated with (5.1) is constructed, which is Differentiating the Lagrangian with respect to the reaction forcer n+1 leads to which apparently is the displacement boundary condition (2.3b). For dynamic problems, the imposition of displacement boundary conditions is achieved in the same manner.

Reformulation as a SOCP Problem
The aforementioned optimization problem [i.e. (4.8)] can be reformulated as a standard second-order cone programming (SOCP) problem in the form of where x (x 1 , x 2 , …, x n ) T is the vector of optimization variables, a, b, and c are the matrix and vectors of factors, and K is a tensor product of second-order cones such that K K 1 × · · · × K l . The second-order cones can be of the following types: • Quadratic cone Specifically, the SOCP problem equivalent to the maximization problem (4.8) is where x n+1 is a vector consisting of all the optimization variables [i.e. (5.9)].

Remark 1
The SOCP problem (5.6) is obtained by first converting the quadratic term 1 2 σ T C σ in the objective function of (4.8) into the minimization of a scalar variable [i.e., variable m in (5.6a)] subject to linear equalities and rotated quadratic constraints [i.e. (5.6c)]. The same operation is applied to the quadratic term 1 2 t 2γ T n+1 Dγ n+1 in the objective function, resulting in the variable s in (5.6a) and the constraints in (5.6d). Contact constraints in (4.8d-f) are reformulated as those in (5.6e), in which quadratic constraints are enforced. The yield criterion F σ j n+1 ≤ 0 is equivalent to the constraints in (5.6f).

Solution Using MOSEK
The resulting SOCP problem (5.6) can be resolved using the IPM, which is a robust solution scheme available in MOSEK (MOSEK ApS 2019). In this section, the program submitted to MOSEK for the solution is presented in detail. All information related to the program is stored in an object called "prob" in MOSEK. Specifically, the optimization variables of the SOCP problem (5.6) are as follows prob · bl c prob · buc f ; d; −C 1 2σ n ; 1; 1; 0; 0; 0; 0; 0 (5.10c)

Remark 2
The linear equality constraint (5.4b) is implemented as the lower and upper bounds, namely prob·blc ≤ prob·ax ≤ prob·buc.
The inequality constraints in (5.6) are in the form of quadratic or rotated quadratic cones. There are a total of n c + n σ + 2 conic cones according to the expressions in (5.6c-f). Specifying the type of each cone and the index of its members in x, the conic constraints can be defined straightforwardly in "prob." A detailed solution scheme for the SOCP-FEM analysis is provided in Algorithm 1 for reference.

Algorithm 1. Solution scheme
For the ith incremental analysis:

1.
Update the status of variables such as the velocity, the acceleration and the stress; A, B, C, D (1) Calculate g 0 of potential contact nodes;

Remark 3
The optimization solver returns the structure type of the solutions named "res" at the end of each analysis step, and the solutions of x are contained in res.sol.itr.xx. The incremental displacement û is the dual variable ofσ n+1 and is stored in res.sol.itr.y.

Particle Finite Element Technique
Originating in fluid mechanics, the PFEM has demonstrated its ability to tackle issues such as free surface evolution and mesh distortion. Some challenging fluid mechanics problems that have been solved successfully include, but are not limited to, modeling of free surface flows and their interaction with solid structures, wave breaking, multiphase flows, etc. (Idelsohn et al. 2003Oñate et al. 2004Oñate et al. , 2008Oñate et al. , 2011. The key idea behind the PFEM is the treatment of mesh nodes as free particles that can move and even separate from the computational domain to which they originally belonged. In a given time interval [t n , t n+1 ], the basic steps of the PFEM are as follows (see also Fig. 3): 1. Erase the mesh topology and update the position of mesh nodes based on the solved incremental displacement to obtain a cloud of particles, C n+1 (Fig. 3a, b); 2. Use the α-shape method (Edelsbrunner and Mücke 1994) to identify the new computational domain Ω n+1 based on C n+1 (Fig. 3c); 3. Remesh the domain Ω n+1 to obtain a new mesh M n+1 (Fig. 3d); 4. Map the history variables from the old mesh M n to the new mesh M n+1 using the unique element method introduced by Hu and Randolph (1998); 5. Solve the equations using Algorithm 1 based on the new mesh M n+1 .
It is worth noting that the governing equations used in this version of PFEM are under the assumption of infinitesimal strain. At the end of each incremental analysis, Fig. 3 Basic steps of particle finite element method the configuration is updated according to the solved displacement. This assumption may lead to several errors for large-deformation analysis. However, practically, the error induced due to the infinitesimal strain assumption is relatively minor, because the used incremental step is very small. The present version of the PFEM with smallstrain theory has been validated against several benchmarks such as the modeling of a water dam break, granular column collapse, underwater granular flows and the related induced waves, and non-Newtonian flows in an annular viscometer by Zhang et al. (2019), where satisfactory agreement between the PFEM simulation results and available experimental data and/or analytical results was obtained. For the PFEM derived from the concepts of large-strain plasticity and its application to granular material flows which are closely related to landslide propagations, we refer the reader to Dávalos et al. (2015).

Slope Stability Analysis by SOCP-FEM
Practically, the stability analysis of a slope is carried out using the strength reduction method (SRM) to identify the critical state of the slope by gradually reducing the strength of the soil. The critical state is indicated by the factor of safety (FOS), defined as the ratio of the actual soil shear strength to the minimum shear strength required to prevent failure (Bishop 1955); For example, when the Mohr-Coulomb model is used and the cohesion and internal friction angle are represented by c and ϕ, then according to the SRM (Dawson et al. 1999), these parameters are reduced by a reduction factor (RF), viz. The RF is the FOS of the slope when the used c and tan ϕ are at the minimum values required to prevent failure. In other words, a slope with a FOS greater than one is identified as stable; otherwise, it is unstable.
To verify the accuracy of the SOCP-FEM approach for slope stability analysis, a homogeneous soil slope that has been analyzed using the limit equilibrium method (LEM) and the finite element method by Cheng et al. (2007) is reexamined in this section. The initial geometry of the example is illustrated in Fig. 4. The density, elastic modulus, and Poisson's ratio of the soil are 2000 kg/m 3 , 14 MPa, and 0.3, respectively. Simulations with both the associated (ψ ϕ) and nonassociated (ψ 0 • ) flow rules are conducted using SOCP-FEM under quasistatic assumptions with material parameters in line with those in Cheng et al. (2007). For the case of the nonassociated flow rule, the strategy presented in Krabbenhoft et al. (2012) is utilized. A total of 6095 triangular elements are used to discretize the slope in our simulations.
The binary search algorithm illustrated in Fig. 5 is employed to calculate the FOS. The critical failure state of a slope is determined when the optimization solver is not feasible or the maximum incremental displacement is higher than a given threshold. An initial range [RF1 0.2, RF2 10] is set to trigger the search process, and when the tolerance |RF2 − RF1|/RF1 < 10 −3 is achieved, the calculated reduction factor RF1 is the FOS of the slope.
Comparative studies regarding the FOS obtained from our SOCP-FEM analysis, the limit equilibrium method, and the traditional finite element method are presented in Table 1, revealing that satisfactory agreement is achieved, which demonstrates the correctness of the proposed SOCP-FEM for slope stability analysis.

post-failure Analysis by PFEM
In the framework of the PFEM, not only the FOS of a slope can be determined but also the post-failure processes and final runout distances for different reduction factors. This can be achieved through the dynamic analysis shown in Sect. 3.2. To illustrate  Table 1 Factor of safety determined by limit equilibrium method (LEM), strength reduction method with nonassociated flow rule (SRM1) and associated flow rule (SRM2), and SOCP-FEM analysis with nonassociated flow rule (SOCP1) and associated flow rule (SOCP2). Diff1 is the percentage difference between the FOSs obtained from SRM1 and SOCP1 in absolute value, while Diff2 is the same but between SRM2 and SOCP2 The LEM, SRM1, and SRM2 data are from Cheng et al. (2007) this capability, case 12 in Table 1 is reanalyzed with the nonassociated flow rule. Four different reduction factors are used, and the corresponding final deposits from the simulations are illustrated in Fig. 6. As displayed, for RF 2.24 and 2.3, clear sliding surfaces are identified from the PFEM simulation, which coincide with the slip surfaces determined by other Fig. 6 Final deposits obtained from PFEM simulation for reduction factor (RF) of a 2.24, b 2.3, c 2.4, and d 2.9. Contours refer to displacement of soils (unit: meters); slip surfaces from the analysis by SRM1, SRM2, and LEM extracted from Cheng et al. (2007) approaches (Fig. 6a, b). However, soil elements in these two cases experience very small displacements. This is because these reduction factors are very close to the factor of safety of the slope, implying that the slope will tend to become stable again after limited deformation. For RF 2.4 (Fig. 6c), the deformation is well identified and a maximum displacement of around 0.7 m is experienced by soil nodes along the slip surface. Further increase of the reduction factor leads to greater deformations; For example, when RF 2.9, the sliding mass has a moving distance exceeding 2 m, as well predicted by our approach.

Landslide Propagation
Computing the propagation of landslides requires an accurate description of the rheological behavior of the involved geomaterials. In geoscience, the propagation process of a landslide is usually simulated by solving equations describing the fluid behavior of the sliding mass based on Eulerian approaches. In the present PFEM framework, more sophisticated material constitutive models can be implemented for landslide propagation. It is of interest to compare the simulation results of a real landslide case obtained from the present PFEM method and those from depth-averaged models.
To this end, a historical event, i.e., the Tangjiashan mass failure in Sichuan Province, China, is addressed in this section (Fig. 7). The failure was triggered by the 2008 Wenchuan earthquake, and the sliding mass crashed rapidly into the Jianjiang River, causing a landslide dam with a death toll of up to 84 (Hu et al. 2009;Xu et al. 2009). The slide moved along a complex topography with maximum slope angle of 40 • , and the motion lasted around 30 s according to Hu et al. (2009). The data about the landslide body and the deposit are extracted from the field survey conducted by For the time integration, the backward Euler scheme, i.e., θ 1 θ 2 1, is used in the PFEM simulation. In addition to the PFEM, the landslide motion was also simulated by using the depth-averaged model (Xia and Liang 2018). The depth-averaged equations (DAEs) are derived from the Navier-Stokes equations by using the long-wave approximation and the Mohr-Coulomb rheology law (Savage and Hutter 1989) with topography modifications suitable for a global Cartesian coordinate system. The numerical code was developed according to the Nessyahu-Tadmor scheme (Nessyahu and Tadmor 1990;Tai et al. 2002) and has been validated with typical benchmarks and a historical landslide event (Wang et al. 2019). In contrast to the present PFEM model, the adopted depth-averaged model only considers the basal Mohr-Coulomb friction law, ignoring the effects coming from the internal friction angle and the cohesion of the sliding mass. This means that the material parameters required for the depth-averaged analysis consist of the density ρ 2000 kg/m 3 and the frictional coefficient of the basal surface μ 0.27. In the simulation by the depth-average model, the common used earth pressure coefficient (Savage and Hutter 1989) is assumed to be unity, like in other recent works (e.g., Russell et al. 2019).
The positions of the sliding mass obtained from the PFEM analysis at four different time instants are shown in Fig. 8, where the velocity magnitude contour is depicted and the dashed line represents the profile obtained from the depth-averaged model. At t 6 s (Fig. 8a), the landslide moves as a whole with velocity of around 20 m/s while the depth-averaged model provides a slightly faster motion. Later, velocity distinctions  (Fig. 8b). Impeded by the antislope, the front of the landslide decelerates rapidly while the rear of the landslide still moves with a velocity greater than 25 m/s (Fig. 8c). The sliding mass stops at t 24 s in the PFEM simulation, as shown in Fig. 8d, whereas the landslide process lasts longer in the depth-averaged simulation. In particular, the rear sector pushes the front mass and then moves back in the depth-averaged simulation, implying that the sliding mass behaves more like a fluid if the DAEs are used. Figure 9 compares the topographies from the observational data (Hu et al. 2009) and the simulation by SPH (Huang et al. 2012), the present PFEM, and the DAEs. All simulation results agree well with the observation data. The mass variations (represented by the ratio of the current total mass to the initial mass of the sliding body) of the PFEM simulation are shown in Fig. 9b, where a maximum of 1.5% mass change is observed. It is thus concluded that, although reidentification of computational domains is conducted in the PFEM, mass conservation can still be maintained well. Considering the landslide as a whole, the mean velocity profile of the landslide is shown in Fig. 9c. In the PFEM simulation, the landslide accelerates to around 26 m/s in the first 10 s then decelerates to still. In the depth-averaged simulation, the main mean velocity profile is quite similar, while an additional velocity profile is produced due to the back motion of the landslide body. In general, the present two simulations are nearly consistent in capturing the main motion of the Tangjiashan landslide. Nevertheless, it should be emphasized that the DAEs cannot be used to predict the triggering of the landslide (for example, slope stability analysis) due to its nature, whereas the present optimization-based PFEM is capable of both failure and post-failure predictions.

Conclusions
This paper focuses on the numerical implementation of a version of PFEM in which the finite element formulation is transformed into a second-order cone programming (SOCP-FEM), and on its application to failure and post-failure analysis of landslides.
Based on the Hellinger-Reissner variational principle, the finite element formulation for static/dynamic elastoplastic analyses with interactions between deformable bodies and a rigid surface can be cast into equivalent optimization programs and resolved efficiently using available advanced optimization engines. This paper provides a detailed introduction to the numerical implementation of the optimizationbased finite element method, particularly the SOCP-FEM. It is shown that the implementation can be achieved straightforwardly by constructing the corresponding matrices, significantly releasing the researcher from programming. The present details can serve as a reference for those who are willing to develop their own code with modifications for application to different geoscience problems.
Meanwhile, the version of PFEM, which is a mixture of the SOCP-FEM and the particle technique, is applied to the modeling of landslides with detailed comparative studies of both failure and post-failure analyses. Associated with the shear strength reduction technique and the binary search algorithm, the method is tested against the stability analysis benchmark published by Cheng et al. (2007). Additionally, it is shown that the present model is capable of simulating both the failure and postfailure stages of the landslide in a single simulation. Moreover, the propagation of the historical event of the 2008 Tangjiashan landslide, China is investigated by the PFEM model and a depth-averaged model (Xia and Liang 2018). Overall, satisfactory agreement is obtained between the results from the PFEM simulation and the depth-averaged modeling, demonstrating the capability of the PFEM for modeling landslide propagation, although the sliding mass behaves more like a fluid according to the depth-averaged simulation results.