Distribution Pattern of Coseismic Landslides Triggered by the 2017 Jiuzhaigou Ms 7.0 Earthquake of China: Control of Seismic Landslide Susceptibility

Abstract

On 8 August 2017 an earthquake (M_S7.0) occurred within Jiuzhaigou County, Northern Aba Prefecture, Sichuan Province, China, triggering 4834 landslides with an individual area greater than 7.8 m² over a more than 400 km² region. Instead of correlating geological and topographic factors with the coseismic landslide distribution pattern, this study has attempted to reveal the control from seismic landslide susceptibility mapping, which relies on the calculation of critical acceleration values using a simplified Newmark block model. We calculated the average critical acceleration for each cell of the gridded study area (1 km×1 km), which represented the seismic landslide susceptibility of the cell. An index of the potential landslide area generation rate was defined, i.e., the possible landsliding area in each grid cell. In combination with PGA (peak ground acceleration) distribution, we calculated such indexes for each cell to predict the possible landslide hazard under seismic ground shaking. Results show that seismic landslide susceptibility plays an important role in determining the coseismic landslide pattern. The places with high seismic landslide susceptibility tends to host many landslides. Additionally, the areas with high potential landslide area generation rates have high real landslide occurrence rates, consistent with dominant small-medium scale landslides by this earthquake. This approach can aid assessment of seismic landslide hazards at a preliminary stage. Additionally, it forms a foundation for further research, such as the rapid evaluation of post-earthquake landslides and identifying highly impacted areas to help decision makers prioritize disaster relief efforts.

1. Introduction

Earthquake-triggered landslides in mountainous areas are one of the most common geologic hazards. Large-scale landslides can not only cause serious casualties and damage, but also block rivers and form barrier lakes in some cases, which threaten the safety of downstream areas. For example, the 2008 Wenchuan M7.9 earthquake created up to 828 barrier lakes [1], their breaches further worsening the consequences of landslides. Even small-scale landslides may block roads and affect quick rescue processes after the earthquake. For example, one of the significant disasters of landslides induced by the 2013 Lushan earthquake was the blockage to traffic. Since the 2008 Ms 8.0 Wenchuan earthquake, several strong earthquakes have occurred along the eastern and southeastern margins of the Tibetan Plateau, including the 2013 Ms 7.0 Lushan earthquake, the 2014 Ms 6.5 Ludian earthquake, and the 2017 Ms 7.0 Jiuzhaigou earthquake [2,3,4,5,6,7,8,9,10,11]. These earthquakes occurred on lesser-known active faults in the eastern and southeastern Tibetan Plateau, demonstrating the likely occurrence of future high magnitude earthquakes in the region. If we can determine the location and severity of the landslides rapidly after a major earthquake, it will undoubtedly provide very useful help for the rescue arrangement after the event. To this end, what controls coseismic landslide distribution is an important issue to address.

At present, while there are various statistical methods addressing the relationship between controlling factors and earthquake-triggered landslide distribution, many endeavors have been focused on pursuing the causes of the landslide distribution pattern. Compared with events-based methods, physical-based models are more objective and scientific due to the principles they follow. Among these models, the Newmark model is a physics-based model which has been widely and successfully applied to seismic landslide hazard assessment [12,13]. However, when Newmark’s method is used, people are generally interested in the ultimate prediction results, e.g., the permanent displacement, which is an index of landslide occurrence, while neglecting the middle result, e.g., critical acceleration, which is an index that portrays slope stability under seismic loading, namely seismic landslide susceptibility [14]. Unfortunately, the equations used for the permanent displacement predictions are generally calibrated using data from specific regions and applying them to other regions with different geological or climatic conditions will increase the uncertainty of the results. Contrary to the permanent displacement, a seismic landslide susceptibility map can provide useful information about where slopes are likely to fail, and how the intensity of shaking controls whether or not a landslide occurs. This knowledge is more important for assessing seismic landslide hazards at a preliminary stage. Additionally, it forms a foundation for further research, such as the rapid evaluation of post-earthquake landslides and identifying highly impacted areas to help decision makers prioritize disaster relief efforts.

