An Infinite Slope Model Considering unloading joints for Spatial Evaluation of Coseismic Landslide Hazard: A Case Study of the 2013 Lushan Earthquake

1. Introduction

Earthquakes are one of the primary triggers for landslides, often resulting in widespread coseismic landslides within specific regions [1]. In recent years, there has been a rise in strong earthquakes globally, leading to frequent coseismic landslide hazards and causing significant loss of life and property damage. Consequently, the evaluation of regional coseismic landslide hazards holds paramount importance for pre-earthquake urban planning and post-earthquake reconstruction efforts.

Early attempts to assess the coseismic stability of slopes utilized the pseudo-static analysis method introduced by Terzaghi [2] and the finite-element modeling technique pioneered by Clough and Chopra [3]. Additionally, Newmark [4] introduced a straightforward and practical method, which is still widely employed today, to estimate the coseismic permanent displacements of slopes [5]. Franklin and Chang [6] conducted permanent displacement calculations on earth and rock dam slopes using the Newmark method, incorporating 179 actual strong motion records and 10 artificially fitted accelerograms. This extensive dataset greatly expanded upon the information presented in Newmark's 1965 report, significantly bolstering the credibility of the finite sliding displacement methodology. To enhance the accuracy of estimating earthquake-induced permanent displacement, decoupled analysis was introduced by Makdisi and Seed [7]. The decoupled analysis accounts for the deformation of the sliding block during the sliding process. In the decoupled method, the block's deformation is calculated independently, without considering sliding. Subsequently, this deformation is treated as a new acceleration history to be applied to the sliding system. Bray and Rathje [8] introduced a streamlined seismic displacement method by utilizing a fully nonlinear decoupled one-dimensional dynamic analysis technique in conjunction with Newmark's rigid-plastic sliding block analysis. Building upon this work, Bary and Travasarou [9] proposed a simplified approach to predict earthquake-induced permanent displacement. Their method, based on a nonlinear fully coupled sliding block model, incorporates critical acceleration, dominant sliding block period, and seismic spectral acceleration, providing a more efficient and accurate framework for seismic displacement prediction. Rathje and Bray [10,11] introduced the coupled analysis that incorporate the deformation of the sliding block during the sliding process. In the coupled method, deformation and displacement are considered simultaneously, offering a more realistic portrayal of slope response under seismic forces.

Jibson et al. [12,13] conducted an extensive investigation concerning the Northridge earthquake of 1994, which struck in close proximity to Los Angeles, California. Their research focused on the Oak Mountain region, situated to the north of the earthquake's epicenter. In a systematic manner, they elucidated the computational process of the Newmark method for the regional assessment of seismic landslide susceptibility. Additionally, they presented a comprehensive seismic landslide susceptibility zoning map for the study area. This pioneering work significantly advanced the utilization and refinement of the Newmark method within the realm of regional seismic landslide susceptibility analysis. Subsequently, numerous scholars have adopted and applied the Newmark method in their assessments of regional seismic landslide susceptibility. Yuan et al. [14] applied the basic methodology outlined by Jibson [12,13] to derive the Newmark displacement of landslides based on strong ground-motion recordings during the 2013 Lushan Ms 7.0 earthquake. Different empirical equations with new fitting coefficients for estimating Newmark displacement are developed for comparative analysis. Zhang et al. [15] have introduced a novel model for evaluating seismic slope stability, building upon the foundation established by Jibson [16] in the Newmark model. In this enhanced model, the influence of velocity pulse seismic motion on permanent slope displacement was meticulously considered. The study meticulously incorporates the effects of near-fault pulse characteristics, thereby presenting an advanced framework for assessing seismic slope stability under the specific influence of these pulses. Zhang et al. [17] delved into the nuanced topic of force direction applied to the sliding block within the framework of the Newmark model. Their study advocates for the adoption of a pragmatic approach by suggesting a more realistic horizontal base direction for the force exerted on the sliding block.

