An improved method of Newmark analysis for mapping hazards of coseismic landslides

Coseismic landslides can destroy buildings, dislocate roads, sever pipelines, and cause heavy casualties. It is thus important but challenging to accurately map the hazards posed by coseismic landslides. Newmark's method is widely applied to assess the permanent displacement along a potential slide surface and model the coseismic response of slopes. This paper proposes an improved Newmark analysis for mapping the hazards of coseismic landslides by considering the roughness and effect of the size of the potential slide surfaces. This method is verified by data from a case study on the 2014 M_w 6.1 (the United States Geological Survey) Ludian earthquake in Yunnan Province, China. Permanent displacements due to the earthquake ranged from 0 to 122 cm. The predicted displacements were compared with a comprehensive inventory of landslides triggered by the Ludian earthquake to map the spatial variation in the hazards of coseismic landslides using the certainty factor model. The confidence levels of coseismic landslides indicated by the certainty factors ranged from −1 to 0.95. A hazard map of the coseismic landslide was generated based on the spatial distribution of values of the certainty factor. A regression curve relating the predicted displacement and the certainty factor was drawn, and can be applied to predict the hazards of coseismic landslides for any seismic scenario of interest. The area under the curve was used to compare the improved and the conventional Newmark analyses, and revealed the improved performance of the former. This mapping procedure can be used to predict the hazards posed by coseismic landslides, and provide guidelines for decisions regarding the development of infrastructure and post-earthquake reconstruction.

1 Introduction

Earthquakes are recognized as one of the major causes of landslides (Keefer, 1984). Hazards caused by coseismic landslides have drawn increasing attention in recent years (e.g., Jibson et al., 1998, 2000; Khazai and Sitar, 2003; Qi et al., 2010, 2011, 2012; Chen et al., 2012; Xu et al., 2013; Yuan et al., 2014). The damage caused by seismically triggered landslides is sometimes more severe than the direct damage caused by the earthquake (Keefer, 1984). Estimating where a specific shaking is likely to induce a slope failure plays an important role in the regional assessment of coseismic landslides.

Pseudo-static analysis formalized by Terzaghi (1950), and finite-element modeling applied by Clough and Chopra (1966) have been employed to assess the seismic stability of slopes in early efforts (Jibson, 2011). Newmark (1965) first introduced a relatively simple and practical method, which is still commonly used nowadays to estimate the coseismic permanent displacements of slopes (Jibson, 2011). Studies have shown that Newmark's method yields reasonable and practical results when modeling the dynamic performance of natural slopes (Wilson and Keefer, 1983; Wieczorek et al., 1985; Jibson et al., 1998, 2000; Pradel et al., 2005). Rathje and Antonakos (2011) recently presented a unified framework for predicting coseismic permanent sliding displacement based on Newmark's method. Chen et al. (2018) used Newmark's method to calculate the minimum accelerations required for coseismic landslides in the region affected by the 2014 Ludian earthquake. Chen et al. (2019) subsequently developed an easy operation mapping method to assess hazards posed by coseismic landslides in the zone struck by the 2014 Ludian earthquake using Newmark's method.

Such applications generally start from an analysis of the dynamic stability of slopes, which is quantified as the critical acceleration. Barton model (Barton, 1973) has been widely used in rock mechanics and engineering to predict the shear strength of rock joints, which plays a crucial role in the calculation of critical acceleration. However, researchers have not adequately attended to the shear strength of rock joints during the assessment of coseismic landslides. To better estimate the dynamic stability of slopes, in this paper, we introduce the Barton model (Barton, 1973) to Newmark analysis to develop an improved modeling method for mapping the hazards posed by coseismic landslides using data from the 2014 Ludian earthquake in Yunnan Province in southwestern China. As predictions of coseismic landslides are not based on exact results, i.e., the computed permanent displacements, but are also mingled with unformalized expertise, i.e., the interpreted landslides, we present a model of inexact reasoning, i.e., the certainty factor model (CFM), that defies analysis, as an application of sets of inference rules that are expressed in predicate logic (Shortliffe and Buchanan, 1975) to produce a map of the hazards posed by coseismic landslides.

This paper briefly introduces the characteristics and spatial distribution of landslides triggered at the chosen site, describes the method of modeling used for the analysis of the stability of seismic slopes, presents the mapping procedure of the confidence level of seismic slope failure, and finally discusses the results of the assessment of seismic hazard, as well as a comparison with the conventional Newmark analysis.

