Modelling of Critical Acceleration for Regional Seismic Landslide Hazard Assessments by Finite Element Limit Analysis

The critical acceleration model plays an important role in seismic slope stability and determines the predictive accuracy of regional seismic landslide hazard assessments. Recently, the critical acceleration model based on the limit equilibrium method has been used to evaluate the seismic stability of regional slopes. However, when the Hoek-Brown failure criterion is used to evaluate the seismic stability of slopes with angles greater than 60°, the results obtained is unconservative by limit equilibrium method. Therefore, based on the simulation of a typical slope model with finite element limit analysis, prediction equations of the critical acceleration are established. The corresponding results are compared with the prediction results from the limit equilibrium method. This comparison shows that the proposed critical acceleration model has higher predictive accuracy than the limit equilibrium method, especially when the Hoek-Brown failure criterion is used to evaluate the slopes with angles greater than 60°. The proposed model is applicable to the global scope and can be effectively applied to regional seismic landslide hazard assessments.


INTRODUCTION
Earthquakes have induced a large number of landslides, which buried towns and farmland and damaged infrastructure, causing substantial economic losses and casualties in the disaster area (Qi et al., 2010;Xu et al., 2015;Chen et al., 2018;Fan et al., 2018). Therefore, seismic landslides have been widely studied by geoscience experts. Relevant researchers have made beneficial exploration in the evaluation of slope stability and deformation (Yang et al., 2019a;Yang et al., 2019b;Wu et al., 2020;Yang et al., 2021), the compilation of landslide database (Xu et al., 2012; and landslide hazard assessment (Caccavale et al., 2017;Chen et al., 2018;Wang et al., 2018;, which provided important theoretical scientific support for exploring the mechanisms driving seismic landslides and reducing seismic landslide disasters. Regional seismic landslide hazard assessment is one of the most important research directions, and the results from these assessments can provide a reference for reconstructing post-earthquake disaster areas and predicting future earthquakeinduced landslides. Traditionally, there are two quantitative methods to evaluate the hazards of regional seismic landslides, including statistical analysis methods and mechanical model methods based on physical slope models. The statistical analysis method mainly uses seismic landslides that have occurred, establishes a mathematical statistical model of these landslides and the related landslide influence factors, and then applies this model to the whole earthquake area. The advantages of this method are that it is based on actual landslides and that the evaluation results are objective; however, this method lacks an understanding of the mechanism of seismic landslides (Xu et al., 2012;. In the mechanical model method, Newmark's sliding block analysis (Newmark, 1965) (referred to hereinafter as the Newmark method) is commonly used to estimate regional seismic landslide hazards (Jibson et al., 2000;Jibson, 2007;Rodríguez-Peces et al., 2011;Bozzano et al., 2013;Dreyfus et al., 2013;Rodríguez-Peces et al., 2014;Caccavale et al., 2017;Liu et al., 2017;Shinoda and Miyata 2017;Chen et al., 2018;Wang et al., 2018;. In the Newmark method, the slope model should be assumed to be an infinite plane (i.e., an infinite slope), and the failure depth of the slope should be known in advance (usually assumed to be less than 3 m) (Jibson et al., 2000). However, according to previous studies (Wilson and Keefer 1983;Keefer 1984;Huang et al., 2011;Wartman et al., 2013;, seismic landslides are widely developed in soil and rock units, and there are many failure modes. This shows that the assumption of an infinite slope may be too simplistic in many cases, thereby underestimating the risk area of sliding. Some scholars have realized the limitations of the Newmark method in regional seismic landslide hazard assessment and put forward measures to improve the predictive accuracy of the Newmark method (Shinoda and Miyata 2017;Jin et al., 2019;Zang et al., 2020). Although the predictive accuracy of the Newmark method has been improved through some measures (e.g., modified strength parameters), the limitations of the Newmark method are still difficult to overcome (Shinoda and Miyata 2017;Jin et al., 2019;Zang et al., 2020). Therefore, it is urgent and important to develop new mechanical model methods for regional seismic landslide hazard assessments.
Some scholars have recently proposed new mechanical models and used these models to carry out regional seismic landslide hazard assessments (Saade et al., 2016;Yuqiao et al., 2019;Tsai et al., 2019). Compared with the traditional Newmark model, the mechanical models proposed by these scholars have made the following improvements. 1) The modified models are not limited to a single shallow sliding mode. 2) The failure depth of the slope does not need to be assumed. 3) The strength parameters do not need to be reduced to improve the predictive accuracy of the model. 4) Moreover, the modified models have obvious advantages in predicting the potential sliding area. Although these mechanical models provide a new perspective for regional earthquake landslide hazard assessments, they are still not sufficiently comprehensive as they have the following shortcomings. 1) For Saade et al., 2016 andYuqiao et al., 2019, obtaining the critical acceleration of the slope with the limit equilibrium method (LEM) and the finite element method (FEM) requires many iterative calculations, which are estimated by trial and error and may therefore affect the calculation efficiency. 2) Saade et al., 2016, the safety factor and critical acceleration of the slope calculated by the Hoek-Brown (HB) criterion will be overestimated for steep slopes (Li et al., 2009;Li et al., 2019). 3) For Yuqiao et al., 2019, the mechanical model requires a slope angle <45°. 4) For Tsai et al., 2019, 12 equations and six steps are needed to evaluate the permanent displacement of landslides caused by earthquakes, which makes it difficult for engineers to quickly evaluate earthquake-induced landslide hazards.
A modified model of the critical acceleration for regional seismic landslide hazard assessments is proposed based on the work of Saade et al. (2016) to overcome the limitations of the above mechanical models. In this model, finite element limit analysis (FELA) is used to build a simplified slope model, and the material strength parameter and slope angle are taken as variables to conduct simulations of the slopes. Through two-step regression analysis, the prediction equations for the strength parameters, slope angle and critical acceleration are obtained. The predicted results of the proposed model are compared with those of the Newmark method and the LEM (Saade et al., 2016). The validity of the proposed model is verified through this comparison. In addition, the proposed model can be effectively applied to regional seismic landslide hazard assessments.  Saade et al. (2016) conducted a series of numerical simulations to establish the predicted equations of linear function between the slope angle β, the unit weight γ, strength parameters and the critical acceleration (A c ). In addition, Saade et al. (2016) used the Mohr-Coulomb (MC) failure criterion and the HB failure criterion to calculate the critical acceleration when the slope angle <45°and when the slope angle ≥45°, respectively. The prediction equations of two linear functions for critical acceleration is as follows: Where A c(MC) is the critical acceleration for the MC failure criterion, c is the cohesion, γ is the unit weight, H is the slope height, C1 (MC) is the coefficient related to slope angle β and friction angle φ for the MC failure criterion, and C2 (MC) is the coefficient related to slope angle β for the MC failure criterion, A c(HB) is the critical acceleration for the HB failure criterion, σ ci is unconfined compressive strength, and C1 (HB) and C2 (HB) are the coefficient related to slope angle β for the HB failure criterion.