Natural slopes often develop a series of shallow unloading joints parallel to the surface due to valley incisions [18,19]. Research indicates that during strong earthquakes, rock slopes exhibit collapsing and sliding failures along these shallow unloading joints, with approximately 90% of coseismic landslides being shallow falls and slides [20,21,22,23]. Unstable rock blocks are frequently fragmented and mobilized along these joints [24]. However, prior studies have not sufficiently considered the shear strength of these rock joints in the assessment of coseismic landslides. In addressing this gap, Zang et al. [25] introduced the Barton model [26] into the Newmark analysis to enhance the estimation of slopes' dynamic stability. This method was initially employed to assess the hazards posed by coseismic landslides, utilizing data from the 2014 Ludian earthquake in Yunnan Province, southwestern China. The seismogenic fault associated with the 2014 Ludian earthquake is identified as the Baogunao-Xiaohe fault, characterized by left-lateral strike-slip movement [27]. In this study, we critically review the initial methodology and refine it through recalibration, utilizing the coseismic landslide dataset from the 2013 Lushan earthquake. The Lushan earthquake is categorized as a blind reverse-fault events [28].

2. Study Area

The epicenter of the 2013 $M_{\mathrm{w}}$ 6.6 Lushan earthquake is located at the southern segment of the Longmenshan fault zone, on the eastern margin of the Tibetan Plateau [29,30]. A rectangular containing dense concentrations of earthquake-induced landslides was chosen as the study area (Figure 1). The area is located about 12 km northwest of the epicenter and lies on the border of two earthquake-hit counties, i.e., Baoxing county and Lushan county in Sichuan Province, China, covering 138 km² (Figure 1). The study area is characterized by deeply incised valleys and high mountains, with an elevation ranged from 1013 to 3308 m above sea level (Figure 1). Exposures of the geological units have an age that varies from the Sinian to the Triassic, including carbonate, diorite, dolomite, limestone, phyllite, sandstone and shale (Figure 2).

There have been 17 earthquakes with magnitudes ≥ 4.7 in the history in the Lushan earthquake zone [31]. Strong tectonic activity provides favorable conditions that prone to landslides. The Lushan earthquake has triggered thousands of landslides [32]. An inventory of 308 coseismic landslides was compiled by comparisons between pre-earthquake satellite images and post-earthquake aerial photographs in the study area. The total area of these interpreted landslides is 5 km². A majority of landslides developed along the deeply incised river valleys.

3. Methodology

Newmark in 1965 [33] proposed an infinite slope model for the analysis of the dynamic stability of slopes. Newmark’s method [33] simulates a landslide as a rigid plastic friction block lying on an inclined plane with a known critical acceleration (Figure 3). It calculates the cumulative permanent displacement of the block relative to its base as it is subjected to the effects of an earthquake acceleration-time history [12,13]. The Newmark displacement is a valuable index for evaluating the dynamic performance of slopes.

Once an acceleration-time history has been selected, those portions of the record that exceed the critical acceleration (Figure 4a) are integrated once to derive a velocity profile (Figure 4b); the velocity-time history is then integrated a second time to obtain the cumulative displacement of the block (Figure 4c) [12,13,34].

From Figure 4, we can see that the critical acceleration provides a base for predicting the cumulative permanent displacement of a slope. Newmark [33] showed that the critical acceleration of a potential landslide block can be expressed as a simple function of the static factor of safety and the geometry of the slope [12,13]:

$$a_{\mathrm{c}}=\left(F_{\mathrm{S}}-1\right)g\mathrm{s}\mathrm{i}\mathrm{n}\alpha ,$$

(1)

where $a_{\mathrm{c}}$ is the critical acceleration in terms of g, $F_{\mathrm{S}}$ is the static factor of safety, and $\alpha$ is the angle from the horizontal at which the center of the slide block moves when displacement first occurs [12,13]. For a planar slip surface parallel to the slope, this angle generally approximates to the angle of the slope [34].

Natural rock slopes are often cut by a group of shallow unloading joints due to valley incisions, forming a unloading zone on the surface of slopes [18,19]. Slopes behave as collapsing and sliding failures of shallow unloading joints under strong earthquakes [34], and 90% of coseismic landslides are shallow falls and slides [20,22,23,35]. Unstable rock blocks are often activated along the rock joints during seismic shaking (Figure 3). Therefore, the static factor of safety is related to the shear strength of these rock joints [34].

Zang et al. [34] derived the static factor of safety of a slope based on the Barton [26] shear strength criterion for rock joints in another research:

$$F_{\mathrm{S}}=\frac{\mathrm{R}\mathrm{e}\mathrm{s}\mathrm{i}\mathrm{s}\mathrm{t}\mathrm{i}\mathrm{n}\mathrm{g}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{c}\mathrm{e}}{\mathrm{D}\mathrm{r}\mathrm{i}\mathrm{v}\mathrm{i}\mathrm{n}\mathrm{g}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{c}\mathrm{e}}=\frac{\tau }{\gamma t\mathrm{s}\mathrm{i}\mathrm{n}\alpha }=\frac{{\sigma }_{\mathrm{n}}\mathrm{t}\mathrm{a}\mathrm{n}[{\mathrm{J}\mathrm{R}\mathrm{C}}_n{\log }_{10}\left(\frac{{\mathrm{J}\mathrm{C}\mathrm{S}}_n}{{\sigma }_{\mathrm{n}}}\right)+{\varphi }_{\mathrm{b}}]}{\gamma t\mathrm{s}\mathrm{i}\mathrm{n}\alpha }=\frac{\gamma t\mathrm{c}\mathrm{o}\mathrm{s}\alpha \mathrm{t}\mathrm{a}\mathrm{n}[{\mathrm{J}\mathrm{R}\mathrm{C}}_n{\log }_{10}\left(\frac{{\mathrm{J}\mathrm{C}\mathrm{S}}_n}{\gamma t\mathrm{c}\mathrm{o}\mathrm{s}\alpha }\right)+{\varphi }_{\mathrm{b}}]}{\gamma t\mathrm{s}\mathrm{i}\mathrm{n}\alpha }=\frac{\mathrm{t}\mathrm{a}\mathrm{n}[{\mathrm{J}\mathrm{R}\mathrm{C}}_n{\log }_{10}\left(\frac{{\mathrm{J}\mathrm{C}\mathrm{S}}_n}{\gamma t\mathrm{c}\mathrm{o}\mathrm{s}\alpha }\right)+{\varphi }_{\mathrm{b}}]}{\mathrm{t}\mathrm{a}\mathrm{n}\alpha },$$

(2)

where $\tau$ is the peak shear strength of the rock joint, $\gamma$ is the unit weight of the rock mass, $t$ is the thickness of the failure rock block, ${\sigma }_{\mathrm{n}}$ is the effective normal stress, ${\mathrm{J}\mathrm{R}\mathrm{C}}_n$ is the joint roughness coefficient in the in-situ scale, ${\mathrm{J}\mathrm{C}\mathrm{S}}_n$ is the joint wall compressive strength in the in-situ scale, and ${\varphi }_{\mathrm{b}}$ is the basic friction angle.

Considering the impact of joint size, the ${\mathrm{J}\mathrm{R}\mathrm{C}}_n$ and ${\mathrm{J}\mathrm{C}\mathrm{S}}_n$ can be defined as [36]:

$${\mathrm{J}\mathrm{R}\mathrm{C}}_n={\mathrm{J}\mathrm{R}\mathrm{C}}_0{\left(\frac{L_n}{L_0}\right)}^{-0.02{\mathrm{J}\mathrm{R}\mathrm{C}}_0},$$

(3)

$${\mathrm{J}\mathrm{C}\mathrm{S}}_n={{\mathrm{J}\mathrm{C}\mathrm{S}}_0\left(\frac{L_n}{L_0}\right)}^{-0.03{\mathrm{J}\mathrm{R}\mathrm{C}}_0},$$

(4)

where the nomenclature adopted incorporates (0) and (n) for values of the laboratory scale and the in-situ scale, respectively.

It is impractical to conduct a rigorous Newmark’ method during the regional analysis. Therefore, researchers have proposed different empirical regressions to estimate Newmark displacement as a function of the critical acceleration and ground motion parameters [16,37,38,39,40,41]. In this study, we chose a vector model that developed according to more than 2000 strong motion records [40]:

$$\mathrm{l}\mathrm{n}D=4.89-4.85\left(\frac{a_{\mathrm{c}}}{\mathrm{P}\mathrm{G}\mathrm{A}}\right)-19.64{\left(\frac{a_{\mathrm{c}}}{\mathrm{P}\mathrm{G}\mathrm{A}}\right)}^2+42.49{\left(\frac{a_{\mathrm{c}}}{\mathrm{P}\mathrm{G}\mathrm{A}}\right)}^3-29.06{\left(\frac{a_{\mathrm{c}}}{\mathrm{P}\mathrm{G}\mathrm{A}}\right)}^4+0.72\mathrm{l}\mathrm{n}\left(\mathrm{P}\mathrm{G}\mathrm{A}\right)+0.89\left(M_{\mathrm{w}}-6\right),$$