2 Study area

The epicenter of the 2014 M_w 6.1 (the United States Geological Survey) Ludian earthquake was located on the southeastern margin of the Tibetan Plateau. A rectangular area lying immediately around the epicenter containing dense concentrations of the induced landslides was chosen for study (Fig. 1). The elevation of the area ranged from 785 to 3085 m above sea level. Three rivers – the Niulanjiang River, Shaba River, and Longquan River – pass through the study area (Fig. 1). The topography ranges from flat in the river valleys to nearly vertical in the slopes on the banks of the rivers. According to Chen et al. (2015), Niulanjiang River flows from the southeast (SE) to the northwest (NW), and incises to a depth between 1200 and 3300 m, resulting in about 80 % of the slopes having angles greater than 40^∘ distributed along the banks. The predominant geological units of the study area have an age that varies from the Proterozoic to the Mesozoic, including dolomite, limestone, shale, sandstone, basalt, and slate (Fig. 2).

Figure 1 Map of the study area showing the inventoried landslides.

Figure 2 Geological map of the study area showing lithology and faults.

An inventory of 1416 landslides triggered by the 2014 Ludian earthquake (Fig. 1) was compiled by visual inspection through comparisons between pre-earthquake satellite images obtained from Google Earth (30 January 2014) and 0.2 m high-resolution post-earthquake aerial images (7 August 2014; data provided by the Digital Mountain and Remote Sensing Applications Center, Institute of Mountain Hazards and Environment, Chinese Academy of Sciences, and Beijing Anxiang Power Technology Co., ltd.). A majority of landslides triggered by the earthquake were shallow, flow-like landslides (shallower than 3 m), developing in particularly dense concentrations along steeply incised river valleys. The total area of these interpreted landslides was 7.01 km² within a study area of 705 km². A detailed study showed that 846 of the mapped landslides were greater than 1000 m² in area, occupying 6.74 km² and accounting for 96.1 % of the total area of landslides, of which 279 were greater in area than 5000 m², occupying 5.37 km² and accounting for 76.6 % of the total landslide area.

3 Methodology

3.1 Modeling method

In the context of the analysis of the dynamic stability of a slope, Newmark (1965) proposed a permanent displacement analysis that bridges the gap between simplistic pseudo-static analysis and sophisticated but generally impractical finite-element modeling (Jibson, 1993). Newmark's method simulates a landslide as a rigid plastic friction block with a known critical acceleration on an inclined plane (Fig. 3) and calculates the cumulative permanent displacement of the block, as it is subjected to an acceleration time history of an earthquake. Newmark (1965) showed that the dynamic stability of a slope is related to the critical acceleration of a potential landslide block and can be expressed as a simple function of the static factor of safety and the geometry of the landslide (Jibson et al., 1998, 2000):

$$\begin{array}{}\text{(1)}& a_{\mathrm{c}}=\left(F_{\mathrm{S}}-\mathrm{1}\right)g\sin \mathit{\alpha },\end{array}$$

where a_c is the critical acceleration in terms of g, i.e., the acceleration due to the Earth's gravity, F_S is the static factor of safety, and α is the angle from the horizontal at which the center of the slide block moves when displacement first occurs (Jibson et al., 1998, 2000). For a planar slip surface parallel to the slope, this angle generally approximates to the angle of the slope.

Figure 3 Conceptual sliding-block model of Newmark analysis.