This work used the landslides triggered by the 2017 Ms7.0 Jiuzhaigou, China earthquake to refine the Newmark method and develop a system for rapidly assessing the potential likeliness for a given area. Instead of correlating geologic and topographic factors with the coseismic landslide distribution pattern [3,15,16], we related the observed landslide distribution to existing seismic landslide susceptibility models. Then, combined with the distribution of peak ground acceleration (PGA), a potential coseismic landslide hazard map was prepared. Additionally the limitations and influence factors of the model are discussed.

2. Landslides Triggered by the Ms7.0 Jiuzhaigou Earthquake

On 8 August 2017, a Ms 7.0 earthquake (33.20°N, 103.82°E) occurred within Jiuzhaigou County, where the fault system splays into a horsetail structure composed of the Huya fault (F2 in Figure 1), Tazang fault (F3 in Figure 1), and the Minjiang fault (F1 in Figure 1) [17]. Focal-mechanism solutions showed a steeply dipping left-lateral strike-slip fault within the horsetail structure, while there was no obvious ground surface rupture by the event that was found [17]. The lack of such evidence raises the question about which fault within the Kunlun fault zone produced the Jiuzhaigou earthquake. The southern branch of the Tazang fault or an extension of the Huya fault have been considered as the most probable causative faults [18]. Historically, the Jiuzhaigou area has produced more than 50 events with Ms ≥ 5 in the past century, some with Ms ≥ 7 [19], and the faults here have different levels of activity. The Tazang fault is located along the Northeastern margin of the Kunlun system. It trends NW-SE and has a left-lateral strike slip rate of ~3.0 mm/yr, and a recurrence interval of ~2300 years [20]. On the western side of the horsetail structure, the Minjiang fault trends N-S and slips at an average rate of 0.37~0.53 mm/yr [21]. The Huya fault is located within the center of the Kunlun horsetail structure and has been active since the Late Quaternary. Historical events within the southern section of the Huya fault include an Ms7.2 and 6.7 in 1976 [22].

After the event, using 0.5m-resolution Geoeye-1 post-seismic images (shot on 14 August, 2017) and pre-seismic Google Earth (GE) images, the team led by Dr. Xu Chong interpreted 4834 coseismic landslides by the quake (personal communication), which was the most exhaustive inventory reported so far (Figure 1). The planar areas of these landslides are several to hundred thousand m², with the smallest being 7.8 m² and maximum 236,338 m², covering 9.64 km² [16]. Post-earthquake field investigations show that shaking triggered slope failures along steep riverbanks near the epicenter and along steep road cuts within the Jiuzhaigou parkland. The size of the failures ranged from a few cubic meters rock fall to large rock avalanches estimated to be a million cubic meters. Most of the landslides are distributed along the northern extension of the Huya fault in the NW direction and as much as 10 km from the mapped fault zone (Figure 1), which is consistent with the characteristics of landslides caused by earthquakes on strike slip faults. Two examples are the 2002 Mw 7.9 Denali earthquake, Alaska, USA and the 2010 Ms 7.0 Yushu earthquake, China, in which the majority of landslides were located close to the fault and concentrated on both sides of the fault within ~10 km [23,24]. From distances to the possible causative fault, the distribution pattern of the coseismic landslides suggests the causative fault is of strike-slip, likely the northwestward extension of the Huya fault (black dot line in Figure 1), although this needs further confirmation.

3. Methods and Data

According to Newmark [25], landslide potential can be simply modeled as a rigid block on an inclined plane with a known critical acceleration (a_c) necessary to overcome the shear resistance at its base. The cumulative permanent displacement of the block is modeled relative to the base of the block as it is subject to the effect of an earthquake’s acceleration and is used to predict the slope behavior during a shaking event. When an earthquake occurs, if a slope experiences ground motion which surpasses its critical acceleration, then the slope may fail during shaking. Therefore, the critical acceleration can be used to describe the stability of a slope under seismic shaking, i.e., seismic landslide susceptibility [12]. If a map showing the distribution of slope critical accelerations in a region can be prepared, in combination with a map of real seismic ground acceleration, it is possible to rapidly predict a landslide hazard after a major earthquake.

