The study area, a UNESCO world heritage site known as the Jiuzhaigou National Geopark, is situated in the eastern edge of the Tibetan Plateau. The site with an area of approximately 653 km² is located at 32°54′21″ ∼33°16′9″ N and 103°46′24″∼104°3′54″ E in the northern part of Sichuan province, southwestern China (see Figure 1A). It is a tectonically active hilly mountainous region characterized by narrow and steep valleys with elevations ranging from 1429 to 4865 m a.s.l. (above the sea level). The annual precipitation in the region is about 704.3 mm and peaks from June to September, accounting for approximately 70% of the yearly precipitation. This region lies between two principal faults: the nearly NS-trending Minjiang Fault and the NWW-trending left-lateral strike-slip Tazang Fault (Fan et al., 2018b; Wu et al., 2018) (Figure 1B). All of these active structures have the potential for triggering earthquakes of magnitude M _L ≥ 7.5 (Li et al., 2016). According to historical earthquake records, this region was affected by seven strong earthquake events in the last century, including the 1933 Diexi M _w 7.3 earthquake, 1960 Songpan M _w 6.3 earthquake, 1973 Songpan M _w 6.1 earthquake, 1974 Songpan M _w 5.7 earthquake, and 1976 Songpan–Pingwu earthquake swarm (M _w 6.9 occurred on 16 August, M _w 6.4 on 21 August, and M _w 6.7 on 23 August) (see Figure 2 in Luo et al., 2021). These events may have caused serious damage to the rock masses, increasing the probability of slope failure. More specific information about the study area can be found in the studies by Luo et al. (2021) and Wang et al. (2018), as well as references therein.

The geological map at a scale of 1:50,000 was provided by the Jiuzhaigou Management Bureau and modified from the published research studies (Cui et al., 2005; Hu et al., 2019). Outcrops in the study area can be grouped into 16 lithological units, and loose Quaternary materials locally covered the Triassic (T), Permian (P), Carboniferous–Permian (Cp), Carboniferous(C), and Devonian(D) rocks (see Table 1 and Figure 1C).

It is a crucial pre-requisite phase to build a trustworthy and accurate landslide inventory map for landslide analysis (Guzzetti et al., 1999; Li et al., 2022; Tanyas and Lombardo, 2020). In this study, the imagery for interpretation was chosen based on the year of acquisition, coverage, minimum of clouds, and resolution to detect multi-temporal landslides. In particular, for co-seismic landslides, it is rather necessary to collect the pre-seismic images that are closer to the earthquake event occurrence date in order to exclude the pre-earthquake landslides and misjudged zones (Harp et al., 2011). This study mapped multi-temporal landslides via a set of remote sensing images from Spot-5, Gaofen, Landsat, unmanned aerial vehicle (UAV), and Google Earth platform, followed by field investigation for the temporal span from August 2017 to September 2020. The spatial resolutions of images range from 0.2 to 3.0 m. We refined landslide signatures through the systematic visual interpretation in the geographic information system (GIS) environment. Table 2 summarizes the sources and characteristics of data used for landslide interpretation and susceptibility modeling phases. Satellite images around Panda Lake used for multi-temporal landslide identification and some typical landslide photographs in the study area are shown in Figure 2.

As a result, we created detailed polygon-shaped multi-temporal landslide maps that record the location, surface area, and types of pre-, co-, and post-seismic landslides. Unfortunately, landslide initiations were not separated from the runout and deposition areas. In general, it requires conducting a polygon-to-point conversion for each landslide to build a statistical susceptibility assessment model. Previous studies offer three categories of sampling strategies (Hussin et al., 2016), and extracting the points populating the scar on a buffer distance is a common choice. Specifically, we first extracted the highest and lowest elevations of each landslide polygon and then calculated the 20% of the height difference as the buffer distance to the most elevated location of the scar, producing landslide identification points (LIPs, as shown in Figure 1B). The types of landslides are mainly small-medium scale shallow debris/rock slides, rock avalanches, and rockfalls (Hungr et al., 2014).

The size distribution and morphometric characteristics of the mapped landslides are illustrated in Figure 3. The pre-EQ landslide dataset contains landslides covering an area of 3.83 × 106 m² in total, and the size varied from A _Lmin = 940.94 m² to A _Lmax = 30.98 × 104 m² (µ = 2.49 × 104 m², σ = 4.24 × 104 m², see Figure 3A), spanning the scale from small to moderate. Subsequently, we mapped 1022 co-seismic landslides with a total covering area of 3.88 × 106 m², ranging in size from A _Lmin = 10.21 m² to A _Lmax = 6.47 × 104 m² (µ = 0.37 × 104 m², σ = 0.62 × 104 m², Figure 3B). The 364 landslides (post-EQ, covering an area of 2.11 × 106 m²) range in size from 3.89 × 102 m² to 5.39 × 104 m² (µ = 0.58 × 104 m², σ = 0.62 × 104 m², Figure 3C) after the earthquake occurred to the present. Subsequently, for each landslide dataset, we generated a corresponding temporal model with a set of covariates in the Binary Logistic Regression Model.