Natural slopes often develop a group of shallow unloading joints (Fig. 4) parallel to the surface due to valley incisions (Gu, 1979; Hoek and Bray, 1981). Studies have shown that rock slopes behave as collapsing and sliding failures of shallow unloading joints under strong earthquakes, and 90 % of coseismic landslides are shallow falls and slides (Harp and Jibson, 1996; Khazai and Sitar, 2003; Dai et al., 2011; Tang et al., 2015). According to Qi et al. (2012), two typical kinds of landslides are triggered by earthquakes: (a) shallow, flow-like landslides with a depth of less than 3 m in general and (b) rockfalls thrown by the shaking caused by the earthquake that usually occur at the crest of the slope. For both types, unstable blocks of rock are often cut and activated along the rock joints. Therefore, the static factor of safety in terms of the critical acceleration in these conditions is related to the peak shear strength of the rock joints. For the purpose of regional analysis, we use a limit equilibrium model of an infinite slope (Fig. 3) by referring to the simplification of Newmark's method by of Jibson et al. (1998, 2000). The value of the static factor of safety against sliding given by the ratio of resistance to the driving forces is determined by conventional analysis without considering accelerations, expressed as follows:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{F}_{\mathrm{S}}& ={\displaystyle \frac{\mathrm{resisting}\mathrm{force}}{\mathrm{driving}\mathrm{force}}}={\displaystyle \frac{\mathit{\tau }L}{mg\sin \mathit{\alpha }}}={\displaystyle \frac{\mathit{\tau }L}{\mathit{\gamma }Lt\sin \mathit{\alpha }}}\\ & ={\displaystyle \frac{\mathit{\tau }}{\mathit{\gamma }t\sin \mathit{\alpha }}},\end{array}\end{array}$$

where τ is the peak shear strength of the rock joint, L is the length of the rock joint, m is the mass of the failure rock block, γ is the unit weight of the rock mass, and t is the thickness of the failure rock block.

Figure 4 A schematic diagram showing shadow unloading joints in the slope.

For a Newmark analysis, it is customary to describe the shear strength of rocks instead of rock joints in terms of Coulomb's constants, i.e., friction angle (φ) and cohesion (c). However, both are not only stress dependent but also scale dependent (Barton and Choubey, 1977). According to Barton (1973), a more satisfactory empirical relationship for predicting the peak shear strength of a joint can be written as follows:

$$\begin{array}{}\text{(3)}& \mathit{\tau }={\mathit{\sigma }}_{\mathrm{n}}\tan \left[\mathrm{JRC}{\log }_{\mathrm{10}}\left({\displaystyle \frac{\mathrm{JCS}}{{\mathit{\sigma }}_{\mathrm{n}}}}\right)+{\mathit{\varphi }}_{\mathrm{b}}\right],\end{array}$$

where σ_n is the effective normal stress, JRC is the joint roughness coefficient, JCS is the joint wall compressive strength, and ϕ_b is the basic friction angle, i.e., the angle of frictional sliding resistance between rock joints, which can be obtained from residual shear tests on natural joints (Barton, 1973).

The effective normal stress (σ_n) generated by gravity acting on the rock block is as follows:

$$\begin{array}{}\text{(4)}& {\mathit{\sigma }}_{\mathrm{n}}={\displaystyle \frac{mg\cos \mathit{\alpha }}{L}}={\displaystyle \frac{\mathit{\gamma }Lt\cos \mathit{\alpha }}{L}}=\mathit{\gamma }t\cos \mathit{\alpha }.\end{array}$$

Considering the impact of size on JRC and JCS, the formulations developed by Barton and Bandis (1982) are shown as follows:

$$\begin{array}{}\text{(5)}& {\displaystyle }{\mathrm{JRC}}_n={\mathrm{JRC}}_{\mathrm{0}}{\left({\displaystyle \frac{L_n}{L_{\mathrm{0}}}}\right)}^{-\mathrm{0.02}{\mathrm{JRC}}_{\mathrm{0}}},\text{(6)}& {\displaystyle }{\mathrm{JCS}}_n={\mathrm{JCS}}_{\mathrm{0}}{\left({\displaystyle \frac{L_n}{L_{\mathrm{0}}}}\right)}^{-\mathrm{0.03}{\mathrm{JRC}}_{\mathrm{0}}},\end{array}$$

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

Hence, the static factor of safety (F_S) of a slope can be written as follows:

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{F}_{\mathrm{S}}& ={\displaystyle \frac{\mathit{\tau }}{\mathit{\gamma }t\sin \mathit{\alpha }}}={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{n}}\tan \left[{\mathrm{JRC}}_n{\log }_{\mathrm{10}}\left(\frac{{\mathrm{JCS}}_n}{{\mathit{\sigma }}_{\mathrm{n}}}\right)+{\mathit{\varphi }}_{\mathrm{b}}\right]}{\mathit{\gamma }t\sin \mathit{\alpha }}},\\ & ={\displaystyle \frac{\mathit{\gamma }t\cos \mathit{\alpha }\tan \left[{\mathrm{JRC}}_n{\log }_{\mathrm{10}}\left(\frac{{\mathrm{JCS}}_n}{\mathit{\gamma }t\cos \mathit{\alpha }}\right)+{\mathit{\varphi }}_{\mathrm{b}}\right]}{\mathit{\gamma }t\sin \mathit{\alpha }}},\\ & ={\displaystyle \frac{\tan \left[{\mathrm{JRC}}_n{\log }_{\mathrm{10}}\left(\frac{{\mathrm{JCS}}_n}{\mathit{\gamma }t\cos \mathit{\alpha }}\right)+{\mathit{\varphi }}_{\mathrm{b}}\right]}{\tan \mathit{\alpha }}}.\end{array}\end{array}$$

