Enhancing post-seismic landslide susceptibility modeling in China through a time-variant approach: a spatio-temporal analysis

ABSTRACT

The pre-phase landslides will have a legacy effect on future landslides, changing the landslide susceptibility. This study attempts to establish a reasonable post-seismic landslide susceptibility model and analyze the spatio-temporal characteristics of landslide susceptibility in the Jiuzhaigou M_S7.0 earthquake-struck region. Firstly, an integrated ‘space-ground’ monitoring technology is used to establish a multi-temporal post-seismic landslide dataset. Then, the buffer analysis method documents the spatio-temporal characteristics of post-seismic landslides. Thirdly, the distance is selected as an indicator to quantify the legacy effect. An improved time-variant model is established to evaluate the post-seismic landslide susceptibility. Finally, the spatio-temporal characteristics of landslide susceptibility are generalized, to sum up the changing law. Our results show that the post-seismic landslides are gradually closer to the pre-phase landslides with time. Distance is a critical factor in measuring the impact of pre-phase landslides on future landslides, which can improve the assessment accuracy of the post-seismic landslide susceptibility model. After a large seismic event, the correlation between landslide susceptibility and earthquakes gradually weakens. Post-seismic landslide prevention should focus on the pre-phase landslide expansion triggered by rainfall. Moreover, it should clean up the landslide deposits in time and reasonably dredge the debris flows to avoid secondary geological disasters.

KEYWORDS:

• Multi-temporal landslide dataset
• post-seismic time-variant landslide susceptibility model
• spatio-temporal characteristics
• legacy effect
• distance

1. Introduction

Strong earthquakes occurring in the mountainous region can not only trigger a large number of landslides but also cause severe vegetation destruction, surface deformation and cracks, geomorphic changes, and a large amount of loose material (Lin et al. 2009; Wang et al. 2014; Xu et al. 2016; Guo, Guo, and Feng 2021a; Carrión-Mero et al. 2021). This study uses the expression coseismic landslides to refer to landslides triggered within seconds to minutes of strong ground shaking and post-seismic landslides to refer to the landslide activities that occurred after a significant seismic event (Yegian, Marciano, and Ghahraman 1991; Wang 2004; Kritikos, Robinson, and Davies 2015; Adolfo et al. 2019; Fan et al. 2019a; Serey et al. 2019; Guo et al. 2021b). Previous studies show that the post-seismic landslides might be active over ten years, and then they will gradually become stable in strong earthquake-struck regions (e.g. Tang 2010; Cui et al. 2011; Huang and Li 2014; Chen et al. 2020). Therefore, evaluating the landslide susceptibility and summarizing the variation characteristics of post-seismic landslide susceptibility are critical for post-seismic landslide treatment and ecological restoration.

Currently, the conventional estimation models of landslide susceptibility can be divided into qualitative and quantitative ones. The qualitative model mainly relies on expert experiences to estimate susceptibility, and the quantitative model relies on data acquisition and analysis (e.g. Cao et al. 2019; Guzzetti et al. 1999; Shao et al. 2018; Caccavale, Matano, and Sacchi 2017; Zang et al. 2020; Adolfo 2021a; Arroyo-Solórzano, Quesada-Román, and Barrantes-Castillo 2022). The statistical models belong to quantitative models. Statistical models involve the statistical determination of combinations of variables that have led to landslide occurrence in the past. Then, they can predict the susceptibility of areas currently free of landslides but where similar conditions exist (Dai et al. 2001). Statistical models often have high accuracy and have been successfully used in many studies (e.g. García-Rodríguez et al. 2008; Budimir, Atkinson, and Lewis 2015; Lombardo and Mai 2018; Reichenbach et al. 2018; Adolfo 2021b).

Traditionally, landslide susceptibility is assumed time-invariant based on the assumption that the likelihood of landslide occurrence is time-independent for a given region and trigger event (Samia et al. 2020). The temporal component is usually investigated separately from the susceptibility in the above models and some early case studies (e.g. Aleotti 2004; Guzzetti et al. 1999; Wu et al. 2015). Because of this, few examples exist where the spatial and temporal components of the landslide hazard are combined (e.g. Ghosh et al. 2012; Fan et al. 2020). However, dynamic variations of landslide susceptibility in time are very common, and they exist in not only rainfall-induced landslides (e.g. Wang et al. 2022; Lin et al. 2022; Steger et al. 2022) but also earthquake-induced landslides (Roberts, Jones, and Boulton 2021). It is well-described that climatic, tectonic, and anthropogenic drivers may vary with time after the earthquake (Jones et al. 2021). Especially the pre-phase coseismic landslides could influence the landscape evolution and future landslide occurrences, known as the legacy effect. Therefore post-seismic landslide susceptibility also can be varied with time.

Recent studies have proposed the concept of time-variant landslide susceptibility based on the multi-temporal landslide dataset. It should be noted that time-variant landslide susceptibility still meets the definition of susceptibility, which does not include triggering factors. Fan et al. (2020) calculated the post-seismic time-variant landslide susceptibility by updating the multi-temporal post-seismic landslide dataset inputted into the susceptibility model. However, this method ignores the legacy effect between landslides of different phases. Samia et al. (2017a) accounted for landslide legacy effects across time (Lombardo and Tanyas 2021). Subsequently, researchers further described and quantified this phenomenon (Samia et al. 2017b). Based on the above research, Samia et al. (2018) put the legacy effect into the landslide susceptibility model, and the results showed that compared with conventional landslide susceptibility models, models considering the landslide legacy effect could significantly improve the evaluation accuracy of landslide susceptibility (Samia et al. 2020). However, landslide susceptibility models based on path dependency theory generally need high-quality datasets to quantify the legacy effect of the pre-phase landslides. They need to use the accurate geometric attributes of each landslide to establish the detailed spatial association of multi-temporal post-seismic landslides. Establishing a multi-temporal landslide dataset in a large region requires many resources and time (Reichenbach et al. 2018), and there is insufficient data for this type of research in many areas. The existing post-seismic time-variant susceptibility models put forward high requirements for establishing landslide datasets and processing multi-temporal geometric relationships. Globally, few multi-temporal landslide datasets can be used for this method. From this point of view, this method is challenging to implement and is currently only carried out in several study cases with small areas (Samia et al. 2017b; Roberts, Jones, and Boulton 2021).