The aspect ratio (also called the length–width ratio, L/W) describes the planar shape of a single landslide (Tian et al., 2017). According to values of this parameter, the landslides were categorized as 1) the transverse landslide (LA1, L/W ≤ 0.8), 2) the isometric landslide (LA2, 0.8 < L/W ≤ 1.2), 3) the longitudinal landslide (LA3, 1.2 < L/W ≤ 3), and 4) the elongated landslide (LA4, L/W > 3). An aspect ratio close to 1 means a typical rotational slide, and with the increase of the aspect ratio, it is a much more elongated shape typical of flow-type landslides. According to the statistics, the minimum and maximum aspect ratios for pre-EQ landslides (Figure 3D), co-seismic landslides (Figure 3E), and post-EQ landslides (Figure 3F) are 2.09, 0.29, and 2.11 and 12.03, 84.08, and 11.36, respectively. The average value equals 3.66, 5.04, and 3.34, respectively. Moreover, it is similar that all landslides for three datasets are mainly composed of longitudinal (LA3) and elongated landslides (LA4). Specifically, LA3 accounted for 42.86%, 34.92%, and 54.67% and LA4 for 57.14%, 54.87%, and 45.33% for pre-, co-, and post-seismic landslides, respectively. This elongation resulted, in general, from moderately long run-out distances down steep slopes below landslide source areas.

In the vast literature reports (e.g., Carrara et al., 1995; Guzzetti et al., 1999; Reichenbach et al., 2018), the most frequent mapping units for large-scale or regional landslide susceptibility assessments are grid cells and slope units (SUs, hereafter). SUs are the basic units of landslide occurrence generated according to the terrain hydrology unit bounded by drainage and ridgelines (Carrara et al., 1995). Compared to the grid cells, SUs offered two main advantages: one for a smaller computational burden and the other for the geomorphologically oriented interpretation (Steger et al., 2016). Because it has higher internal homogeneity and between-unit heterogeneity, SUs have been proven to be reliable, as tested in some cases by Alvioli et al. (2016), Amato et al. (2019), and Lombardo et al. (2020). In this study, the r.slopeunits module (Alvioli et al., 2016) was used to automatically delineate SUs from an input DEM, and the study area was partitioned into 1234 SUs. Figure 4A portrays the size distribution and statistic characteristics of the generated slope-units. Their areas range from a minimum of 0.05 km² to a maximum of 3.25 km² with a mean of 0.53 km², and a standard deviation 0.43 km², indicating that the local topography is fairly rugged.

For each landslide inventory, we counted the LIPs per SU and transformed the results to binary presence–absence to reflect the presence–absence of the landslide distribution over the study area at the SU level. If an SU included at least one LIP, it was given a positive landslide status of “1,” and if not, it was set as “0” for an absence status. The presence/absence status assignment of the partitioned 1234 SUs for each landslide inventory is shown in Figure 4B. Consequently, 96 SUs (i.e., ∼ 7.78% of the total 1234 SUs) were assigned as “1,” while the remaining 1138 SUs were allocated as “0” for the pre-EQ landslide dataset. Regarding the co-seismic landslide inventory, 173 SUs accounting for approximately 14.02% of the total SUs were classified as presence, while the rest were classified as absence, which accounts for 86% of all SUs. For the post-EQ landslides, 163 SUs were assigned as “1,” accounting for 13.21% of the total (N = 1234) SUs.

In this work, we opted for a probabilistic model as the binary formulation of the logistic regression model. It is a generalized linear model belonging to the statistical family and widely used in the vast majority of slope failure studies for a given mapping unit to explain the distribution of landslides over space (Brenning, 2005; Reichenbach et al., 2018). More specifically, it aims at modeling a linear relationship between covariates and landslide occurrences based on a set of predictors, which can correspond either to continuous or categorical variables in the model. The BLR model assumes a Bernoulli probability distribution as the underlying stochastic process, and it can be stated as follows:pt η ( x ) = In [ π ( x ) 1 − π ( x ) ] = β 0 + β 1 x 1 + ⋯ + β j x j , (1) where η is the logic link, π(x) is the conditional probability of the landslide occurrence, given the predictors, β ₀ is the constant intercept, and β ₀, … …, β _j are regression coefficients measuring the effects of each covariate. Once the model has estimated the linear predictor, the following logic link η is used to obtain the probability of the landslide occurrence P: P = e β 0 + β 1 x 1 + ⋯ + β j x j 1 + e β 0 + β 1 x 1 + ⋯ + β j x j . (2)