After calculating the angle of the slope and static factor of safety, the critical acceleration of the slope can be determined. Once the time history of the earthquake' acceleration has been selected, portions of the record lying above the critical acceleration a_c (Fig. 5a) are integrated once to derive a velocity profile (Fig. 5b); the time history of velocity is then integrated a second time to obtain the profile of cumulative displacement of the block (Fig. 5c). Users finally determine the dynamic performance of the slope based on the magnitude of the Newmark displacement (Jibson et al., 1998, 2000; Jibson, 2011). The detailed procedure for conducting a Newmark analysis with the Barton model is discussed in the following sections.

Figure 5 Demonstration of the Newmark analysis algorithm (adapted from Wilson and Keefer, 1983; Jibson et al., 1998, 2000).

3.2 Static factor-of-safety map

Considering that the mapped landslides greater than 1000 m² in area occupied 96.1 % of the total landslide area, we selected a 30 m × 30 m digital elevation model (DEM) from the ASTER Global Digital Elevation Model (https://doi.org/10.5067/ASTER/ASTGTM.002; NASA/METI/AIST/Japan Spacesystems and U.S./Japan ASTER Science Team, 2009), which facilitated the subsequent hazard analysis. A basic slope algorithm was applied to the DEM to produce a slope map (Fig. 6), where the slope was identified as the steepest downhill descent from a cell to its neighbors (Burrough and McDonell, 1998). The slopes ranged from greater than 60^∘ along the banks of the Niulanjiang River, Shaba River, and Longquan River to less than 20^∘ in low mountains and hills in the north and east.

Figure 6 Slope map derived from the DEM of the study area.

According to Jibson et al. (1998, 2000), slopes steeper than 60^∘ remain unstable even at high strengths. We assume that Newmark's rigid plastic block is unsuitable for such a steep sliding surface. In this case, sliding occurs along a plane at an angle (α) of $\mathrm{45}{}^{\circ }+{\mathit{\varphi }}_{\mathrm{b}}/\mathrm{2}$ with the horizon (Fig. 7). Therefore, we assigned an angle (α) of $\mathrm{45}{}^{\circ }+{\mathit{\varphi }}_{\mathrm{b}}/\mathrm{2}$ to slopes steeper than 60^∘ to avoid too small a value of F_S in the Newmark analysis.

Figure 7 Schematic map showing the angle (α) for slopes steeper than 60^∘. σ_1f and σ₃ are the major and minor principal stress in the state of limit equilibrium, respectively. ϕ_b is the basic friction angle.

Table 1 Shear strengths assigned to rock types in the study area.

Friction angle (φ), cohesion (c), and unit weight (γ) were derived from the Geological Engineering Handbook (Geological Engineering Handbook Editorial Committee, 2018).

Figure 8 JRC_n component of shear strength assigned to rock types in the study area.

The digital geological map from the China Geological Survey (CGS) was rasterized at a 30 m grid spacing to assign material properties throughout the study area. According to the literature, JRC₀ and JCS₀ depend strongly on lithology (Coulson, 1972; Barton and Choubey, 1977; Bandis et al., 1983; Bilgin and Pasamehmetoglu, 1990; Priest, 1993; Singh et al., 2012; Alejano et al., 2012, 2014; Giusepone and da Silva, 2014; Yong et al., 2018). Representative values of γ, JRC₀, JCS₀, and ϕ_b assigned to each rock type exposed in the study area were estimated using the test data listed in Table 1. The selected values were near the middle of the ranges represented in the references. These JRC₀ and JCS₀ values were considered in a laboratory scale for a length of 100 mm as L₀. For each grid cell in the regional analysis, the length of the engineering dimension, L_n, can generally be set as a 10-fold range of L₀. This is because the value of JRC_n∕JRC₀ (JCS_n∕JCS₀) is nearly constant when the value of L_n∕L₀ is greater than 10 (Bandis et al., 1981). The values of JRC_n and JCS_n were then calculated by inserting the values of JRC₀ and JCS₀, and L₀, and L_n into Eqs. (5) and (6), respectively. Figures 8 and 9 show the spatial distributions of JRC_n and JCS_n, respectively. The basic friction angle (ϕ_b) map and unit weight (γ) map are shown in Figs. 10 and 11, respectively.

Figure 9 JCS_n component of shear strength assigned to rock types in the study area.

Figure 10 Basic friction angle (ϕ_b) component of shear strength assigned to rock types in the study area.

Figure 11 Unit weight (γ) assigned to rock types in the study area.

For the sake of simplicity, the thickness of the modeled block t was taken to be 3 m, which reflected the typical slope failures of the Ludian earthquake. The static factor-of-safety map was produced by combining these data layers (α, JRC_n, JCS_n, ϕ_b, and γ) in Eq. (7). In the initial iteration of the calculation, grid cells in steep areas with static factors of safety smaller than 1 indicated that the slopes were statically unstable, but this did not necessarily mean that they were moving under shaking induced by the earthquake. In this condition, to avoid conservative results, we neither increased the strengths of the rock types with statically unstable cells nor adjusted the strengths of other rock types to preserve the differences in relative strength between them (as in Jibson et al., 1998, 2000). Instead, we assigned a minimal static factor of safety of 1.01, merely above limit equilibrium (Jibson et al., 1998, 2000), to these slopes to avoid a negative value of the critical acceleration a_c. According to Keefer (1984), most landslides triggered by earthquakes occur with a slope of at least 5^∘. The static factors of safety resulting from slopes of angles smaller than 5^∘ were very high. These slopes were unlikely to fail under the Ludian earthquake and did not produce a statistically significant sample in the analysis. Therefore, slopes less steep than 5^∘ were not analyzed in the second iteration. After the adjustment, the static factors of safety ranged from 1 to 17.4, as shown in Fig. 12.

Figure 12 Static factor-of-safety map of the study area.

3.3 Critical acceleration map

According to Newmark (1965), a pseudo-static analysis in terms of the static factor of safety and the slope angle was employed to calculate the critical acceleration of a potential landslide. The map of critical acceleration (Fig. 13) was generated by combining the static factor of safety and the slope angle in Eq. (1). The critical accelerations were derived from the intrinsic properties of the slope (topography and lithology), regardless of the given shaking. Therefore, the map of critical acceleration indicated the susceptibility of coseismic landslides (Jibson et al., 1998, 2000). The calculated critical accelerations ranged from nearly zero in areas that were more susceptible to coseismic landslides to greater than 1 g in areas that were less susceptible.

Figure 13 Map showing critical accelerations in the study area.

3.4 Shake map

There were 23 strong motion stations within 100 km of the epicenter of the Ludian earthquake (Fig. 14). Each station's record contained the three components of the peak ground acceleration (PGA), south–north direction, east–west direction, and up–down direction, as listed in Table 2 (the dataset was provided by the China Earthquake Data Center, http://data.earthquake.cn, last access: 16 June 2016). We calculated the average PGA of the two horizontal components of each strong motion recording and plotted a contour map (Fig. 15) using an inverse distance-weighted (IDW) interpolation algorithm. It determined the cell values using a linearly weighted combination of a set of sample stations with weights inversely proportional to distance (Watson and Philip, 1985). In addition, given that input stations far from the epicenter, where the prediction was made, might have had poor or no spatial correlation, we eliminated the input stations beyond 100 km from the epicenter from the calculation.

Figure 14 Locations of strong motion stations.

Table 2 Station records of three components of peak ground acceleration.

Figure 15 Contour map of peak ground acceleration (PGA) produced by the Ludian earthquake in the study area. PGA values shown are in g.

3.5 Newmark displacement map

In case of a landslide, in practice it is impossible to conduct a rigorous Newmark analysis when accelerometer records are unavailable. It is also impractical and time consuming to produce a displacement in each cell during the regional analysis. Therefore, empirical regressions (Ambraseys and Menu, 1988; Bray and Travasarou, 2007; Jibson, 2007; Saygili and Rathje, 2008; Rathje and Saygili, 2009; Hsieh and Lee, 2011) have been proposed to estimate Newmark displacement as a function of the critical acceleration and peak ground acceleration or the Arias intensity. Rathje and Saygili (2009) developed a vector model for displacement in terms of the critical acceleration (a_c), peak ground acceleration (PGA), and moment magnitude (M_w) based on an analysis of over 2000 strong motion recordings:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}\ln D& =\mathrm{4.89}-\mathrm{4.85}\left({\displaystyle \frac{a_{\mathrm{c}}}{\mathrm{PGA}}}\right)-\mathrm{19.64}{\left({\displaystyle \frac{a_{\mathrm{c}}}{\mathrm{PGA}}}\right)}^{\mathrm{2}}\\ & +\mathrm{42.49}{\left({\displaystyle \frac{a_{\mathrm{c}}}{\mathrm{PGA}}}\right)}^{\mathrm{3}}-\mathrm{29.06}{\left({\displaystyle \frac{a_{\mathrm{c}}}{\mathrm{PGA}}}\right)}^{\mathrm{4}}\\ & +\mathrm{0.72}\ln \left(\mathrm{PGA}\right)+\mathrm{0.89}\left(M_{\mathrm{w}}-\mathrm{6}\right),\end{array}\end{array}$$