The 2017 M_S7.0 Jiuzhaigou earthquake provides an excellent example to explore a reasonable post-seismic landslide susceptibility model and analyze the spatio-temporal characteristics of susceptibility. The multi-temporal landslide dataset is established using an integrated ‘space-ground’ monitoring technology. This study quantifies the legacy effect and attempts to develop an improved post-seismic time-variant landslide susceptibility model. An improved post-seismic time-variant landslide susceptibility model is simple and can be applied to more earthquake cases with a larger research area globally. This study also tries to obtain the changing laws of post-seismic landslide susceptibility by analyzing their spatio-temporal characteristics of them. It has important implications for predicting post-seismic landslides, defining hazard prevention zones, and evaluating post-seismic recovery in similar geological environments.

2. Material and methods

2.1. Study area and material

The Jiuzhaigou natural heritage site is located in the northeastern Tibetan Plateau. The terrain is high in the South and low in the North. It is the topographic transition zone from the Tibetan Plateau to Sichuan Basin. The Devonian to Triassic sedimentary rocks are prominent in the study area (e.g. Dong et al. 2020; Zhang et al. 2020). The most representative lithology is a thick sequence of deep marine deposits, including limestone, flysch complex, and sandstone. Quaternary loose deposits are distributed in river valleys (Dai et al. 2017). Deep-incision gullies and high mountains with slopes greater than 30° characterize the landforms. The maximum peak is over 4700 m a.s.l., and the minimum elevation is about 1200 m a.s.l. (Fan et al. 2018; Guo et al. 2021c) (Figure 1).

Figure 1. Jiuzhaigou M_S7.0 Earthquake. (a) The map shows the location of Jiuzhaigou, (b) Tectonic setting of the 2017 Jiuzhaigou M_S7.0 Earthquake (modified by Guo et al. 2021b). MJF: Minjiang Fault, NHYF: North Huya Fault, TZF: Tazang Fault.

At 21:19:46 on 8 August 2017, an M_S7.0 earthquake occurred in Jiuzhaigou County, Aba Autonomous Prefecture, Sichuan Province, China, with an affected area of more than 4000 km² (VII degrees seismic intensity and above). Due to the high magnitude, shallow source, high intensity of the earthquake, and fragile geological environment (high slope gradient, fragile lithology), this earthquake triggered thousands of landslides (e.g. Wu et al. 2018; Dong et al. 2020). Due to the dense vegetation in the area, the damage caused by landslides is significant and easy to identify in optical remote sensing images. Thus, the Jiuzhaigou area can provide a rare natural laboratory for studying post-seismic landslide susceptibility. Moreover, as a United Nations Educational, Scientific and Cultural Organization (UNESCO) natural heritage site, monitoring and evaluating landslides in Jiuzhaigou are essential and urgent for tourism safety and ecological environment protection and management.

This study used the pre-seismic and post-seismic Sentinel-2 images (10 m), post-seismic aerial images (1 and 0.16 m), and a series of Google Earth pre-seismic images and post-seismic images (about 0.5∼1 m) to identify landslides (Table 1) and establish multi-temporal post-seismic landslide dataset. The remote sensing images with multi-spatial resolution can ensure the accuracy and completeness of the landslide dataset and avoid false extraction. This study selects Sentinel-2 images that covered the entire study area from a similar time each year (29 July 2017, 7 September 2017, 27 September 2019, and 27 August 2020) with low cloud and good data quality. The quality of remote sensing images in 2018 was poor, and the data were not used. The ortho-level aerial remote sensing image of emergency observation (10 August 2017) is provided by the Bureau of Surveying and Mapping. The post-seismic ortho-level unmanned aerial vehicle (UAV) image acquired by the UAV (CW-15) was taken on 10 June 2020. The ortho-level aerial remote sensing image of emergency observation and the UAV images cover all the key areas. The Google Earth images were taken on 13 August 2013, 21 October 2015, and 27 September 2019, covering the entire study area. These Google images with high spatial resolution were used due to without aerial images in 2013, 2017, and 2019. Moreover, all of the images are registered before using them. This processing method can also help ensure the completeness and consistency of the landslide dataset.

Table 1. Remote sensing data used in this study.

Data type	Data sources	Data acquisition time	Spatial resolution (m)
Pre-seismic data	Google Earth images	2013/08/13 2015/10/21
	Sentinel-2 images	2017/07/29	10
Coseismic data	Image of emergency observation	2017/08/10	1
	Sentinel-2 images	2017/09/07	10
Post-seismic data	Sentinel-2 images	2019/09/27 2020/08/27	10
	Google Earth images	2019/09/27	0.5-1
	CW-15 UAV image	2020/06/10	0.16