3.1. Slope Critical Acceleration of the Study Area as a Basic Map

First, we supposed that the seismic stability of slopes can be evaluated in terms of cumulative permanent displacement rather than a traditional minimum factor of safety (F_S). The critical acceleration, which is a simple function between the static factor of safety (F_S) and slope geometry (Equation (1)), represents a measure of intrinsic slope properties that are independent of any ground-shaking scenario. It is a connection between static and dynamic slope stability analysis and portrays seismic landslide susceptibility [12,26].

$$a_{\mathrm{c}}=\left(F_{\mathrm{s}}-1\right)g\times sina$$

(1)

where F_S is the static factor of safety and can be expressed as:

$$F_S=\frac{c^{’}}{\gamma tsin\alpha }+\frac{tan\phi ’}{tan\alpha }-\frac{m{\gamma }_{w}tan\phi ’}{\gamma tan\alpha }$$

(2)

The variables in Equation (2) are related to the slope material characteristics and slope geometry. When Fs is greater than 1, the slope is stable. Otherwise the slope is unstable.

Table 1. Variables in Equations (1) and (2).

Variable	Description
ac	critical acceleration in terms of g
g	acceleration of Earth’s gravity
FS	static factor of safety
α	angle from the horizontal direction
$\phi ’$	effective frictional angle
c´	effective cohesion
γ	material unit weight
γw	water unit weight
t	slope-normal failure slab thickness
m	proportion of the slab thickness saturated

shows the variables in Equations (1) and (2) and their descriptions.

From Equation (2), primary influence factors of stability Fs are material properties, the angle and saturation of the slope. Of them, the third item is related with seasons, which is ignored for non-raining seasons. As the Jiuzhaigou event did not occur in a rainy season, so the proportion of the slab thickness saturated is not considered, i.e., m = 0 in Equation (2). Meanwhile, field investigations show that the cosesimic landslides by the Jiuzhaigou shock are mostly shallow slope failures, thus the slope-normal failure slab thickness is assumed to be t = 2.0m.

In general, Fs of slopes before earthquakes should be equal or greater than 1. However, due to the effect of input parameters in calculating Fs, especially the uncertainty of rock properties, Fs may be less than 1 in the actual calculation. Some researchers use Fs = 1.1 to invert possible rock parameters [27]. However, such a way is unrealistic in a regional calculation. The purpose of this paper was to analyze the relationship between seismic landslide susceptibility and distribution of coseismic landslides, which does no need to determine whether landslides occur or not. Therefore, this work still uses representative rock parameters [12] for the slopes. Even if Fs<1, leading a_c<0, it does not influence the comparison of the slopes’ relative magnitudes, i.e., the stability of the slope with a_c = –0.3, which is less than slope with a_c = –0.2.

Geologic units from 1:200,000-scale geologic maps [28] were used to estimate material properties for the study region, including exposures of the strata from the Devonian to Quaternary periods, while lacking Jurassic, Cretaceous to Tertiary sequences (Figure 2).

Table 2. Strata and lithology of the study area.

Stratum	Lithology
Quaternary (Q)	Alluvial deposits, intercalated with mud, sand
Triassic (T)	Green-grey metamorphic tuffaceous sandstone, siltstone
Permian (P)	Massive limestone intercalated with carbonate rocks, siltstone
Carboniferous (C)	Grey-white thick bedded limestone, dolomitic limestone
Devonian (D)	Quartz sandstone, siltstone intercalated with mudstone

shows the strata and lithology of the study area. Almost all bedrocks appearing in this region were weathered and deformed due to intense tectonic movements. Although precise material strength parameters play an important role in slope stability analysis, it is not practical and beyond the scope of this study to test the parameters across such a large region. Therefore, we used a set of representative shear-strength values for each of the geologic units. First, the rocks in this study region were classified into four types, which included a hard rock group (type I), a moderately hard rock group (type II), a soft rock group (type III) and a second soft rock group (type IV) (Figure 2). Then, rock shear-strength parameters were assigned to each type based on the “Standard for engineering classification of rock masses” [29] and some other relevant references [12,27,30,31].