where D is the predicted displacement in centimeters, and a_c and PGA are in units of g.

This model is a preferred displacement model at a site where acceleration time recordings are not available. Incorporating multiple parameters of ground motion into the analysis typically results in less variation in the prediction of displacement (Rathje and Saygili, 2009).

The Newmark displacement of each cell was calculated by combining the corresponding values of the critical acceleration, peak ground acceleration, and moment magnitude in Eq. (8). The predicted displacements ranged from 0 to 122 cm, as shown in Fig. 16.

Figure 16 Map showing predicted displacements in the study area.

3.6 Coseismic landslide hazard map

According to Jibson et al. (1998, 2000), predicted displacements provide an index of the seismic performance of slopes, where larger predicted displacements relate to a greater incidence of slope failures. But the displacements do not correspond directly to measurable slope movements in the field. To produce a coseismic landslide hazard map, we chose a model of inexact reasoning, the certainty factor model (CFM) created by Shortliffe and Buchanan (1975) and improved by Heckerman (1986), to explore the relationship between the occurrences of landslides and their predicted displacements. The CFM was created as a numerical method, initially used in MAYCIN, a backward-chaining expert system used in medical contexts (Shortliffe and Buchanan, 1975), for managing uncertainty in a rule-based system. In this model, the certainty factor (CF) represents the net confidence in a hypothesis H based on the evidence E (Heckerman, 1986). Certainty factors range between −1 and 1. A CF with a value of −1 means a total lack of confidence, whereas a CF with a value of 1 means total confidence. Values greater than zero favor the hypothesis while those less than zero favor its negation. According to Heckerman (1986), the probabilistic interpretation of CF is as follows:

$$\begin{array}{}\text{(9)}& CF=\left\{\begin{array}{ll}\frac{p\left(H\mathrm{|}E\right)-p\left(H\right)}{p\left(H\mathrm{|}E\right)\left[\mathrm{1}-p\left(H\right)\right]},& p\left(H\mathrm{|}E\right)>p\left(H\right)\\ \frac{p\left(H\mathrm{|}E\right)-p\left(H\right)}{p\left(H\right)\left[\mathrm{1}-p\left(H\mathrm{|}E\right)\right]},& p\left(H\mathrm{|}E\right)<p\left(H\right)\end{array}\right.,\end{array}$$

where CF is the certainty factor, p(H|E) denotes the conditional probability for a posterior hypothesis that relies on evidence, the posterior probability, and p(H) is the prior probability before any evidence is known. In the displacement analysis, p(H|E) was defined as the proportion of the area of the landslide within a specific displacement area and p(H) was defined as the proportion of the landslide area within the entire study area, excluding slopes less steep than 5^∘. In this way, the values of CF represented the confidence level for coseismic landslides. Positive values corresponded to an increase in the confidence level in slope failure while negative quantities corresponded to a decrease in this confidence. Higher positive values indicated higher confidence levels for coseismic landslides.

Given the above definition, we produced a coseismic landslide hazard map in terms of the certainty factors. First, displacement cells every 1 cm were grouped into bins such that all cells with displacements between 0 and 1 cm were grouped into the first bin, those with displacements between 1 and 2 cm were grouped into the second bin, and so on. The displacements were grouped into 123 bins, from 0 to 122 cm. We then calculated the proportion of cells occupied by areas of landslides in each bin. This proportion was considered the posterior probability of each bin as previously defined. The prior probability calculated by dividing the entire landslide area by the entire study area was the same in each bin. Finally, the values of CF were computed in each bin by using Eq. (9) to combine the corresponding values of the posterior and prior probabilities. The certainty factors ranged from −1 to 0.95. The values of CF indicated the confidence level of the occurrence of a landslide for each bin in the study area and provided the basis for producing the coseismic landslide hazard map.

As shown in the hazard map (Fig. 17), 73.2 % of landslides triggered by the Ludian earthquake were in areas with higher confidence levels with CF values greater than 0.6. The interpreted landslides were covered on the map to demonstrate their goodness of fit for the predicted confidence levels for coseismic landslides (Fig. 17).

Figure 17 Map showing confidence levels of coseismic landslides during the Ludian earthquake using the proposed method. Confidence levels are portrayed in terms of values of CF.

4 Results and discussion

The predicted displacements represent the cumulative sliding displacements for a given time history of acceleration. Based on the statistically significant sizes of the areas, displacements less than 60 cm, which was around the middle of the range of displacement, occupied about 80 % of the study area, while displacements greater than 80 cm occupied a very small area. Jibson et al. (1998, 2000) assumed that shallow falls and slides in brittle, weakly cemented materials fail following a relatively small displacement, whereas slumps and block slides in more compliant materials likely fail following a larger displacement. That is to say, the study area was more susceptible to rockfalls and shallow, disrupted slides that fail following a relatively small displacement. By contrast, it was subjected at a lower probability to coherent, deep-seated slides that would fail following a larger displacement. Indeed, the majority of landslides triggered by the Ludian earthquake were shallow, disrupted slides and rockfalls (Zhou et al., 2016). Although a few catastrophic rock avalanches, such as the Hongshiyan landslide (Chang et al., 2017), occurred in the field, they did not produce statistically significant samples that could meaningfully contribute to the model, which is consistent with the statistical results as discussed above. Therefore, the model should relate well to typical kinds of earthquake-induced landslides in the study area, thus demonstrating its usefulness in predicting the probability of other types of landslides.

According to Jibson et al. (1998, 2000), a function of CF and Newmark displacement would make it possible to predict the spatial variation in coseismic landslides in any scenario of interest involving the ground shaking. As mentioned above, 80 % of the study area featured predicted displacements of less than 60 cm. The numbers of the Newmark displacement cells were uneven. There were more cells in 1 cm bins for smaller displacements and fewer cells in 1 cm bins for larger ones. This might have affected the statistical significance of the function of CF and Newmark displacement. Therefore, the predicted displacement cells were grouped into bins based on quantile statistics. The breakpoints were 0, 10, 30, 39, 46, 51, 55, 59, 63, and 122. In this way, the number of cells in each bin was equal. Figure 18 shows, for each bin, the CF value of the Newmark displacement as plotted as a dot. As CF values ranged from −1 to 1, and not from 0 to 1, the Weibull (1939) curve developed by Jaeger and Cook (1969) is unsuitable here. Therefore, we modified the functional form as follows:

$$\begin{array}{}\text{(10)}& \mathrm{CF}=\mathrm{2}k\left[\mathrm{1}-\exp \left(-a{D}^b\right)\right]-\mathrm{1},\end{array}$$

where CF is the certainty factor, k is the maximum CF value represented by the data, D is the predicted displacement, and a and b are the regression constants. The regression curve based on data from the Ludian earthquake is

$$\begin{array}{}\text{(11)}& \mathrm{CF}=\mathrm{1.837}\left[\mathrm{1}-\exp \left(-\mathrm{0.073}{D}^{\mathrm{0.821}}\right)\right]-\mathrm{1}.\end{array}$$

From the curve shown in Fig. 18, when the predicted displacement increased, the value of CF increased monotonically, meaning that the confidence level for slope failure grew and landslide would probably occur. Such a procedure is consistent with the interpretation of certainty factor theory. Therefore, we were able to obtain estimates of the hazard different from the one used in this study using the same procedure described here.

Figure 18 CF as a function of Newmark displacement. A dot shows the CF value of a Newmark displacement bin; the red line is the fitting curve of the data using a modified Weibull function.

When fitting the results of shear tests using Coulomb's linear relation, the shear strengths varied widely from high normal stress in the laboratory to low normal stress in the field (Barton, 1973). We introduced the Barton model to the Newmark analysis to reduce the variation in shear strength in terms of Coulomb's constants. We also considered the impact of scale effects by using Eqs. (5) and (6) to prevent Newmark's method from underestimating the shear strength of geological units in regional analysis. In addition, for the Barton model, the joint roughness coefficient (JRC) was estimated from tilt tests, or by matching Barton's joint standard roughness profiles as defined by the International Society for Rock Mechanics (ISRM, 1978). The joint wall compressive strength (JCS) was estimated by Schmidt hammer index tests. These tests helped make a quick estimate of the shear strength in situ, which can facilitate the use of Newmark's method in an emergency hazard and risk assessment after an earthquake.

It is difficult for a statically stable slope to fail under an earthquake. Earthquakes usually cause slopes to fail in the state of limit equilibrium. For this reason, it is important to characterize the shear strength of the slope accurately. The shear strengths were assigned to the geological units using the results of hundreds of shear tests reported in the references provided in Table 1. We assigned the original shear strengths to the geological units, instead of increasing them to render the cells statically stable, as was done by Jibson et al. (1998, 2000). This would have changed the statically stable level of the entire study area, especially the slopes in the state of limit equilibrium. In addition, we considered the size effect of the potential slide surface, which could yield a lower F_S and, in turn, a higher displacement. However, the inventory of landslides was used to calibrate the predicted displacements, and the confidence levels indicated by the certainty factors fitted well with the spatial distribution of coseismic landslides, as shown in the hazard map (Fig. 17).

We also ran a conventional Newmark analysis using the assigned strengths, such as friction angle (φ) and cohesion (c), as shown in Table 1. The predicted displacements calculated by the conventional Newmark analysis ranged from 0 to 121 cm, compared with 0 to 122 cm as obtained by the proposed method. Figure 19 shows the hazard map produced using conventional Newmark analysis. The CFs ranged from −1 to 0.94, indicating a very similar result to that of the proposed method above. However, there were large differences along the Shaba River and upstream of the Niulanjiang River between the methods. By comparing Fig. 17 with Fig. 19, we see that the confidence levels of the proposed method fitted the data better than those of the conventional method, especially near upstream portions of the Niulanjiang River. The area under the curve (AUC) was employed to compare the performance of the methods. To create an AUC plot, the cumulative area of CFs within each interval of the calculated values, from the maximum to the minimum, was determined as a proportion of the total study area (x axis) and plotted against the proportion of cumulative landslides falling within those CFs (y axis) (Miles and Keefer, 2009). A value of 0.5 of the AUC indicates that performance is not better than a random guess and that of 1 indicates perfect performance (Miles and Keefer, 2009). Figure 20 shows the results of the AUC analysis of both methods. The calculated value for the proposed method was 0.58 while that for the conventional Newmark's method was 0.53. That is to say, the method introduced here yielded better results and is an improvement over the conventional Newmark analysis.

Figure 19 Map showing confidence levels of coseismic landslides during the Ludian earthquake using a conventional Newmark analysis. Confidence levels are portrayed in terms of values of CF.

Figure 20 Plots of area under the curve comparing the proposed method with the conventional Newmark's method.

5 Conclusion

Newmark's method is a useful physical model to estimate the seismic stability of natural slopes. The mapping procedure for data on the 2014 Ludian earthquake shows the feasibility of a Newmark analysis combined with Barton's shear strength criterion. Such a method has practical applications in the assessment of regional seismic hazard. We also considered the size effect of parameters of shear strength, such as the joint roughness coefficient (JRC) and the joint wall compressive strength (JCS), in regional analyses. Moreover, linking the Newmark displacements to the certainty factor model improved the utility of Newmark's method to predict the hazard posed by coseismic landslides. Finally, the results of an AUC analysis indicate that the proposed method is more reliable than the conventional Newmark method.