(5)

where $D$ is the predicted Newmark displacement, $\mathrm{P}\mathrm{G}\mathrm{A}$ is the peak ground acceleration, and $M_{\mathrm{w}}$ is the moment magnitude. $D$ is in centimeters, and $a_{\mathrm{c}}$ and $\mathrm{P}\mathrm{G}\mathrm{A}$ are in units of g.

Figure 5 is a flowchart showing the sequential steps of the spatial evaluation of coseismic landslides. After calculating the Newmark displacement, the predicted displacement were then compared with an inventory of landslides triggered by the Lushan earthquake to map the spatial variation of coseismic landslide hazards with the help of the certainty factor model [42]. The certainty factor model (CFM) is a model of inexact reasoning that created by Shortliffe and Buchanan [42] and improved by Heckerman [43], to explore the relationship between the predicted displacement and the occurrence of landslides. In the CFM, the certainty factor (CF) represents the net confidence in a hypothesis $H$ based on the evidence $E$ [43], ranging from -1 to 1. A value of -1 means a total lack of confidence, whereas a value of 1 means total confidence. Values greater than zero favor the hypothesis while those less than zero favor its negation. The probabilistic interpretation of CF can be expressed as follows [43]:

$$\mathrm{C}\mathrm{F}=\left\{\begin{array}{c}\frac{p\left(H|E\right)-p\left(H\right)}{p\left(H|E\right)\left[1-p\left(H\right)\right]}, p\left(H|E\right)>p\left(H\right)\\ \\ \frac{p\left(H|E\right)-p\left(H\right)}{p\left(H\right)\left[1-p\left(H|E\right)\right]}, p\left(H|E\right)<p\left(H\right)\end{array}\right.,$$

(6)

where CF is the certainty factor, $p\left(H|E\right)$ is the posterior probability that denotes the conditional probability for a posterior hypothesis that relies on evidence, and $p\left(H\right)$ is the prior probability without any evidence. For the spatial evaluation of coseismic landslides, $p\left(H|E\right)$ was defined as the proportion of the area of the landslide within a specific displacement area, and $p\left(H\right)$ was defined as the proportion of the landslide area within the entire study area [34]. In this way, values of CF represented the confidence level for occurrences of coseismic landslides. Positive values corresponded to an increase in the confidence level in slope failure while negative quantities corresponded to its negation [34].

4. Mapping of Coseismic Landside Hazard

In this section, each step involved in the mapping of coseismic landslide hazard is discussed in detail.

4.1. Static Factor of Safety Map

We selected a 30 m digital elevation model (DEM) to produce a slope map (Figure 6) by applying a basic slope algorithm, where the slop was identified as the steepest downhill descent among adjacent cells [44]. The derived slopes ranged from 0° to 75°. The northwestern and northeastern portions of the study area are high mountains and deeply incised valleys with slopes greater than 50°, and the central and southeastern parts of the area are relatively flat (Figure 6).

We assigned physical and mechanical parameters to the rock types using a digital geological map (Figure 2). The map was rasterized at 30 m grid spacing for following calculation after the parameter assignment. Studies have shown that ${\mathrm{J}\mathrm{R}\mathrm{C}}_0$ and ${\mathrm{J}\mathrm{C}\mathrm{S}}_0$ strongly depend on lithology [45,46,47,48,49,50,51,52,53,54,55]. Values of ${\mathrm{J}\mathrm{R}\mathrm{C}}_0$ and ${\mathrm{J}\mathrm{C}\mathrm{S}}_0$ assigned to each rock types were estimated based on test data from the references listed in Table 1. These values were obtained from laboratory tests with a sample length of $L_0$=100 mm. As we used a projected coordinate system for regional analysis and the grid spacing was 30 m, the engineering dimension, $L_n$, in a grid cell was equal to 30 m/ cos$\alpha$. We can calculate ${\mathrm{J}\mathrm{R}\mathrm{C}}_n$ and ${\mathrm{J}\mathrm{C}\mathrm{S}}_n$ by substituting ${\mathrm{J}\mathrm{R}\mathrm{C}}_0$, ${\mathrm{J}\mathrm{C}\mathrm{S}}_0$, $L_0$ and $L_n$ into Equations (3) and (4), respectively. And the spatial distribution of ${\mathrm{J}\mathrm{R}\mathrm{C}}_n$ and ${\mathrm{J}\mathrm{C}\mathrm{S}}_n$ are shown in Figure 7a,b, respectively. Other parameters such as $\gamma$ and ${\varphi }_{\mathrm{b}}$ are listed in Table 1.

For the sake of simplicity, the thickness of the failure rock block (Figure 3) $t$ was set to 3 m according to the filed investigation on typical slope failures of the Lushan earthquake. We produced the map of static factor of safety by combing the data layer of $\alpha$, ${\mathrm{J}\mathrm{R}\mathrm{C}}_n$, ${\mathrm{J}\mathrm{C}\mathrm{S}}_n$, $\gamma$ and ${\varphi }_{\mathrm{b}}$ in Equation (2). The calculated static factors of safety ranged from 0.2 to 133.8. Grid cells with static factors of safety less than 1 indicated that the slopes were statically unstable, but this did not necessarily mean that they were sliding under seismic shaking [34]. However, according to Equation (1), to avoid a negative value of the critical acceleration, we assigned a minimal static factor of safety of 1.01 for these cells. Keefer [56] showed that the minimum slope angle for earthquake-induced landslides was 5°. In addition, the calculated static factors of safety for the cells with a slope angle less than 5° were very high, and these cells did not produce a statistically significant sample in the analysis [34]. Therefore, slopes gentler than 5° were not analyzed. The adjusted static factors of safety range from 1 to 10.4, as shown in Figure 9. Statically unstable slopes with a static factor of safety less than 1.2 are mainly distributed in the northwestern, northeastern, and southeastern parts of the study area (Figure 8).

4.2. Critical Acceleration Map

After the adjustment of static factors of safety, we produced the map of critical acceleration using Equation (1) to combine the slope angle with the static factors of safety (Figure 9). The critical acceleration represents the intrinsic properties of slopes independent of any shaking scenario; thus, the critical acceleration map indicates the susceptibility of coseismic landslides [12,13]. A smaller value of critical acceleration denotes a lower intensity of shaking needed to overcome the stability of slope. Therefore, coseismic sliding is more likely to develop on the slope. From Figure 9, dynamically unstable slopes with a critical acceleration less than 0.3 g are mainly distributed in the northwestern part of the study area, and only a small portion is present in the southern part of the study area.

4.3. Predicted Dispalcement Map

We downloaded the contour map of peak ground acceleration (PGA) of the study area (Figure 10) from the United States Geological Survey. The Newmark displacement of each cell was calculated by employing Equation (5) to combine corresponding grid values of critical acceleration, PGA, and moment magnitude. The predicted displacements range from 0 to 76 cm (Figure 11). Slopes with a predicted displacement ranging between 0 cm and 4 cm are distributed in the northwestern portion of the study area, slopes with a predicted displacement ranging between 4 cm and 6 cm are distributed in the middle portion of the study area, slopes with a predicted displacement ranging between 6 cm and 15 cm are scattered in the area with a displacement between 4 cm and 6 cm, and slopes with a predicted displacement greater than 15 cm are only distributed in the southeastern portion of the study area. According to the statistical sizes of the areas, displacements less than 4 cm occupy 90.5 % of the study area, displacements between 4 cm and 6 cm occupy 5.5 % of the study area, and displacements between 6 cm and 15 cm occupy 3.6 % of the study area. Displacements greater than 15 cm occupy a very small area.

4.4. Coseismic Landslides Hazard Map

Jibson et al. [12,13] showed that Newmark displacements provide an index of the seismic performance of slopes but do not correspond directly to measurable slope movements in the field. Therefore, predicted displacements needed to be compared with an inventory of landslides triggered by earthquakes to be used in a predictive scenario of coseismic landslide hazards. We employed Equation (6) to produce a coseismic landslide hazard map in terms of certainty factors.

The Newmark displacement cells every 1 cm were grouped into bins, such that all cells with displacements between 0 cm and 1 cm were grouped into the first bin; those with displacements between 1 cm and 2 cm were grouped into the second bin, and so on [12,13,34]. The predicted displacements were grouped into 75 bins. For each bin, the proportion of cells occupied by landslide areas was calculated. This proportion represents the posterior probability of the bin according to Equation (6). The proportion of the entire landslide area within the entire study area was calculated to obtain the prior probability, which was same in each bin. We used Equation (6) to compute values of CF by combing the corresponding values of the posterior and prior probabilities. The calculated certainty factors ranged from −1 to 0.99. The CFM provides a necessary linkage between the Newmark displacements and confidence levels of coseismic landslide occurrence in the study area. We produced the coseismic landslide hazard map for the ground shaking scenario of the Lushan earthquake based on the values of CF, as shown in Figure 12. The landslides triggered by the Lushan earthquake were shown to demonstrate the goodness of fit for the predicted confidence levels of coseismic landslides.

5. Discussion

The predicted displacement is the cumulative sliding displacement of slopes for a given history of acceleration. As depicted in Figure 11, the predicted displacements range from 0 to 76 cm. Analysis of the statistical size of each displacement area indicates that 90.5% of the study area experiences displacements less than 4 cm, while 99.6% of the study area experiences displacements less than 15 cm. Displacements greater than 15 cm are observed in a very small area. Jibson et al. suggested that shallow, disrupted rock falls and rock slides in relatively brittle, poorly cemented sediments tend to collapse at small displacements [12,13]. Thus, the study area is more susceptible to shallow falls and slides. According to high-resolution aerial photographs and field investigations, most of the landslides triggered by the Lushan earthquake were rock falls and relatively shallow, disrupted slides, typically 1–5 m in depth [23,32,57]. Therefore, the proposed model is particularly useful for predicting the spatial distribution of the typical coseismic landslides in the study area.

Predicted displacements do not accurately reflect the actual slope movement in the field. Instead, modeled displacements serve as an indicator that can be correlated with field performance [12,13]. By using a function that incorporates Newmark displacement and CF, confidence levels of slope failure in the field can be estimated, and the spatial variation in coseismic landslides can be predicted under any set of ground-shaking conditions [12,13,34]. The number of Newmark displacement cells per 1 cm was uneven, so the predicted displacement cells were grouped into bins based on quantile statistics to obtain a more rational regression curve relating the predicted displacement and CF. The breakpoints were set at 0, 2, 4, 6, 15, 25, 44, 60, and 76, with an equal number of cells in each bin. For each bin, the proportion of cells in landslide areas was calculated, and the CF value of Newmark displacement was plotted as a dot. To fit the data, a modified Weibull [58] curve proposed by Zang et al. [8] was employed, with the following functional form:

$$\mathrm{C}\mathrm{F}=m\left[1-\exp \left(-a{D}^b\right)\right]-1,$$

(7)

where CF is the certainty factor, $(m-1)$ is the maximum CF value, $D$ is the predicted displacement, and $a$ and $b$ are the regression constants. The regression curve based on the Lushan data is:

$$\mathrm{C}\mathrm{F}=1.254\left[1-\exp \left(-0.669D^{0.682}\right)\right]-1.$$

(8)

The data fits the curve well with an R-squared value of 92%. The CF value can be used to directly predict the confidence level of slope failure as a function of Newmark displacement. Figure 13 shows, the CF value increases monotonically with increasing Newmark displacement. The CF value increases rapidly in the first 10 cm of Newmark displacement and then levels off abruptly after 15 cm range at a CF value of about 0.25.

Equation (8) takes the same forms as our previously published equation [34]. But the equation proposed in this study has different regression constants owing to a different data. Once calibrated, the curve and corresponding equation can be used in any set of ground-shaking conditions to predict the confidence level of slope failure as a function of predicted Newmark displacement [12,13]. However, a problem exposed by the earlier study should be noted. The maximum CF value is equal to $(m-1)$ not $m/2$ as indicated by Zang et al. [34]. The values of $m$, $a$, and $b$ in Equation (8) may vary in other regions if the lithology, topography, and ground motion conditions significantly differ from those observed in the study area. Additionally, in regions where the predominant type of slope failure differs, the shape of the fitting curve in Figure 13 would probably exhibit some variations [12,13]. If rock falls and rock slides in more brittle materials were predominant, the fitting curve would likely be more steep and might flatten out at a smaller maximum displacement value [12,13,34].

Figure 9 shows the distribution of critical acceleration in the study area. The critical acceleration describes seismic susceptibility of a slope. Newmark [33] showed that the critical acceleration is a function of static factor of safety and slope angle (Equation (1)). Compare Figure 6 with Figure 8 and Figure 9, respectively, we found that the modeling procedure is heavily slope-driven, which is consistent with the findings from Jibson et al. [12,13]. However, compared Figure 11 with Figure 10, we observed that the distribution pattern of Newmark displacements bears a resemblance to that of PGA. Areas with displacements between 2–4 cm, 4–6 cm, and 6–15 cm are distributed in bands in the study area, indicating the dominant role of seismic ground motion in triggering coseismic landslides during the Lushan earthquake.

In Section 3, it was mentioned that during seismic activity, unstable rock blocks tend to slide along shallow un-loading joints. The static factor of safety is closely related to the shear strength of these rock joints. The conventional Newmark analysis traditionally describes the shear strength of rocks using Coulomb's constants, such as friction angle and cohesion, but these values vary significantly between laboratory conditions with high normal stress and field conditions with low normal stress [26]. Additionally, these values are also dependent on the scale of the analysis[46]. To address these challenges, we introduced the Barton model into the Newmark analysis. This model enables us to predict the shear strength of rock joints and reduce the variability associated with Coulomb's constants. Furthermore, we considered the impact of scale effects by incorporating Equation (3) and Equation (4) to prevent overestimation of the shear strength of geological units in regional analysis. We also conducted a conventional Newmark analysis using assigned strengths, such as friction angle and cohesion, as shown in Table 1. The predicted displacements obtained through the conventional Newmark's method ranged from 0 to 78 cm, whereas the proposed method yielded values ranging from 0 to 76 cm. Subsequently, the calculation of the CF using the conventional Newmark analysis produced a range of -1 to 0.74, while the proposed method resulted in a range of -1 to 0.99. To evaluate the performance of the two methods, we utilized the area under the curve (AUC) analysis. The AUC plot illustrates the cumulative area of CFs within each interval of calculated values, representing the proportion of the total study area (x-axis), and the proportion of cumulative landslides falling within those CFs (y-axis) [59]. An AUC value of 0.5 suggests performance no better than random guessing, while a value of 1 indicates perfect performance [59]. Figure 14 depicts the results of the AUC analysis for both methods, with a calculated value of 0.57 for the proposed method and 0.53 for the conventional Newmark's method. Therefore, the introduced method yields improved results compared to the conventional Newmark analysis.

However, it is noteworthy that both methods yield relatively not large AUC values. This outcome can be attributed to the seismic landslide interpretation data employed in the analysis. As is well-known, landslides can be categorized into three distinct areas: the source area, the transport area, and the deposition area. The interpretation data utilized for the calculations encompasses all three areas. The initiation of landslides occurs in the source area, while the deposition area is where the landslide material comes to rest. The source area is typically characterized by steep slopes, leading to relatively large Newmark displacements that align with the model's predictions. Conversely, the deposition area tends to exhibit flatter terrain with smaller slopes, resulting in relatively smaller Newmark displacements. Nevertheless, following the occurrence of a landslide, a substantial amount of material accumulates in the deposition area. Consequently, a significant portion of the landslide area in the interpretation data falls within regions associated with smaller Newmark displacements. This leads to deviations in the predictive outcomes of the model. To achieve more accurate results, it is crucial to precisely delineate the source area of seismic landslides and employ it as the primary input for hazard assessment. This aspect warrants considerable attention in future assessments of seismic landslide hazards.

6. Conclusions

This study presents a pioneering method to assess the coseismic landslide hazards triggered by the 2013 Lushan earthquake, incorporating the roughness of potential sliding surfaces in the inner slope. The primary findings are outlined below:

(1) The process of mapping data from the 2013 Lushan earthquake demonstrates the viability of integrating Newmark analysis with Barton’s shear strength criterion. This approach holds practical utility in evaluating regional seismic hazards.

(2) By associating Newmark displacements with the certainty factor model, the effectiveness of Newmark’s method in predicting the hazard posed by coseismic landslides was enhanced.

(3) The outcomes of the AUC analysis demonstrate the higher reliability of the proposed method compared to the conventional Newmark approach.