Table 3. Rock parameters in the study area.

Rock Type	c´(Mpa)	$\mathit{\phi }’(°)$	γ(kN/m3)
I	0.022	35	27.0
II	0.018	30	25.0
III	0.012	25	23.0
IV	0.010	23	20.0

shows material strength values for the rock types in the study region.

For this study, the slope angle was estimated using a slope map derived from 30 m Shuttle Radar Topography Mission (SRTM) elevation data [32] (Figure 3). For the assessment of slope stability under regional seismic shaking of the study area, such resolution of the data can meet the requirement of the analysis based on the gridding method. Additionally, for post-seismic rapid assessment of landslide hazard, the data SRTM30 was easy to obtain.

3.2. Acquisition of PGA Distribution of the Seismic Event

In general, seismic acceleration is easy to acquire for earthquakes. Here, we generated a PGA map using data from 13 seismic monitoring stations that were within 150 km to the epicenter and recorded ground motions of the 2017 Ms7.0 Jiuzhaigou earthquake. The nearest seismic station was 35 km from the epicenter and outside of the areas affected by co-seismic landslides (Figure 4, labeled as JZB). The three-component instrument recorded EW, NS, and vertical PGAs of 129.5, 185.0, and 124.7cm/s/s, respectively. We created a map showing the distribution of the PGA by interpolating the point data between sensors. To make up for insufficient data at the earthquake epicenter, we assigned a PGA of 260 gal based on calculations from the USGS website [33].

3.3. Analysis of Potential Seismic Landslide Development Degree

Based on distributions of slope’s critical acceleration and PGA, we evaluated the coseismic landslide potential by analyzing the difference between the seismic landslide susceptibility of a particular slope and the peak ground acceleration at that site following the Ms7.0 Jiuzhaigou earthquake. The study area was gridded by 1 km × 1 km. For each grid cell, its seismic landslide susceptibility was defined as:

$$S=\raisebox{1ex}{${{\displaystyle \sum }}_{i=1}^{i=N}{a}_i$}\!\left/ \!\raisebox{-1ex}{$N$}\right.$$

(3)

where a_i is the critical acceleration of the slope at the ith point calculated by Equation (1), and N is the number of calculation points of the grid, which depends on the resolution of DEM. In this case, we used SRTM30 to calculate the slope degree, therefore there were around 1000 dots in a 1 km × 1 km cell. Larger S values mean lower seismic landslide susceptibility, i.e., only a larger external force can make the slope failure, otherwise the slope can keep stable during seismic shaking. In contrast, smaller S values indicate low stability of the slope, where landsliding can easily occur.

We used the following method to describe the development degree of seismic landslides: For each grid cell, its potential landslide generation rate is defined as:

$$P_{ls}=\frac{\sum A_{site}(a_c<PGA)}{A_{cell}}$$

(4)

where A_cell is the area of the cell, and A_site is the area with PGA greater than the critical acceleration. Larger P_ls means larger area affected by landslides and higher hazard.

4. Results

Using the topographic data and geologic maps, we calculated and mapped critical accelerations for slopes throughout the study area using the Newmark method (Figure 5). In this analysis, slopes less than 10° were considered stable, and were therefore excluded from the calculation [12].

Consequently, a critical acceleration map based on cells was created to display the distribution of seismic landslide susceptibility of the study area (Figure 5). On this map, smaller critical accelerations (high seismic landslide susceptibility represented by red shades) signify slopes that are more likely to fail under seismic load, while blue colors indicate relatively stable slopes. We observed a good spatial correlation between the modeled seismic landslide susceptibility and the actual landslides within the extent of coseismic landslides (blue dot line in Figure 5).