FELA
In this study, FELA can be used to obtain the collapse load of slope (Lyamin et al., 2005;Sloan 2013;Ali et al., 2016). By establishing strict upper and lower bounds with a certain level of error, it is possible to accurately estimate the collapse load. For a given problem, L is the lower limit of the collapse load and U is the upper limit of the collapse load. Then, the estimate of the exact solution is simply taken as the average of the upper and lower limits: In addition, the exact solution is denoted e. Since L ≤ e ≤ U, the following expressions can be derived: where ε is the relative error. Previous experience shows that when appropriate upper and lower bounds are used, the actual error of the average value is usually less than the estimated value. Hence, the average value is an excellent estimate of the exact solution. Therefore, we calculate the upper-and lowerlimit solutions of the critical acceleration in this paper. The average bound solution of the critical acceleration is obtained according to the upper-and lower-limit solutions of the critical acceleration. The critical accelerations mentioned hereinafter are all average bound solutions. Detailed information about the technologies and solution algorithms used in FELA can be found in the literature (Sloan 1988;Sloan 1989;Lyamin and Sloan 2002a;Lyamin and Sloan 2002b;Lyamin et al., 2005;Sloan 2013;Krabbenhoft and Lyamin, 2015;Li et al., 2019).

Calculation Settings and Input Parameters
At present, the design using pseudostatic analysis is usually based on the horizontal seismic coefficient (k h ). The critical acceleration mentioned hereafter corresponds to the critical horizontal seismic coefficient (k c ) multiplied by the gravitational acceleration. FELA can directly solve the ultimate load of a slope subjected to an earthquake (such as the critical horizontal acceleration) without requiring iterative calculations, thereby significantly reducing the computational effort (Sloan 2013;Utili and Abd 2016). A typical slope model (slope height H is 30 m; unit weight γ is 22 kN/m 3 ) is applied for parametric study (Figure 1). To obtain the minimum critical acceleration for the MC failure criterion, DH is set as 100 m when the slope height is 30 m (Loukidis et al., 2003;Saade et al., 2016). Besides, to eliminate the influence of the lateral boundary on the sliding surface of the slope, the lateral boundary is set away from the slope surface. The bottom boundary of the slope model is fixed. The numerical model uses the adaptive mesh generation method, in which the number of adaptive iterations is 3, the initial number of elements is 2,000, and the maximum number of elements is 4,000. The control variable of adaptive meshing is shear dissipation. The mesh refinement factor is 0.25, and the mesh coarsening factor is 1.50.
Based on Saade et al., 2016, we also use the MC failure criterion and HB failure criterion to calculate the critical acceleration of the slope when the slope is less than 45°and when the slope is greater than or equal to 45°, respectively. Table 1 and Table 2 show the slope angle and strength parameters for the two failure criteria. For each combination described in Table 1 and Table 2, two slope models (including the upper-limit solution and lower-limit solution) are established. While keeping the gravity load constant, the horizontal load multiplier is maximized (lowerbound simulation) and minimized (upper-bound simulation), and the critical horizontal acceleration is calculated. The results of finite element discretization of the upper-and lower-limit solutions of the MC failure criterion and HB failure criterion are shown in Figure 2.