The large covariates’ hyperspace we have built may posess the potential multi-collinearity issue, which would have a negative impact on our statistical models (Kelava et al., 2008). To avoid this, we performed a multi-collinearity test and utilized Pearson’s correlation coefficients r to measure how closely the variables are linearly related. In accordance to Moore et al. (2006), |r| ≤ 0.8 stands for an extremely weak collinearity that can be considered with barely no linear dependency between the pairs of covariates. Figure 6 is the Pearson’s correlation coefficient matrix obtained in the software R. Overall, not even one case exceeded the thresholds (|r| = 0.8) for significant dependency although some moderate pair-wises existed, including these moderate positively correlated pairs: ELσ-TWIσ (r = 0.71), ELµ-Dis2roadsµ (r = 0.63), ELσ-SLµ (r = 0.61), ELσ-Dis2roadsσ (r = 0.62), PRCµ-TWIµ (r = 0.68), RSPµ-RSPσ (r = 0.62), Dis2roadsσ-Dis2riversσ(r = 0.63), Dis2roadsσ-Dis2seisfaultsσ (r = 0.61), and Dis2faultsσ-Dis2seisfaultsσ (r = 0.61) and negatively correlated pairs including PRCµ-RSPσ (r = -0.64) and TWIµ-RSPσ (r = -0.65). Ultimately, no conclusive estimation was pointed out for excluding the covariates originally reported in Table 3. In other words, we incorporated all covariates into the Eqs 1, 2 for the following landslide susceptibility analyses.

Receiver operating characteristic (ROC) curves and their integrated area under curves (AUCs) can be used to assess model performance, including the fitting and prediction skill (Hosmer and Lemeshow, 2000). The ROC is obtained by plotting the true positive rate (TPR, Eq. 3) against the false positive rate (FPR, Eq. 4). The AUC value spans from 0.5 (random guess) to 1 (ideal fit or predictions), and the overall model performance can be recognized as a function of AUC values and classed as follows (Hosmer and Lemeshow, 2000): 1) 0.7 < AUC ≤0.8, acceptable; 2) 0.8 < AUC ≤ 0.9, excellent; and 3) AUC>0.9 outstanding performance. T P R = T P T P + F N , (3) F P R = F P T N + F P , (4) where false positive (FP)/false negative (FN) represents the number of misclassified landslide/non-landslide samples, and true positive (TP)/true negative (TN) represents the number of landslide/non-landslide classified correctly.

The regression coefficients of those examined covariates obtained for each model we built with the multi-temporal landslide inventory, that is, pre-seismic Mod1, co-seismic Mod2, and post-seismic Mod3 are shown in Figure 7. A total of 12 variables were significant in at least one of the three models out of all the 26 continuous covariates (see Table 3). Due to each of the continuous variables being rescaled with mean zero and unit variance, the covariate effects reported in Figure 7 are all on the same scale and directly comparable. The regression coefficient values above the zero line indicate a significant positive contribution to landslide occurrence, and symmetry values below the zero line correspond to significant negative effects.

As shown in Figure 7, the continuous covariates’ (see Figure 7A) and categorical fixed effects (see Figure 7B) of each model are reported by the mean regression coefficient and the associated 95% credible interval. Overall, some covariates play a dominant role in the model, and the majority of covariates have a limited credible interval. Specifically, for Mod1, the most relevant covariate is the mean slope, which positively contributes to the landslide occurrence (we utilized the mean regression coefficient β to describe the covariates’ effects, hereafter, β _SLµ = 0.665). In contrast, Dis2roadsμ shows a negative contribution (β = −0.809). For the co-seismic Mod2, the mean of PGA appears to dominate the susceptibility pattern with the contribution (positive, β = 3.199) in the absolute value much larger than any other covariate, whereas Dis2roadsµ (β = −1.574) shows a similarly negative contribution to the landslide with Mod1. With the post-seismic Mod3, the mean Dis2seisfault (negative, β = −1.366) and elevation (negative, β = −1.417), in the absolute value, show the greatest degree of negative influence for the slope instability. In all three models, some covariates appear to play the same role, such as the regression coefficient of the mean elevation per SU appears to be negative. The possible reason is that most of the landslides were recognized in the lower portions of the topography in the study area (see Figure 1B). The SLµ positively contributes to the slope instability (β _Mod1 = 0.665, β _Mod2 = 0.713, and β _Mod3 = 0.526); this is consistent with the general knowledge that landslides generally trigger in the steep terrain. The mean of northness appears to have a similar negative effect on the landslide occurrence, which may be due to the decrease in the solar radiation exposure with the increase in northness (β _Mod1 = −0.305, β _Mod2 = −0.246, and β _Mod3 = −0.309). In terms of the mean Dis2faults (β _Mod1 = −0.517, β _Mod2 = −0.736, and β _Mod3 = −0.073), it shows a negative contribution to the landslide occurrence because a greater extent of weathered and fractured rock masses around the faults would decrease the slope stability. Ultimately, RSP, TWI, and PRC showed no significance in all models.