Combining the map of seismic landslide susceptibility and PGA distribution permits the prediction of a potential coseismic landslide hazard after a major earthquake (Figure 6). The result is not only related with seismic landslide susceptibility before the event, but also affected by the PGA during seismic shaking. In this work, there were 278 grid cells within the limit of the landslide affected area (within blue dot line in Figure 6). Among them, there were 223 cells where coseismic landslides occurred and was determined by this work, accounting for 80% of the affected area. Comparison shows that the predicted (Pls) and actual landslide occurrence rates were well consistent expressed by cell numbers (Figure 7).

5. Discussion

As a basic map, seismic landslide susceptibility is commonly used in the research of earthquake-triggered landslides. For instance, in the rapid assessment of post-seismic landslides, with PGA data available, it permits one to compare slope critical acceleration and PGA thus predicting the area with a high landslide hazard. Additionally, in combination with zoning maps of seismic shaking parameters, it can aid medium-and long-term assessment of a seismic landslide hazard under varied exceeding probabilities. Brabb (1984) gave the following definition “a landslide susceptibility map depicts areas likely to have landslides in the future correlating some of the principal factors that contribute to landsliding, such as steep slopes and weak geological units, with the past distribution of landslides” [34]. While this paper followed the definition of seismic landslide susceptibility by Jisbon et al, [12], which highlights the possibility of induced landslides under seismic shaking. It is not difficult to find that some geologic factors in the traditional definition on landslide susceptibility, such as slope angles and geologic units, are actually contained in the rock mechanical parameters of the Newmark method. Additionally, critical acceleration can serve as an index to portray a measure of intrinsic slope properties independent of any ground-shaking scenario [12].

In this work, we used the grid cell as a basic unit to analyze the control of seismic landslide susceptibility on the distribution of coseismic landslides. Instead of using the minimum critical acceleration value (e.g., the most unstable) to assess a grid cell, we calculated the average critical acceleration for all the points in a grid cell to describe the cell’s seismic landslide susceptibility. Compared with using the minimum critical acceleration value as the index, although the critical accelerations of some points may be underestimated in such an averaged calculation, the overall feature of the grid cell can be better depicted. While using the minimum critical acceleration value as an index, the hazard may be overestimated, which would be a conservative assessment. In addition, this work defines an assessment index of potential landslide development area rate, which shows the feature of the Jiuzhaigou earthquake that small scale coseismic landslides are dominant. For example, dense distributions in the northwest and southeast to the epicenter were consistent with the places with high susceptibility predicted by this work. Moreover, in practice, the expression way of gridding can display the results more intuitively and clearly.

Some studies have discussed the limitation of the Newmark model. For example, it works fairly well for shallow landslides, but not for deeply seated landslides. A major limitation of the methods is that the accuracy of the analysis is subject to the accuracy of the inputs [12,27,30,35]. In this paper, for convenience, we merged the strata with similar lithology and assigned them same mechanical parameters. Such a treatment can directly influence the accuracy of the slope critical acceleration. As it is known, the genesis and distribution pattern of seismic landslides are influenced by multiple factors, and just because of this, various combinations of these factors create various types of landslides in different geologic and geomorphic environments. For the rigid block model, rock mechanical parameters are important in characterizing slopes, which can express the influence of slope materials on landslide distribution. In this study, we used a representative shear-strength for each rock type, even though each type was composed of several different geologic units (Figure 2). When the spatial variation of shear strength of geologic units is not considered, the prediction accuracy of the model prediction is mainly controlled by slope angles. For example, in the Quaternary stratum of the study area in this work (lower left of Figure 5), because this area has same material parameters, seismic landslide susceptibility is therefore controlled by slope angles, and landslide distribution is consistent with slope angles. Besides, inappropriate rock mechanical parameters also result in incorrect prediction. In a small range of the southeastern study area (Figure 5), the calculated critical accelerations were relatively bigger, meaning lower seismic landslide susceptibility, but many coseismic landslides occurred there. The possible reason for this error may be attributed to larger mechanical parameters assigned to the rock of this place, which does not accord with the real situation, leading to errors of prediction.