HB Failure Criterion
Based on Saade et al., 2016, we also use the same linear function, as shown in Eq 2. Figure 5A shows the relationship between the critical acceleration and non-dimensional term (σ ci /γH) when the GSI is 20 and m i is 25. For each GSI (20,30,40,50) and m i (5,10,15,20,25), similar curves are generated. According to different functions of slope angle β, the C1 (HB) and C2 (HB) are obtained respectively, as shown in Figures 5B,C. In addition, the optimal fit equations of the C1 (HB) and C2 (HB) for other GSI-m i conditions are considered in this study, as shown in Table 4. The C1 (HB) and C2 (HB) of any other conditions can be calculated by interpolation method. The fit index (R2 adjusted) is greater than 0.92 for any conditions. Figure 6 shows the predicted and calculated results of critical acceleration using the LEM, the Newmark method and FELA when the GSI is 20, m i is 25, σ ci is 5 MPa and 45°≤β ≤ 80°. As shown in Figure 6, the results of the Newmark method substantially overestimate the critical acceleration: A cNewmark = 1.8819A c(HB) +0.0006 (Saade et al., 2016). However, the results when using FELA are in good agreement with the critical acceleration: A cFELA = 1.02454A c(HB) −0.00134. This is very consistent with the results reported by Saade et al. (2016) using the LEM: A cLEM = 1.0218A c(HB) +0.0006.  Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 10 | Article 830371 5

Comparison of the LEM and FELA Results
From Section 3.2, specific strength parameters are selected to compare the critical acceleration values predicted by the LEM and FELA. This shows that the two methods have good consistency. However, to extend the predicted regression equation to additional conditions, the results of other strength parameters and geometric parameters need to be further tested to assess their validity, especially when used the HB failure criterion. Notably, since the FELA results are closer to the real solution than the LEM results, the calculated critical accelerations adopt the results from the FELA simulation (Li et al., 2009;Sloan 2013;Li et al., 2019) as the Y-axis in Figure 7 and Figure 8. Figure 7 compares the results of predicted A c and FELA simulation for c = 20kPa, and φ = 25°and φ = 45°at variable β when use MC failure criterion. This shows that the regression equations from the LEM and FELA are in good agreement under different strength parameters. This is consistent with the conclusion of Loukidis et al. (2003). Figure 8 compares the Ac values predicted and calculated with FELA for m i = 25, σ ci = 5 and 25 MPa, 45°≤β ≤ 80°and GSI = 20, 30, 40 and 50 when using the HB failure criterion. Figure 8 shows that the A c(HB) values predicted and calculated with FELA are very consistent under any combination of parameters (the relationship is almost one to one). However, the A c(HB) values predicted with the LEM and FELA are not consistent, and the difference between these values increases gradually as the strength parameter and slope angle increase. Taking     Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 10 | Article 830371 6 m i = 25, and σ ci = 25 MPa, the regression equation from Saade et al. (2016) is y = 1.9568x-1.1260. When the slope angle β is 70°, the predicted critical acceleration obtained by the LEM regression equation is 29.76-86.58% higher than that obtained by the FELA regression equation. In addition, the FELA regression equation can accurately estimate the critical acceleration of slopes with 60°<β ≤ 80°, which is difficult to determine with the regression equation from Saade et al. (2016).