For all three models, lithology does not appear to be significant because the zero line crosses the distribution of each categorical class (Figure 7B). For Mod1, two classes show the slightly significant for the landslides, that is, Cpd₂ chert limestone, bioclastic rock (negatively, β = −0.179) and Cm₂ bioclastic limestone, dolomitic limestone with chert limestone (positively, β = 0.102), while the other classes to the susceptibility are negligible. For Mod2, two classes show the slightly significant for the landslides, Quaternary eluvium, alluvial and diluvial deposits, colluvium and diluvium deposits (negatively, β = −0.193), and Pds₂ limestone intercalated with metamorphic siltstone and slate (positively, β = 0.170). For Mod3, Cpd₁ limestone mixed with clastic limestone (positively, β = 0.277) and Dcy bioclastic limestone and chert-bearing limestone with gravel limestone (negatively, β = −0.349) show the most relevant to the landslides out of these outcropping lithology classes.

Table.4 and Figures 8A–F show the derived landslide susceptibility results (mean probability and the uncertainty). The equal spacing method was used to reclassify the mean susceptibility into five classes: very low (0∼0.2), low (0.2∼0.4), moderate (0.4∼0.6), high (0.6∼0.8), and very high (0.8∼1). According to Table 4, six SUs are highly susceptible to landslide occurrence, whereas 3 years after the Jiuzhaigou earthquake, 57 SUs (accounting for 4.62% of the total SUs) were highly susceptible to instability. The number of slopes with very high susceptibility increased from only one (pre-seismic) to 21 (post-seismic). For the pre- and post-earthquake periods, the very high landslide susceptibility areas are mainly distributed along the valleys, especially the Shuzheng valley and Zechawa valley for the pre-earthquake period, whereas the Rize valley and Panda Lake, which are closer to the epicenter, for the post-earthquake period. In contrast, the co-seismic susceptibility map shows clustered around the epicenter of the earthquake and a southward decreasing trend. Moreover, the pre- and post-seismic uncertainty maps follow a similar pattern to the mean susceptibility maps. In the case of co-seismic landslides, the uncertainty associated with the mean susceptibility appears to be relatively small in the southernmost sector of the study area, although it has a considerably greater spread to the north.

The error plot (Figures 8G–I) shows the mean susceptibility against its associated uncertainty measured in a 95% confidence interval (CI) per SU of each space-time model. It is crucial for determining if the estimates of constructed models are reasonably acceptable in the landslide susceptibility studies (Rossi et al., 2010). The left and right tails of the mean probability distribution should be associated with a very limited uncertainty for an ideal model (Reichenbach et al., 2018). Compared to the pre-seismic (Figure 8G) uncertainty map, the error plots of co-seismic (Figure 8H) and post-seismic (Figure 8I) maps are much better determined, and this is also proved by the model performance shown in Model Performance.

A performance evaluation step is required for each landslide susceptibility assessment model. Figure 9 shows the ROC curves and AUC values of calibration and 10-fold cross-validation results for each model. Although all models produce excellent and outstanding results according to the Hosmer and Lemeshow (2000) classification, Mod2 appears to be outperforming the other two models both in terms of estimated fitting and cross-validation results. More specifically, the fitness results of AUC values for Mod1, Mod2, and Mod3 are 0.838, 0.956, and 0.910 (Figure 9A), respectively, corresponding to be excellent and the rest two for outstanding goodness-of-fit. Also, the validation scheme (see Figures 9B–D for ROCs and associated AUCs distribution) generally performs worse than the corresponding fitting in Figure 9A; the average AUCs 0.801, 0.942, and 0.880 of Mod1, Mod2, and Mod3, respectively, showed the excellent prediction capability for potential slope failure. Our results show that the BLR model has good applicability in predictive performance with outstanding AUC results both for non-seismic, co-seismic, and post-seismic landslide datasets, which is consistent with the grid-based Jiuzhaigou co-seismic landslide susceptibility assessment outcomes in Fan et al. (2018b) (AUC = 0.851) and Ma et al. (2019) (AUC = 0.89). As the validation database changes, a limited variability between the minimum and maximum AUC in the resulting 10 AUCs is observed, that is, 0.069 (AUC_max = 0.976, AUC_min = 0.907) for Mod2 and 0.15 (AUC_max = 0.880, AUC_min = 0.809) for Mod3, although this cannot be said same for the pre-seismic Mod1 which has a significant variability (0.315, AUC_max = 0.926, AUC_min = 0.611). A minor variability of AUCs indicates a more robust predictive model performance that did not change much with the training and test subset replicates in a cross-validation scheme.