Additionally, this study selected seismic, topographic, geological, and human activity factors as the basic influencing factors of susceptibility, which had been widely used in previous studies (e.g. Huang and Li 2008; Keefer 2000). They included peak ground acceleration (PGA), slope, aspect, elevation, strata, and distances to seismogenic fault, river, and road (Figure 2) (Fan et al. 2020; Talaei 2014; Qi et al. 2011; Zou et al. 2022). Among them, terrain data such as slope, aspect, and elevation were obtained by ALOS DEM (12.5 m) (ASF DAAC n.d.) (http://search.asf.alaska.edu/, accessed May 2, 2022). Strata and seismogenic fault data were from the 1:200000 geological map (Li et al. 2017) (http://dcc.ngac.org.cn/geologicalData, accessed May 2, 2022). The USGS website (USGS Earthquake Hazards Program 2017) (http://earthquake.usgshakemap.gov, accessed May 2, 2022) provided seismic data such as PGA. The data on human activities, such as roads and rivers, came from the national data on roads and rivers (Geospatial Data Cloud 2019) (http://www.gscloud.cn/sources, accessed May 2, 2022).

Figure 2. Influencing factors in Jiuzhaigou (modified by Guo et al. 2021c). (a) Slope, (b) Aspect, (c) Elevation, (d) PGA, (e) Seismogenic fault (Dai et al. 2017; Fan et al. 2020), (f) Strata, (g) River, (h) Road.

2.2. Methods

The steps of this study are listed as follows (Figure 3): (1) The multi-temporal post-seismic landslide dataset is quickly and accurately established by the integrated ‘space-ground’ monitoring technologies, combined with automatic classification technology; (2) The spatio-temporal characteristics of multi-temporal post-seismic landslides are summarized and obtaining the typical factor quantifying the legacy impact of pre-phase landslides on future landslides; (3) The typical influencing factor is put into the susceptibility model, establishing an improved post-seismic time-variant landslide susceptibility model and quickly getting more accurate landslide susceptibility results; (4) The spatio-temporal characteristics of post-seismic landslide susceptibility are analyzed.

Figure 3. Procedural flow chart of this study.

2.2.1. Landslide identification

The post-seismic landslides are new landslides that occur in areas free of landslides in the past (Wang 2004). The post-seismic landslides mainly come from two types of landslides. One is the landslide after the earthquake formed by the expansion of the pre-phase landslide (pre-seismic landslide, coseismic landslide, or post-seismic landslide in the previous phase), which is referred to as the pre-phase landslide expansion. Another one is the single landslide after the earthquake, far away from the pre-phase landslide, referred to as the new single post-seismic landslide (Li, Wang, and Liu 2018; Guo et al. 2021b).

The Jiuzhaigou area has dense vegetation, with a coverage rate of 85.5% (Deng 2011). In most landslide areas of Jiuzhaigou, the vegetation on the slopes was destroyed. So, it is easy to identify the bare land where landslides damaged vegetation. The vegetation in the area above 3700 m a.s.l. is sparse due to frost weathering and permanent snow, and landslides are more challenging to find there. However, previous studies showed that landslides above 3700 m are rare in the Jiuzhaigou area (Wu et al. 2018). Therefore, this study only extracts landslides in areas lower than 3700 m a.s.l. to ensure a more accurate dataset.

In this study, we used an integrated ‘space-ground’ monitoring technology to extract landslides within 2055 km² of the Jiuzhaigou area. Firstly, the support vector machine (SVM) with high accuracy and efficiency (Vapnik 1998) was used to preliminarily identify the bare land (exposed ground without plants, snow, or water) which is larger than half a pixel (10 m/pixel) by Sentinel-2 remote sensing images (Table 1). Then, visual interpretation was used to verify whether the bare land represented the landslide results and eliminate false extraction areas by the aerial remote sensing images (aerial remote sensing image of emergency observation or UAV remote sensing images) or Google Earth images. The false extraction areas mainly include houses, roads, squares, and other impervious surfaces. The impervious surfaces such as houses, roads, and squares have significant geometric features, making them easy to identify and eliminate by visual interpretation. Finally, we conducted field investigations to verify the results for key or uncertain areas.

The multi-temporal landslide dataset is established quickly and accurately. The landslide areas of each phase in the dataset include all landslide areas in the previous phase. For example, the 2017 landslide areas include pre-seismic landslide areas and coseismic landslide areas; the 2019 landslide areas include pre-seismic landslide areas, coseismic landslide areas, and 2017–2019 post-seismic landslide areas (the landslides occurred from August 2017 to September 2019). After obtaining the accurate landslides polygons, this study used the 2017 landslide areas to erase the 2019 landslide areas to get the post-seismic landslide polygons in 2017-2019, and then this study used the 2019 landslide areas to erase the 2020 landslide areas to obtain the post-seismic landslide polygons in 2019-2020. Finally, a multi-temporal post-seismic landslide dataset was established. The sources of remote sensing data and the landslide extraction method in this study can help ensure the completeness and consistency of the landslide dataset.

2.2.2. Landslide spatio-temporal characteristics

Pre-phase landslides have legacy effects, impacting the new post-seismic landslides in the future (Roberts, Jones, and Boulton 2021; Samia et al. 2017a). However, the existing indicator quantifying the legacy effect is not simple enough. This study used buffer analysis to explore the spatio-temporal relationship between post-seismic landslides of each phase and their pre-phase landslides. Then, this study analyzed the spatio-temporal characteristics and summarized the spatio-temporal distribution law of post-seismic landslides. Finally, the distance factor is used to quantify the legacy effect in each phase. It helped to establish an improved post-seismic time-variant landslide susceptibility model.

The buffer analysis of the 2017–2019 post-seismic landslides distribution was taken as an example to explain the principle, as shown in Figure 4. This method starts from the edge of the pre-phase landslides and uses 50 m (or other length) as the buffer interval to generate a multi-level buffer zone. Then, this study counted landslides falling in different buffer zones. Notably, the distance factor refers to how long the landslide of each phase is from the edge of pre-phase landslides. For 2017–2019 post-seismic landslides, buffer zones started from the edge of the landslide in 2017; for 2019–2020 post-seismic landslides, buffer zones developed from the edge of the landslide in 2019.

Figure 4. Example of 2017–2019 post-seismic landslides distribution buffer analysis.

2.2.3. Landslide susceptibility model

Statistical models are frequently used in landslide susceptibility methods. The common statistical models can be divided into single models and coupling models. The common single models include the information (I) model, certainty factor (CF) model, SVM, logistic regression (LR) model, and so on. The common coupling models are coupled by some single models, such as the information and logistic regression (I + LR) coupling model (Guo et al. 2021c). The coupling models are usually much better than the single models which make up them (Zhang et al. 2020). The I + LR coupling model is composed of the I and the LR models and is easy to implement. It has high accuracy and is suitable for landslide susceptibilities in different geological backgrounds, such as the Alpine karst region and the loess region (e.g. Fan et al., 2018; Tian et al. 2020; Guo et al. 2021c). Therefore, this study selected the I + LR model as the primary model to establish an improved post-seismic time-variant landslide susceptibility model. The I model, LR model, and I + LR coupling model are as follows.

In the I model, the landslide probability is represented by the entropy decrease in the landslide prediction process. The average uncertainty reduction caused by the combination of the factors equals the entropy change of the landslide system. The greater the information value is, the greater the landslide possibility is (Yin 2013).

In every factor $x_i$ (i = 1, 2, 3, … , n), the information value of each category $j$ (j = 1, 2, 3, … , m) can be listed as follows (Ba et al. 2017; Pasang and Kubicek 2018; 2020): (1) $I_{ij}=\log {\displaystyle \frac{N_j/S_j}{N/S}}$(1) where $I_{ij}$ is the information value of the category $j$ in the factor $x_i$; $N_j$ is the landslide area in the category $j$; $N$ is the total area of landslides in the whole study area; $S_j$ is the area of category $j$; $S$ is the total area of the whole study area; $n$ is the number of factors considered; $m$ is the number of categories in every factor. Therefore, the total information I, representing landslide possibility, is the sum of the information values of all influencing factors $x_i$. (2) $I={\sum }_{i=1}^{\mathrm{n}}{I}_{ij}={\sum }_{i=1}^{\mathrm{n}}\log {\displaystyle \frac{N_j/S_j}{N/S}}$(2) where $n$ is the number of factors considered. If I < 0, the landslide possibility is lower than the average of the whole area; if I = 0, the landslide possibility is equal to the average of the whole area; if I > 0, the landslide possibility is higher than the average of the whole area. The greater the information value I is, the greater the landslide possibility.

The LR model predicts the probability of events by establishing a regression relationship between a dependent variable and multiple independent variables. In analysis, the dependent variable $Y$ is a binary variable representing whether there is a landslide. The value ‘1’ and ‘0’ of $Y$ indicate the occurrence and non-occurrence of landslides, respectively. An LR model has $n$ independent variables ($x_1$, $x_2,$ $x_3$, … , $x_n$). Influenced by $n$ factors, the conditional probability of landslide is $P$= $P$($Y$ = 1 | $x_1$,$x_2$,$x_3$, … , $x_n$). Then the model can be expressed as (Ayalew and Yamagishi 2005; Mathew, Jha, and Rawat 2009; Conoscenti et al. 2014): (3) $Z=a_0+a_1{x}_1+a_2{x}_2+\dots a_{\mathrm{n}}{x}_{\mathrm{n}}$(3) (4) $P={\displaystyle \frac{\exp (z)}{1+\exp (z)}}$(4) where Z is the linear combination; $x_i$ ($i$ =  1, 2, 3, … , n) is the independent factor; $a_0$ is the regression constant; $a_i$ is the coefficient that measures the contribution of the independent factor $x_i$ ($i$ =  1, 2, 3, … , n); $P$ is the prediction value of landslide probability.

In the conventional I + LR coupling model, the information value $I_{ij}$ in each factor calculated by the formula (1), substitute for the $x_i$ of the formula (3) as the index value of the LR model. We calculate the constant and regression coefficients of the LR model by Statistical Product Service Solutions (SPSS) 25 software (IBM Corp 2017) and then put them into formula (4) to obtain P, the landslide probability.

The multi-temporal post-seismic landslide dataset is updated quickly and accurately by automatic classification technology in this study. Although the post-seismic landslide susceptibility can vary with time by putting the updated multi-temporal landslide dataset into the conventional susceptibility model, the conventional susceptibility model still can not consider the legacy effect. So, the conventional susceptibility model can not accurately calculate the post-seismic time-variant susceptibility.

In this study, the distance factor came from the buffer analysis of the post-seismic landslide distribution in each phase. We put the distance factor, whose information value can represent the legacy effect of landslides, into the conventional landslide susceptibility model, establishing an improved time-variant landslide susceptibility model. The improved landslide susceptibility model not only quantifies the legacy effect but also uses the multi-temporal post-seismic landslide data. In the improved landslide susceptibility model, the information value $I(x_i,H)$ of each factor, including the distance factor in each phase calculated by the formula (1), substitute for the $x_i$ of the formula (3) as the index value of the LR model. The subsequent processing of the improved landslide susceptibility model is the same as that of the conventional I + LR coupling model. We also calculate the constant and regression coefficients of the LR model in each phase by using SPSS 25 software and then put them into formula (4) to obtain P, representing the landslide probability of each phase.

2.2.4. Accuracy verification of susceptibility results

The landslide susceptibility was divided into five classes using Jenks natural breaks optimization: very low, low, moderate, high, and very high. The Jenks natural breaks optimization divides the data into similar classes according to the natural (non-artificial) turning points or breakpoints of the data, making the classification very objective. Many studies have shown that this classification method is reasonable (e.g. Aditian, Kubota, and Shinohara 2018; Dai et al. 2017).

This study quantitatively evaluated the accuracy of the post-seismic time-variant landslide susceptibility model by the receiver operating characteristic (ROC) curves. The ROC curve is regarded as an effective method to evaluate the prediction ability of models (Swets 1988). The abscissa of the curve is the specificity or the false positive rate, which is the probability that the non-landslide is judged as a landslide by the model; the ordinate is the sensitivity or the true positive rate, which is the probability that the landslide is judged as a landslide. The area under the curve (AUC) is a quantitative index of the ROC curve to evaluate the model. The AUC value ranges from 0.5–1. The closer the AUC value is to 1, the better the prediction effect of the model is. If AUC = 0.5, the prediction result of the model is random; if AUC is at 0.5–0.7, the accuracy of the model is low; if AUC is at 0.7–0.9, the accuracy of the model is high; if AUC > 0.9, the accuracy of the model is very high (Fawcett 2006). In this study, the ROC curves and AUC values of each model were obtained by SPSS 25 software.

3. Results

3.1. Landslide inventory

The landslide areas of 2017 (after the Jiuzhaigou earthquake), 2019, and 2020 extracted in this study are 13.7, 14.9, and 16.5 km², respectively. In the landslide identification process, visual interpretation of Google Earth images, aerial remote sensing images, and field investigations were conducted to verify the landslide results, such as landslides near Five Flower Lake, Mirror Lake, and Heye community (Figures 5–7). It should be noted that all of the above landslide areas included the area of coseismic landslides, see Figure 8.

Figure 5. Remote sensing images and field photo around Five Flower Lake. (a) 2013, (b) 2017, (c) 2020, (d) Field photo.

Figure 6. Remote sensing images and field photo around Mirror Lake. (a) 2013, (b) 2017, (c) 2020, (d) Field photo.

Figure 7. Remote sensing images and field photo around Heye community. (a) 2013, (b) 2017, (c) 2020, (d) Field photo.

Figure 8. Landslides area distribution after the 2017 Jiuzhaigou M_S7.0 earthquake. (a) 2017, (b) 2019, (c) 2020.

Then, based on the above landslide areas of 2017 (after the Jiuzhaigou earthquake), 2019, and 2020, this study used the ArcGIS-Analysis-Erase tool to get the real post-seismic landslides. The results are as follows: the 2017–2019 post-seismic landslide area is 3.62 km²; the 2019–2020 post-seismic landslide area is 3.38 km² (Figure 9).

Figure 9. Post-seismic landslides. 2017–2019: (a) The Jiuzhaigou region, (b) The Jiuzhai Paradise region, (c) The Panda Lake region; 2019–2020: (d) The Jiuzhaigou region, (e) The Jiuzhai Paradise region, (f) The Panda Lake region.

3.2. Spatio-temporal characteristics of post-seismic landslide

This study generated buffer polygons with different distances in 2017–2019 and 2019–2020. The buffer polygons of each phase were used to cut all the post-seismic landslide polygons of each phase, not one polygon. Then, this study made statistics on the area of the already cut post-seismic landslide polygons falling in different buffer polygons. The 2017–2019 and 2019–2020 results are shown in Table 2 and Figure 10. Over 85% of the 2017–2019 post-seismic landslides are distributed within 200 m around the 2017 coseismic landslide area. More than 95% of the 2019–2020 post-seismic landslides are distributed within 200 m around the 2017–2019 landslide area. The trends in Figure 10 and our field investigation show that the pre-phase landslide expansion is the primary source of landslides after the Jiuzhaigou earthquake (Figure 11). Moreover, the post-seismic landslides in 2019–2020 are closer to the pre-phase landslides than those in 2017–2019 (Figure 10). More than 90% of the 2019–2020 post-seismic landslides are distributed within 50 m around the 2017–2019 landslide area. It suggests that the probability of new single landslides occurring far away from the pre-phase landslides is reduced with time. The spatio-temporal distribution of landslides shows more significant spatial aggregation characteristics. Therefore, the distance to the pre-phase landslide is a reasonable indicator to quantify the legacy effect of the pre-phase landslide.

Figure 10. Spatio-temporal distribution percentage of post-seismic landslides.

Figure 11. Examples of post-seismic landslides. Five-Flower Lake landslide edge: (a) 2017, (b) 2020; Heye landslide edge: (c) 2017, (d) 2020; Shuzheng landslide edge: (e) 2017, (f) 2020.

Table 2. Spatio-temporal distribution of post-seismic landslides in 2017–2019 and 2019–2020.

Buffer (m)	Post-seismic landslides within the buffer zone in 2017–2019		Post-seismic landslides within the buffer zone in 2019–2020
	Area (km^2)	Percentage of area (%)	Area (km^2)	Percentage of area (%)
50	2.59	71.49	3.15	93.17
100	2.88	79.59	3.25	96.27
150	3.01	83.15	3.28	96.90
200	3.09	85.30	3.29	97.32
250	3.14	86.74	3.30	97.65
300	3.18	87.88	3.31	98.06

3.3. Post-seismic landslide susceptibility

This study generates five training datasets and five testing datasets in 2017–2019 and 2019–2020, respectively, to compare the accuracy and robustness of the improved and the conventional landslide susceptibility model. The method of obtaining each training and testing dataset is as follows: 3000 points are randomly selected in the landslide area by ArcGIS. Then random numbers were assigned to the landslide points by Excel. The random numbers were sorted to select the first 2250 points as training samples and the remaining 750 landslide points as testing samples. The same procedure was repeated in the non-landslide area. Finally, this study can produce a training dataset with 4500 points (2250 landslide points and 2250 non-landslide points) and a testing dataset with 1500 points (750 landslide points and 750 non-landslide points). This method ensures that the points in the training and testing datasets do not overlap.

Five training datasets are used to train the improved landslide susceptibility model and the conventional landslide susceptibility model in each phase, obtaining five improved landslide susceptibility models and five conventional landslide susceptibility models. Then the testing dataset is used to evaluate the accuracy, and the ROC curves are obtained by SPSS 25 software (Figure 12). The results show that the accuracy of the five improved landslide susceptibility models is higher than that of the five conventional landslide susceptibility models in each phase. The improved landslide susceptibility model in this study is robust and highly accurate.

Figure 12. ROC curves and AUC values of landslide susceptibility models. 2017–2019: (a) Testing dataset 1, (c) Testing dataset 2, (e) Testing dataset 3, (g) Testing dataset 4, (i) Testing dataset 5; 2019–2020: (b) Testing dataset 1, (d) Testing dataset 2, (f) Testing dataset 3, (h) Testing dataset 4, (j) Testing dataset 5.

Although selecting the group of models with the best performance in this study can achieve the highest accuracy, other studies generally only have one training dataset and hope to obtain good results. Therefore, it is more realistic for us to randomly select a group instead of selecting the most effective group. Finally, this study randomly selects a group of improved and conventional landslide susceptibility models in each phase to calculate the post-seismic landslide susceptibility of the Jiuzhaigou area, as shown in Figure 13 and Figure 14.

Figure 13. Post-seismic landslide susceptibility calculated by the improved time-variant landslide susceptibility model. (a) 2017–2019, (b) 2019–2020.

Figure 14. Post-seismic landslide susceptibility calculated by conventional landslide susceptibility model. (a) 2017–2019, (b) 2019–2020.

For spatial distribution, the post-seismic landslide susceptibilities obtained by the conventional susceptibility model and the improved landslide susceptibility model have similar spatial distribution, see Figure 13 and Figure 14. But, the high and very high susceptibility areas of the latter are significantly reduced compared with the former. The susceptibility in Zhangzha town, which is the human habitation, is different too. The field investigation shows that the post-seismic landslide activities near Zhangzha town are weak, and no large-scale new landslides exist. Therefore, the spatial distribution of the improved landslide susceptibility model is more reasonable than the conventional susceptibility model.

This study used the coseismic landslide susceptibility in 2017 (Figure 15) (Guo et al. 2021c) as a benchmark for analyzing the changing trend of post-seismic landslide susceptibility. Then, the proportions of each susceptibility class in 2017, 2017–2019, and 2019–2020 landslide susceptibility are counted (Table 3 and Table 4).

Figure 15. Coseismic landslide susceptibility in 2017 (modified by Guo et al. 2021c).

Table 3. Changes of post-seismic landslide susceptibility distribution by the improved landslide susceptibility model.

Susceptibility	Percentage of classification area (%)
	2017	2017–2019	2019–2020
Very low	13.56	51.87	55.57
Low	22.65	23.83	24.71
Moderate	31.18	13.29	10.85
High	21.73	7.37	5.61
Very high	10.88	3.63	3.10

Table 4. Changes of post-seismic landslide susceptibility distribution by the conventional landslide susceptibility model.

Susceptibility	Percentage of classification area (%)
	2017	2017–2019	2019–2020
Very low	13.56	37.94	40.95
Low	22.65	25.87	26.51
Moderate	31.18	16.94	16.30
High	21.73	11.63	10.24
Very high	10.88	7.62	6.00

The moderate landslide susceptibility area in 2017 accounted for the largest proportion, 31.18%. The proportions of the very high and very low landslide susceptibility areas are close, 13.56% and 10.88%, respectively. The post-seismic landslide susceptibility calculated by the improved landslide susceptibility model in 2017–2019 shows that high and very high susceptibility areas decreased rapidly compared with the coseismic landslide susceptibility in 2017. Their total proportion decreases from 32.61% to 11.00%. The very low susceptibility area quickly grows, and its proportion increases from 13.56% to 51.87%. Overall, with the landslide susceptibility from high to low, their proportions are from low to high. The post-seismic landslide susceptibility calculated by the improved landslide susceptibility model in 2019–2020 shows that the high and very high susceptibility areas continue to decrease compared with the post-seismic landslide susceptibility in 2017–2019. Their total proportion decreases from 11.00% to 8.71%. The very low susceptibility area continued to grow, and its proportion increased from 51.87% to 55.57%. Overall, with the landslide susceptibility from high to low, their proportions are from low to high. The changing trend of post-seismic landslide susceptibility calculated by the conventional landslide susceptibility model is similar to that calculated by the improved landslide susceptibility model.

Notably, this study uses the posterior method further to compare the performance of the conventional and improved models. It makes statistics of 2019–2020 post-seismic landslides distributed in every 2017–2019 post-seismic landslide susceptibility classification area (Table 5). As a result of the conventional model, 81.89% of the 2019–2020 post-seismic landslides fell in the high and very high susceptibility areas. As a result of the improved model, 80.07% of the 2019–2020 post-seismic landslides fell in high and very high susceptibility areas. The total area proportion of high and very high susceptibility areas in the improved model is 11.00%, and that in the conventional model is 19.25% (Table 3 and Table 4). Almost the same proportion of post-seismic landslides fell in the high and very high susceptibility areas of the conventional and improved models (Table 5). But, the high and very high susceptibility areas of the conventional model are twice as large as those of the improved model. It means the improved landslide susceptibility model has better-predicting performance than the conventional one.

Table 5. 2019–2020 post-seismic landslides distributed in post-seismic landslide susceptibility classification area of 2017–2019.

Susceptibility	Percentage of classification area (%)
	Conventional model	Improved model
Very low	1.84%	2.52%
Low	5.00%	5.20%
Moderate	11.27%	12.21%
High	21.45%	24.15%
Very high	60.44%	55.92%

3.4. Spatio-temporal characteristics of post-seismic landslide susceptibility

This study analysis spatio-temporal characteristics of post-seismic landslide susceptibility and explores the effects of the different factors on the change characteristics by calculating and ranking the Pearson correlation coefficients of influencing factors (Table 6). Pearson correlation coefficient is also called R-value, R = −1∼1. When P < 0.01, it indicates a significant correlation. The greater the $|R|$ is, the stronger the correlation is, and vice versa.

Table 6. Ranking of Pearson correlation coefficient in improved time-variant landslide susceptibility model.

Ranking	Pearson correlation coefficient
	2017		2017–2019		2019–2020
1	PGA	0.790	Fault	0.708	Distance	0.696
2	Fault	0.777	Distance	0.598	Fault	0.603
3	Strata	0.578	Slope	0.547	Slope	0.593
4	Slope	0.557	Aspect	0.407	PGA	0.428
5	Elevation	0.245	PGA	0.318	Aspect	0.396
6	Aspect	0.163	Strata	0.291	Strata	0.357
7	Road	0.142	Road	0.269	River	0.241
8	River	0.108	Elevation	0.152	Road	0.160
9			River	0.102	Elevation	0.026

The high landslide susceptibility area of 2017 is mainly distributed in areas with high PGA, as shown in Figure 15 and Figure 2 (d). The P-value of every factor is less than 0.01. The R-values of PGA and fault were 0.790 and 0.777, which ranked first and second, respectively. The coseismic landslide susceptibility is significantly related to earthquake-related factors, especially PGA. The high post-seismic landslide susceptibility in 2017–2019 is distributed along both sides of the seismogenic fault, as shown in Figure 13 and Figure 2 (d) (e). In Zhangzha Town, the decline of the post-seismic landslide susceptibility in 2017–2019 was apparent. The P-value of every factor is less than 0.01. The R-values of fault and PGA were 0.708 and 0.318, which ranked first and fifth, respectively. The correlation with PGA had significantly weakened. The R-value of distance was 0.598, which ranked second. The post-seismic landslides susceptibility in 2019–2020 continued to decline. The high and very high landslide susceptibility areas shrunk rapidly towards the concentration area of landslides. Except for the elevation factor, the P-value of other factors is less than 0.01. The R-values of fault and PGA were 0.603 and 0.428, which ranked second and fourth, respectively. The correlation with seismogenic fault also weakened, as shown in Figure 13 and Figure 2 (d) (e). The R-value of distance was 0.696, which ranked first. The correlation between post-seismic landslide susceptibility and distance after the earthquake is the strongest. The result indicates that the correlation between new post-seismic landslides and seismic events gradually weakened with time. The landslide susceptibility was mainly affected by earthquake-related factors such as seismogenic fault and PGA during the earthquake and gradually changed to be primarily influenced by the distance to the pre-phase landslide.

4. Discussion

In this section, we focus first on why there are legacy effects and how to quantify legacy effects. Then, we compare our post-seismic susceptibility model with similar studies in other regions and discuss the implications of our result. We also discuss the replicability of the improved model. In the end, we discuss factors affecting the spatio-temporal characteristics of post-seismic landslide susceptibility.

4.1. Causes and quantification of legacy effects

This study shows that the new post-seismic landslides include the pre-phase landslides expansion and the new single post-seismic landslides. However, the number of new single post-seismic landslides in Jiuzhaigou is not more than 10%, and new single post-seismic landslides distribute in the whole area. Thus, the spatial distribution of new single post-seismic landslides (Figure 16) is scattered. These post-seismic landslides are not very active, and their activities decay rapidly. The pre-phase landslides are activated by rainfall and other factors, and their landslide areas are expanded, which is the main form of post-seismic landslide activity. This form of landslide mainly comes from two types (Figure 17): (1) Since the lower rock mass slides in the pre-phase, the slope gradually forms cracks from bottom to top. The rock mass on the upper part of the slope loses the support of the lower part, and then the rainfall increases the pore water pressure of the cracks. The upper rock mass is easy to slide and forms a retrograde landslide. (2) Coseismic landslide deposits are often loose, fractured, and easy to move under heavy rain in the area. The new landslide also easily occurs at the front edge of the pre-phase landslide. The main increasing form of post-seismic landslides determines the legacy effect of pre-phase landslides on future landslides. The post-earthquake investigations and statistical results show that most of the new post-seismic landslides in each phase are distributed within 200 m away from the pre-phase landslides. Distance is the key factor in quantifying the legacy effect in the area with more pre-phase landslide expansion. This study used the multi-temporal post-seismic landslide and reasonably inputted the distance to the landslide susceptibility model, getting an improved post-seismic time-variant landslide susceptibility model.

Figure 16. Diagram showing the development of the new single post-seismic landslides. (a) Slope without landslide, (b) New single post-seismic landslides (u_w: Pore water pressure).

Figure 17. Diagram showing the development of the pre-phase landslides expansion. (a) Coseismic landslide, (b) Post-seismic landslide caused by cracks propagation, (c) Post-seismic landslide caused by deposits movements, (d) Post-seismic landslides caused by cracks propagation and deposits movements (u_w: Pore water pressure).

4.2. The post-seismic susceptibility model compared with similar studies

Chen et al. (2019) studied the Wenchuan post-seismic landslide susceptibility in 2009 and 2011. They graded the activity intensity of landslides, calculated the susceptibility of landslides with different activity intensities, and then used Spatial Overlay to obtain the final susceptibility results. 70.2% of the post-seismic landslides fell in the high and very high susceptibility areas in 2009, whose area accounts for 32.1% of the total area of the study area. 86.1% of the post-seismic landslides fell in the high and very high susceptibility areas in 2011, whose area accounts for about 50.0% of the total area of the study area. Although over 70% of landslides fell in high and very high susceptibility areas, the high and very high susceptibility areas accounted for a large proportion.

Samia et al. (2017b) made an effort to predict the follow-up landslides as a function of the geometry and location of a pre-phase landslide. A combination of geometric attributes and topographic attributes of landslides used in logistic regression increased the capability of predicting of occurrence of follow-up landslides (AUC = 0.64 ± 0.02). The accurate geometric attributes of pre-phase landslides are very vital in this method. With the deepening of the study, Samia et al. (2020) tried to use a space–time clustering measure derived from Ripley’s space–time K function to quantify path dependency (legacy effect) among landslides. They compared the performance of the conventional, conventional plus path-dependent, and purely path-dependent landslide susceptibility models and found that the best-performing landslide susceptibility model was the conventional plus path-dependent model. The results demonstrated that adding landslide path dependency variables to such an improved conventional landslide susceptibility increased its performance (from an AUC value of 0.771 to an AUC value of 0.801). Although the latest method only uses point position information, the calculation process is still not simple enough.

In this study, the total area proportion of high and very high susceptibility areas in the improved model is far less than that of the conventional model, which is only half that of the conventional model. But, the proportion of post-seismic landslides that fell in the high and very high susceptibility areas of the conventional and improved models exceeds 80%. It means the improved landslide susceptibility model has better-predicting performance than the conventional one. The improved landslide susceptibility model is robust and has higher accuracy.

The improved model in this study quantifies the through simple buffer analysis, which is easier than the existing post-seismic landslide susceptibility methods. Furthermore, the improved landslide susceptibility model increases the number of landslide datasets that can be used. The multi-temporal post-seismic landslide dataset was obtained by automatic classification technology in this study. The improved landslide susceptibility model can use this post-seismic dataset without further processing. Therefore, the improved time-variant susceptibility model suits more post-seismic multi-temporal landslide datasets. It is easier to implement than the existing post-seismic time-variant landslide susceptibility methods, which can be applied to more earthquake cases with a larger research area globally.

4.3. The replicability of the improved landslide susceptibility model

The improved post-seismic landslide susceptibility model and influencing factors can be extended to landslide-prone areas with similar geological conditions, especially in areas where the expansion of pre-phase landslides is the main form of post-seismic landslide activity. To gain confidence in this method, conducting a spatio-temporal characteristics analysis of the multi-temporal landslide inventory in these areas needs to be questioned to find out the conditions for using this method. This study first analyzes the spatio-temporal characteristics of post-seismic landslides using Buffer analysis. One of its purposes is to analyze the main form of post-seismic landslide activity and determine whether the expansion of pre-phase landslides is the main form. Although the Buffer analysis method is simple, it is effective and can quickly analyze the main form of post-seismic landslide activity. Suppose the post-seismic landslides are not concentrated around the pre-phase landslide. In that case, there are more ‘new single landslides’, and the improved model of this study is not suitable for calculating post-seismic landslide susceptibility in this case.

However, the enormous energy released by an earthquake far exceeds the external forces after the earthquake. If an area did not slide during the earthquake, although it may have potential damage and may slide in the future due to external forces, this probability is not very high. In contrast, if the geological conditions of these areas are more similar to those of landslides that slid during earthquakes or nearby landslides, and they suffer greater damage (such as cracks) from the earthquake and the landslides, the probability of landslides is very high in these areas. For example, retrogressive landslides are abundant in earthquake-struck regions. Hence, the post-seismic landslides mainly come from pre-phase landslide expansions. In our other study, we used the 2008–2018 Wenchuan landslide dataset published by Fan et al. (2019b) for statistical analysis, and the spatio-temporal characteristics of post-seismic landslides in Wenchuan are very similar to those in Jiuzhaigou (Guo 2022). The improved model also performs well in the study of Wenchuan post-seismic landslides susceptibility. Therefore, the improved model is also suitable for most post-seismic landslide studies.

4.4. Factors affecting the spatio-temporal characteristics of post-seismic landslide susceptibility

It has been found that seismogenic fault and PGA are the most important influencing factors for coseismic landslides (e.g. Keefer 2000; Qi et al. 2010; Zou et al. 2022). This study confirms this point again. Meanwhile, we noticed that most post-seismic landslides occurred around the pre-phase landslides. This phenomenon is because post-seismic landslides mainly come from the propagation of cracks and the movement of deposits (Figure 16). The cracks and the deposits are significantly related to the pre-phase landslides and distributed around the pre-phase landslide. Therefore, the correlation between post-seismic landslide susceptibility and the distance to the pre-phase landslide is significant. The distance to the edge of the pre-phase landslide plays a vital role in the susceptibility analysis of the post-seismic landslide. Distance gradually becomes the important factor affecting the spatio-temporal characteristics of post-seismic landslide susceptibility and manifests as the high and very high landslide susceptibility areas shrink towards the concentration area of landslides.

Although post-seismic landslide is mainly triggered by precipitation, the definition of susceptibility does not include triggering factors. The fragmented geological environment is still the main influencing factor of the spatio-temporal characteristics of post-seismic landslide susceptibility. The precipitation of the Jiuzhaigou area in 2017, 2019, and 2020 was 786.2, 768.8, and 846 mm, respectively, and that of 2018 was abandoned due to not being updated in the annual hydrological report. After the earthquake, the precipitation in the Jiuzhaigou area was relatively stable. The precipitation in 2020 was the largest, but it was still within the normal precipitation range, and the fluctuation was not more than 5.5% compared with the annual average precipitation. The post-seismic landslide susceptibility of 2019–2020 remains on a decreasing trend. Therefore, although the post-seismic landslides in the Jiuzhaigou area are related to rainfall, the most crucial cause of the post-seismic landslides is still the fragmented geological environment caused by the earthquake. In this study, the distance factor is a proxy for the fragmented surface environment, so it is a crucial factor affecting the spatio-temporal characteristics of post-seismic landslide susceptibility.

In addition, We infer that the correlation between influencing factors, including distance, and the spatio-temporal characteristics of post-seismic landslide susceptibility will change over time. However, due to the limited monitoring time length available to the author, this study did not obtain the changing trend of correlation. It still needs to be further studied in different monitoring time lengths.

5. Conclusions

This study developed a multi-temporal post-seismic landslide dataset and established an improved post-seismic time-variant landslide susceptibility model. The spatio-temporal characteristics of post-seismic landslide susceptibility are generalized, and we draw the following conclusions:

1. The post-seismic landslides are mainly caused by the expansion of pre-phase landslides. The new post-seismic landslides are mainly distributed around the pre-phase landslides in the Jiuzhaigou area. The spatial distribution of new post-seismic landslides significantly became closer to the pre-phase landslide with time.

2. Distance is the critical factor in measuring the legacy effect of the pre-phase landslides on future landslides. It can improve the assessment accuracy of the post-seismic landslide susceptibility model and reasonably calculate the post-seismic landslide susceptibility.

3. The coseismic landslide susceptibility is mainly affected by earthquake-related factors such as seismogenic fault and PGA. The correlation between post-seismic landslides and earthquakes gradually became weaker with time.

In general, post-seismic landslide prevention should focus on the expansion of pre-phase landslides, especially the large landslides triggered by rainfall. Moreover, it should clean up the landslide deposits in time and reasonably dredge the debris flows to avoid the secondary geological disasters.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The data that support the findings of this study are available from the corresponding author, upon reasonable request.

Additional information

Funding

This work was supported by the Second Tibetan Plateau Scientific Expedition and Research Program [grant number STEP 2019QZKK0901]; Integrated ‘Space-Ground’ of Jiuzhai Valley World Natural Heritage Site [grant number E0H20501].