EXAMPLES
From Sections 3.2 and 3.3, we compare the A c values predicted and calculated with FELA when H is 30 m and γ is 22 kN/m 3 . To further test the influence of the slope height and unit weight on the regression equation, Table 5 and Table 6 compares the A c values predicted and calculated with FELA for different slope heights (H is 10 and 50 m) and unit weights (γ is 18 kN/m 3 and 28 kN/m 3 ). From Table 5 and Table 6, at a certain slope height and unit weight, the larger the slope angle is, the larger the error in the critical acceleration predicted by the LEM. The critical acceleration predicted and calculated by FELA are in good agreement. Case Case 1. is a 10-m-tall slope with β = 45°composed of rock with GSI is 40, m i is 25, γ is 22 kN/m 3 , σ ci is 10 MPa.
1) The process for predicted critical acceleration of LEM from the regression coefficient (Saade et al., 2016) and Eq 2 is as follows: Predicted 2) The process for predicted critical acceleration of FELA from the 15th regression coefficient in Table 4 in combination with Eq 2 is as follows:  engineers. Every engineer may obtain different values for the material strength parameter in the same region, and there are also differences in the empirical sliding displacement models in each region. However, the critical acceleration model is an objective and unified factor that can serve global applications and plays an important role in the framework of earthquakeinduced permanent displacements. In addition, critical acceleration can not only be used as a criterion to evaluate the susceptibility of seismic landslides , but also be used to evaluate the hazard of seismic landslides by comparing with known peak ground acceleration (Chen et al., 2018). Therefore, it is of great significance to establish a reliable critical acceleration model for regional seismic landslide hazard assessments.
In previous studies, due to the inherent limitations of the traditional Newmark method and its modified versions, the calculated critical accelerations are substantially different from the numerical simulations. Compared with the Newmark method, the critical acceleration models from Saade et al. (2016), Yuqiao et al., 2019, Tsai et al. (2019 and this study offer obviously higher predictive accuracy, providing strong support for regional seismic landslide hazard assessments. In addition, in the regional scale analysis, the vertical acceleration, slope direction, groundwater conditions, water content, material degradation and rock slope disturbance coefficient may affect the critical acceleration model, which were not discussed in this study. Therefore, the development of a more comprehensive critical acceleration model should be the direction of future research. However, our model can directly estimate strength parameters from a good engineering geological description while still providing an accurate quantitative index, which can quickly measure the behaviour of a seismic slope at a regional scale. On the other hand, Li et al. (2008) found that the difference of safety factor from limit analyses and LEM is less than 4% for steep slopes when used the HB failure criterion. However, Li et al. (2008) reported this conclusion when the safety factor of the slope is near 1. Saade et al. (2016) did not seem to notice this phenomenon. In addition, for lower GSI values (GSI ≤ 50), the stability numbers of slope (the stability number is equal to σ ci /γHF, wherein F is the safety factor of slope) obtained using the LEM are 18-30% lower than the average bound solutions obtained with FELA (Li et al., 2009). Moreover, the σ ci /γHF from the LEM are close to the results of FELA when β ≤ 60°. However, the difference of the σ ci /γHF between FELA and the LEM tends to be substantial when β ≥ 75° (Li et al., 2019). Therefore, according to the research of Li et al. (2009Li et al. ( , 2019, we revised the model from Saade et al. (2016).  Our model improves the limitations of previous mechanical models for regional seismic landslide hazard assessments. Herein, FELA is used to calculate the average bound solution of the critical acceleration of the simplified slope model. The proposed model of critical acceleration is not limited to failures of shallow sliding slopes nor does it require the failure depth to be assumed in advance. The MC failure criterion and the HB failure criterion are used to simulate the seismic behaviour of soil and rock slopes. The prediction equations for the critical horizontal acceleration, material strength and slope gradient are established. The corresponding results are compared with the prediction results from the Newmark sliding analysis and the LEM. The results show that the predictive accuracy of FELA and the LEM is significantly higher than that of the Newmark sliding analysis. When the MC failure criterion is used, the prediction results from the FELA and LEM regression equations are consistent. However, the prediction results from the FELA and LEM regression equations are significantly different, especially when the slope gradient is greater than 60°and uses HB failure criterion. In general, the predictive accuracy of our model is higher than that of the LEM. Therefore, the proposed model can be effectively applied to regional seismic landslide hazard assessments.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material further inquiries can be directed to the corresponding authors.