To better analyze the relationship between rock mechanical parameters and Fs, this work calculated slope angles of each rock type for Fs = 1.0. Figure 8 shows the relationships between slope stability (Fs) and slope angles for four rock types used in this study. Apparently, for Fs=1, under the rock properties defined (Table 3), different rock types corresponded to different slope angles, which were 23 º, 30 º, 38 º, 45 º (IV to I, i.e., soft to hard), respectively. This means that in a stratum, if the slope angle is greater than that of the stability, instability will occur. However, in real cases, many slope angles exceeding the stability angles exist in varied strata while they remain stable. In this work, considering the landslides are dominantly shallow slope failures, related to weathering material, we assigned a relatively small c’ value taking reference to some previous studies (e.g., Jibson et al., 2000), and obtained the slope stable angles in Figure 8. Therefore, how to assign reasonable rock mechanical parameters to slopes needs much further research, which is a key step for the methods based on a physical model to improve the prediction accuracy.

Dreyfus discussed the effect of parameter variation on prediction results in a quantitative manner, suggesting the influence of model parameters on the results was greater than those of the displacement model itself and ground motion parameter [30]. From the case in this work, the slope stability under seismic shaking, especially when slope critical accelerations value are similar, can pose a major effect on assessment results (see Equation (4)). Within 150 km to the epicenter of the Jiuzhaigou earthquake, there were 13 seismic stations deployed non-uniformly, most of which were far from the epicenter. Even if the station closest was 35 km apart, it was outside the limit of the landslides (Figure 4). Therefore, the data from these stations are actually not useful to the landslide quantitative analysis. We can only rely on mathematical interpolation to generate PGA distribution, which attenuates towards outside and around the epicenter without directional differences, not even considering the direction of PGA attenuation with respect to the strike-slip fault [36]. Additionally, the effect of terrain amplification was not considered either, which had obvious effects in some studies [37,38]. Therefore, the hazard assessment results presented here are at a preliminary stage. Inspiringly, more and more studies raise the possibility of improving the accuracy of landslide prediction in a regional scale, such as the integrated Spectral Element Method (SEM)- Newmark model applied in the Hong Kong island and the 2014 M6.5 Ludian earthquake in China [39]. Moreover, estimation of 3D topographic amplification point by point would make up the insufficiency of seismic data and refine the analysis in a future study.

As a triggering factor and a measurement of the magnitude of seismic ground motion, PGA is closely correlated with landslide occurrence. Jibson and Harp estimated the PGA between 0.02–0.04 g as the landslide distance limit for 23 August 2011 Mineral, Virginia, earthquake (Mw 5.8) [40]. In this case, can the real PGA on the boundary of the landslide distribution limit line be represented by the slope critical acceleration on the boundary? Although further study is needed to confirm this speculation, the critical acceleration possibly permits the use of the landslide distribution to obtain PGA distribution in the future work.

6. Conclusions

The slope critical acceleration from Newmark model analysis can well depict slope stability under seismic shaking. The feature that small- and medium-scale landslides are dominant of the Jiuzhaigou earthquake allows the successful use of this model to conduct analysis of seismic landslide susceptibility. This work makes such an analysis using a 1 km × 1 km grid of the study area. Results showed that the distribution of the landslides triggered by the 2017 Jiuzhaigou Ms7.0 earthquake was closely related to seismic landslide susceptibility expressed by the slope critical acceleration in the affected area. Overall, the distribution of the seismic landslide susceptibility values controls the pattern of coseismic landslides, i.e., the places with higher seismic landslide susceptibility have more slope failures. This work defines an index of the potential landslide area development rate to predict the post-seismic landslide hazard. In the study area, landslides are densely distributed in the northwest and southeast of the epicenter, largely consistent with the prediction using this index.

The simplification of slope material parameters input into the Newmark model can influence the accuracy of the calculation of seismic landslide susceptibility. Although relative magnitudes of rock mechanical parameters can be used in the early stage of assessment of regional landslide development degree, more reasonable methods for assigning parameters are required when higher accuracy should be attained. For the Jiuzhaigou Ms7.0 earthquake, because of the limited seismic station data, the accuracy of PGA distribution can also pose some effect on the calculation of the potential landslide area development rate.