Investigation on the spatial distribution of landslides in Sichuan Province, southwest China

Abstract

In this study, a total of 34,893 landslides collected from the Sichuan Province of China are used to reveal their spatial distribution. Correlations between the spatial distribution of landslides and five main factors [i.e. topography, rainfall, lithology, peak ground acceleration (PGA), and human activity] are analyzed based on the Geographical Information System (GIS) platform. Topographic factors exert a significant influence on landslide occurrence, and other factors also sensitive to landslides include lithology (Classes 3 and 4), PGA (≤0.15 g), rainfall (90–140 mm), and road density (≤0.6 km/km²). The geographical detector model is used to detect the relative importance of individual influencing factors and the factor interaction effects on the landslide distribution. The results show that the interacting factors between rainfall and PGA, lithology, river density, and road density are more closely related to the spatial distribution of landslides than every single factor. The results of this study can provide a foundation and insights for hazard reduction schemes in the future.

Keywords:

• Landslides
• spatial analysis
• geographical detector model
• Sichuan Province

1. Introduction

Landslides are one of the widespread geological hazards around the world and can cause numerous casualties, serious economic losses, and substantial environmental damage (Li et al. 2021; Zhao et al. 2021; Franceschini et al. 2022; Tanyaş et al. 2022). Landslides are often accompanied by a series of major hazards in the influenced areas. Examples include transportation disruptions, destruction of agricultural land and factories, river blockages, and even buries villages and people. It is a major constraint on the economic construction and development of local areas (Chen et al. 2021; Zhou et al. 2022). Under the influence of strong tectonic activities, frequent climate changes, and human activities, the losses caused by landslides have an increasing trend, especially in mountainous areas (Huang et al. 2021). With the continuous increase in population and the development of the social economy in mountainous areas, the exploitation of the environment by human beings is increasing, which seriously destroys the geological environment and increases the probability of landslide occurrences (Jiang et al. 2022).

Approximately two-thirds of China’s territory is located in mountainous areas, making it one of the countries with the highest frequency of geological disasters in the world (Lin and Wang 2018). The statistics from the China Geological Environment Monitoring Institute show that the number of geological disasters in China is recorded at 292,290, which resulted in 8672 total deaths, with ∼542 fatalities per year during 2006–2021 (Figure 1). The direct economic loss was ∼9.004 billion US$. In particular, landslides are considered to be the most common type of geological hazard in China’s mountain areas, accounting for 72.67% of all events. Sichuan Province in China is highly exposed to landslide hazards because of its unique geomorphological, geological, and climatic characteristics (Cui et al. 2021; Sun et al. 2021; Yi et al. 2021; Xia et al. 2022). In recent years, many disastrous landslides have occurred in Sichuan Province. For example, on 10 July 1989, a catastrophic landslide occurred in the Huaying Mountains in Guang’an City, Sichuan Province, that had a volume of 0.2 million m³, destroyed four villages, and killed 221 people (Tang et al. 2021). The direct economic loss was ∼0.75 million US$. On 5 September 2004, a large-scale landslide was triggered by heavy rainfall in Tiantai Town of Xuanhan County, Sichuan Province (Huang 2009). The landslide destroyed 2983 houses, making 1255 people homeless and causing a direct economic loss of more than 14.5 million US$. In 2017, a large catastrophic landslide with a volume of >18 million m³ occurred in Xinmo village, Sichuan Province, which claimed 83 human lives (Chang et al. 2022). On 14 August 2019, the Ganluo landslide in the Chengkun Aidai railway station, Ganluo County, Sichuan Province, killed 17 people with an accumulation volume of ∼0.05 million m³ (Zhu et al. 2020). On 21 August 2020, a disastrous landslide with a volume of >5 million m³ occurred in Hanyuan County, Sichuan Province, destroyed eight houses, and caused nine fatalities (He et al. 2021).

Figure 1. Graph showing the number of recorded geohazards, landslides, the amount of economic loss, and the number of fatalities for the period 2006–2021.

The spatial distribution of previous landslides can indicate the likelihood of future landslides within a specific area (Zhuang et al. 2018; Xu et al. 2020; Shao et al. 2023). A considerable number of studies on landslide distribution have been reported in the scientific literature. For example, Peng et al. (2019) investigated the regularity of the regional distributions of loess landslides and the governing factors of their occurrences and found that the Chinese Loess Plateau can be split into eight zones with high potential for landslides that are classified into six genetic categories. In the Bailong River Basin of the NE Qinghai-Tibet Plateau, Qi, Meng, et al. (2021) carried out a systematic analysis of large landslides. They discovered that fault zones control the spatial distribution and geomorphic characteristics of large landslides. Qi, Zhao, et al. (2021) studied the correlation between landslide geographical distribution and causative factors. The results showed that the landslide distribution on various scales was closely correlated with the annual precipitation index, average slope, and vegetation coverage index. Wang, Clague, et al. (2021) discussed the association between landslides and the strength of rock masses. The results revealed that rock mass strength controls the stability of slopes, but the major failure causes are rock uplift and river incision, rather than a decrease in rock strength, as has been suggested in some tectonically passive areas. Understanding the spatial distribution of landslide occurrence is a prerequisite for landslide risk mapping and susceptibility assessment. However, it is challenging to evaluate the spatial distribution of landslides since, based on current landslide classifications, the spatial characteristics of historical landslides are poorly understood. Furthermore, the factors that affect the distribution of landslides in different regions may vary substantially, and thus, scientific and technical efforts are required to bridge the gap.

Over the years, the landslides distributed in Sichuan Province have been extensively explored by several scholars. Research works have mainly focused on the triggering mechanisms, characteristics, susceptibility assessments, and stability analysis of landslides (Cao et al. 2019; Yang et al. 2019; Li et al. 2020; Wang et al. 2020; Wu et al. 2021; Hou et al. 2022; Ji et al. 2022; Peng et al. 2022; Yao et al. 2022). However, these findings are difficult to apply to other landslide disasters in Sichuan Province. Environmental conditions influencing landslide occurrence include topography, geomorphology, lithology, geologic structure, tectonic activity, rainfall, earthquakes, and anthropogenic activity, many of which are specific to a particular area (Pudasaini and Krautblatter 2021; Wang, Schweizer, et al. 2021; Troncone et al. 2022). Identifying the various factors influencing landslide distribution and occurrence not only facilitates a comprehensive understanding of landslide mechanisms but also provides a scientific basis for landslide vulnerability, hazard, and risk assessment. However, there is a lack of systematic research on the statistical relationships between geoenvironmental factors and the spatial occurrence of landslides in the entire Sichuan Province. Geographical information systems (GISs) are excellent at managing and analyzing landslide data from a wide range of sources. GIS in landslide research has the potential to enhance our ability to quantitatively and spatially evaluate landslide disasters.

Characterization of the spatial distribution of landslides and their influential factors is not only meaningful but also fundamental for geomorphological studies and natural hazard evaluation. Documenting landslide occurrence in Sichuan Province and deducing the factors that are related to mass wasting are prerequisites for the assessment of regional hazards and risks in Sichuan Province. In this study, we attempt to analyze the spatial distribution of landslides in the Sichuan Province, with the aim of understanding the relationships between possible influencing factors and the distribution and development of landslides. The results of this study can provide more fundamental information for forecasting potential risks in landslide regions.

2. Regional setting

Sichuan Province, located in southwest China, occupies an area of 486,100 km² and lies between longitudes of 97°21′ to 108°31′ E and latitudes of 26°03′ to 34°19′ N (Figure 2). The Yangtze River and its tributaries (such as the Jinsha River, Yalong River, and Minjiang River) run through Sichuan Province. The geomorphology and topography of the study area are highly diversified and are composed of basins, hills, mountains, and plateaus (Yang et al. 2020). The study area is bordered by the Tibetan Plateau in the west, the Daba Mountains, Min Mountains, and Bayankala Mountains in the north, and the Yunnan-Guizhou Plateau in the south. The altitude of the study area varies from 88 to 7845 m above sea level, with an average altitude of 3750 m. The slope angles in the area range between 0 and 89°, with an average of 44°. The study area is strongly influenced by the Southeast Pacific Ocean monsoon and the Southwest Indian Ocean monsoon. Sichuan Province can be divided into three major climates, including the mid-subtropical humid climate of the Sichuan Basin, the subtropical and semihumid climate of southwestern Sichuan, and the alpine climate of the northwestern Sichuan Plateau. The study area belongs to a temperate continental monsoon climate, being warm in summer and cool in winter. The average temperature in the region ranges from 5 (January) to 27 °C (July). According to the Sichuan Provincial Meteorological Service (http://sc.cma.gov.cn/), the average annual rainfall ranges from 900 to 1200 mm, decreasing from east to west. The major fault zones in the study area include the Longmenshan fault zone, the Jinsha River fault zone, and Daliangshan fault zone, and the Minjiang fault zone. The seismic activity of the study area is relatively strong, and 186 large earthquakes of Ms ≥ 5.0 (such as the 1973 Ms 7.6 Luhuo earthquake, the 2008 Ms 8.0 Wenchuan earthquake, the 2013 Ms 7.0 Lushan earthquake, and the 2017 Ms 7.0 Jiuzhaigou earthquake) have been recorded in this area (Figure 3).

Figure 2. Location map of the study area.

Figure 3. Distribution of large historical earthquakes in Sichuan Province.

3. Data and methods

3.1. Data

The landslide distribution data in Sichuan Province were acquired from the Department of Natural Resources of Sichuan Province, a government institution primarily responsible for the monitoring and investigation of geologic hazards. A total of 34,893 landslides were identified, including occurrence times, latitude and longitude coordinates, landslide names, and areas and volumes (Figure 4). Figure 5 shows eight typical landslides in the study area. Actual landslide occurrence is controlled by a variety of environmental variables, such as topography, lithology, rainfall, earthquakes, and human activity. To analyze the effects of these variables on the spatial distribution of landslides, fourteen factors (Figure 6), including slope angle, elevation, slope aspect, cutting depth, river density, relief amplitude, surface roughness, stream power index (SPI), plan curvature, profile curvature, road density, rainfall, PGA, and lithology (Table 1), were considered for the statistical analysis. These factors can be divided into five categories: topography (slope angle, aspect, elevation, relief amplitude, cutting depth, SPI, plan curvature, profile curvature, river density, and surface roughness), climate (rainfall), geology (lithology), seismic factor (PGA), and human activity (road density). The details of the available data are listed in Table 2. The topography data were acquired by processing the digital elevation model (DEM) from ALOS PALSAR with a 12.5-m resolution. Rainfall data can be provided by the Sichuan Provincial Meteorological Service. Google Earth^© images were used to obtain information on roads. The lithology data were extracted by vectorizing geological maps at a scale of 1:200,000.

Figure 4. The spatial distribution of landslides in Sichuan Province.

Figure 5. Photographs of landslides in the study area.

Figure 6. Landslide influencing factor maps: (a) elevation, (b) slope angle, (c) slope aspect, (d) relief amplitude, (e) cutting depth, (f) river density, (g) plan curvature, (h) surface roughness, (i) profile curvature, (j) SPI, (k) rainfall, (l) lithology, (m) PGA, and (n) road density.

Table 1. Description of geological units of the study area.

Group	Code	Lithology	Geological age
Group 1	Q1^fp, Q2^fp, Q2^g, Q3^fl, Q3^fl, Q4^f	Gravel, sand, sandy gravel, clay	Quaternary
Group 2	N, E, R, γ5	Sandy conglomerate, shale, granite	Tertiary
Group 3	K1, K2, r3 5, γη3 5	Sandy conglomerate, mudstone, siltstone, sandstone, granite	Cretaceous
Group 4	J1, J1 + 2, J2, γη2 5, γδ2 5, J3	Mudstone, sandstone, limestone, granite, diorite	Jurassic
Group 5	γη1 5, ηO2 5, δO1 5, γO1 5, T1–2, T2, T3, V1 5, P2T1, δ1 5, γη1 5, CP, γδ1 5, γδ1 4, γ1 5, P1, P2, εO1 5, ε1 5	Monzonite, diorite, granite, gabbro, peridotite, limestone, sandstone, syenite	Triassic-Permian
Group 6	DC, βμ, δO4, σ4, βμ4, γσ4, D	Dolerite, diorite, peridotite, granite, dolomite, limestone, mudstone, sandstone, shale	Carboniferous-Devonian
Group 7	S, OS, S1 + 2, O, AnS	Slate, phyllite, tuff, shale, marl	Silurian-Ordovician
Group 8	Є, Z, Z1, Z2, Pt1, ZЄ	Limestone, dolomite, sandstone, sandy conglomerate, shale	Cambrian-Sinian
Group 9	Pt2, δ2, γη2, γ2, Ar3-Pt1	Diorite, griotte, gneiss, granite, andesite, slate	Pre-Sinian

Table 2. Details of the dataset used in the study.

Dataset	Factors	Scale/resolution	Type
Landslide inventory^a	Landslide	/	Datasheet
Topographical map^b	Plan curvature	12.5 m	Grid
	Profile curvature	12.5 m	Grid
	Slope angle	12.5 m	Grid
	Slope aspect	12.5 m	Grid
	Relief amplitude	12.5 m	Grid
	Cutting depth	12.5 m	Grid
	River density	12.5 m	Grid
	Surface roughness	12.5 m	Grid
	Elevation	12.5 m	Grid
	SPI	12.5 m	Grid
Geological map^c	Lithology	1:200,000	Grid
Seismic map^d	PGA	1:200,000	Vector
Rainfall map^e	Rainfall	1:200,000	Grid
Road map^f	Road density	1:200,000	Vector

Data source: ^aDepartment of Natural Resources of Sichuan Province, ^bALOS PALSAR DEM, ^cChina Geological Survey, ^dChina Geological Survey, ^eSichuan Provincial Meteorological Service, ^fGoogle Earth^© images.

3.2. Methods

To explore how environmental variables influence the spatial distributions of landslides, the landslide concentration index (LCI) is introduced here (Qi et al. 2010; Zhan et al. 2018). The LCI is calculated according to the following equation (Zhan et al. 2018): (1) $$\mathit{LCI}=\frac{L_i/L}{F_i/F}$$(1) where L_i is the number of landslide cells of class i, L is the total number of landslide cells in the study area, F_i is the number of cells of class i, and F is the number of total pixels in the study region.

The geographical detector model (GDM), first proposed by Wang et al. (2010), is a statistical method developed to measure the spatially stratified heterogeneity of an event and assess the association between the event and potential influencing factors. The general assumption of GDM is the following: if an independent variable (X) controls or contributes to a dependent variable (Y), then the spatial distribution of Y should be similar to that of X (Luo et al. 2016). The GDM was first applied to medical geography (Wang et al. 2010) and then to studies of the socioeconomy (Kim et al. 2018), eco-environment (Liang and Yang 2016), and landslide susceptibility assessment (Luo and Liu 2018). One of the advantages of GDM is that both quantitative and qualitative data can be analyzed. The calculation is as follows (Wang et al. 2010): (2) $$q=1-\frac{1}{N{\sigma }^2}\sum _{i=1}^L{N}_i{\sigma }_i^2$$(2) where q is the power of the determinant, N is the number of sample units (grids) in the entire study region, N_i is the number of samples in zone i, σ² is the variance of the entire study region, and σ_i² represents the variance within stratum i. The q-value ranges from 0 to 1, and the larger the q-value is, the stronger the explanatory power of the determinant. When q = 1, Y is completely determined by X. In contrast, when q = 0, Y is completely unrelated to X.

Another advantage of the GDM is the ability to uncover the interaction effect between any two independent variables. The interactive detector describes whether two influencing factors together weaken or enhance each other or whether they affect the dependent variable independently. More detailed descriptions can be found in these studies (Wang et al. 2010; Ju et al. 2016). In this paper, the dependent variable Y can be quantified by the landslide density of Sichuan Province (Figure 7), which can be calculated with the point density tool in ESRI ArcGIS^® 10.2.

Figure 7. Landslide density map of Sichuan Province: (a) seismic landslides and (b) rainfall landslides.

4. Spatial analysis of landslides with possible influencing factors

4.1. Correlation with topography

Plan curvature represents the shapes of the slopes and influences the concentration of soil moisture. The correlation between plan curvature and landslide distribution is displayed in Figure 8(a). As the plan curvature gets closer to zero, the values of landslide abundance proxies (frequency and LCI) decrease. These results indicate that landslides occur more frequently in concave and convex areas. It was not possible to determine a clear effect of plan curvature on the occurrence of landslides for this case, as both convex and concave terrain were similarly influenced by slope failures. Concave slopes, commonly covered by thick layers of loose weathered materials, are prone to instability under the influence of rainfall. The horizontal compressive stress of the convex slopes along the slope surface is weak, and sometimes even tensile stress occurs. The topographic amplification effect of the convex slope is obvious under earthquake action, which is not conducive to the stability of the slope. Profile curvature can determine the driving and resistance forces in the direction of a potential slide. Fifty percent of all landslides occurred at a profile curvature of >0.001, which covered an area of 50.43% (Figure 8(b)). The LCI value is the largest (1.0036) in the < −0.001 class, followed by >0.001 (1.0007) and −0.001–0.001 (0.6968). For the values of frequency and LCI, profile curvature and plan curvature have similar trends.

Figure 8. The relationship between possible influencing factors and landslide abundance proxies (frequency and LCI).

Gentle and medium terrain gradients dominate Sichuan Province, as 75.57% of the study area has slope gradients <31°. A total of 77.25% of the landslides lie in the slope gradient range <24°, which accounts for 59.73% of the gradients in the area. The tendencies of frequency and LCI with slope gradient are similar. The landslide abundance proxies increase with increasing gradient until the maximum appears in the 9–16° category and then decrease with increasing slope gradient (Figure 8(c)). The minimum frequency and LCI occurred at slope angles higher than 47°, and their values were 0.007 and 0.226, respectively. The low frequency of landslides in steeper slopes is mainly because of the small number of such cells. In addition, the steep natural slopes in the study area are mainly composed of bedrock with strong weathering resistance and thus are not prone to landslides.

Figure 8(d) shows the correlation of landslide distribution with the slope aspect. The area of each class occupies a similar area except for the flat class. The two landslide abundance proxies show similar correlations with the slope aspect. Our analysis indicates that nearly 51.81% of all landslides occurred on the slope aspects of the southeast, south, southwest, and west directions. The values of LCI for these directions are >1.0, indicating a higher probability of landslides. Conversely, the LCI values of the remaining slope aspect are <1, indicating a low probability of landslide occurrences. Particularly, for the slope aspect flat, this value is as low as 0.33. Some landslides occurred in the flat area because the centres of the landslide polygons were used for analysis. The pixel size of the DEM may also be partially responsible for this phenomenon.

Variations in relief amplitude are relatively smaller in the eastern part, generally <118 m, where many landslides are seen (Figures 4 and 6(d)). As shown in Figure 8(e), the two landslide abundance proxies decrease rapidly with increasing terrain relief. Values of LCI are >1 when relief amplitudes are <222 m, indicating a higher probability of landslides. When relief amplitudes are >222 m, LCI values are below 1.0, meaning a weak control of terrain on landslide distribution. 84.58% of the landslides are distributed in the relief amplitude range <222 m, which accounts for 63.21% of gradients in the area. The reason why the landslides are mainly concentrated in areas with lower relief amplitudes is that these areas are relatively flat and suitable for human habitation, and the destruction of human activities on the slope is more serious.

Landslides are preferentially distributed on the cutting depth of <275 m, occupying 93.26% of the total landslide number. The cutting depth of 275–382 m has a larger area: 83,820.7 km², amounting to ∼17.20% of the entire study area (Figure 6(e)). However, only ∼5.90% of the landslides occurred in this area. When the cutting depth increases, the number of landslides increases first and then decreases (Figure 8(f)). On cutting depths >0 m and smaller than 188 m, the LCI values are >1, indicating that landslides are more likely to occur on these slopes. For other cutting depth ranges, the LCI values are <1, indicating that areas in these classes were less prone to landslides. Under the cutting effect of flowing water, the slope in the study area becomes higher and steeper, and the stability becomes worse. However, it is not that landslides are more likely to occur when the cutting depth is greater. The increase in cutting depth promotes the continuous increase in slope gradient, which limits the development of landslides to a certain extent because when the slope exceeds a certain slope, the landslide may transform into other disasters. In addition, hillslopes that are excessively steep hinder runoff accumulation and contrarily restrain landsliding by rainstorms.

The correlation between landslide distribution and river density is shown in Figure 8(g). Rivers drive the presence of landslides in the study area, and the occurrence frequency of landslides increases with increasing river density, except for densities >0.83 km/km². A large number of landslides are mainly concentrated at river densities >0.64 km/km², accounting for 89.05% of the total. The maximum LCI value (1.70) occurs at a river density of >0.83 km/km², and the minimum value (0.62) occurs in the <0.38 km/km² range. The direct proportionality between river density and landslide occurrence observed in the study area is a consequence of the terrain’s susceptibility to undercutting of slopes by river erosion. River incisions and bank erosion changed the three-dimensional stress state of slopes, increased the gravitational potential energies of slopes, and finally led to slope instability near the rivers.

The surface roughness of most of the study area is <1.23. Due to the largest areal coverage (410,128.8 km², 84.43%), 94.64% of the landslides fall in the surface roughness of <1.23 (Figure 8(h)). The values of frequency and LCI show an increasing trend with increasing surface roughness and reach their maximum values in the range of 1.0–1.1 and then decrease rapidly in the range of more than 1.1. The LCI values are >1 within surface roughness categories of 1.0–1.1, indicating a higher probability of landslides. The LCI values of the remaining surface roughness are <1, indicating that such slopes are not prone to landslides. Landslides tend to be concentrated in locations with lower surface roughness because the slopes in these areas have weak weathering resistance and the damage caused by human activities to the slope is relatively serious.

Figure 8(i) shows the variation between the concentration of landslides and elevation. The values of frequency and LCI decrease with increasing elevation. The landslides are mostly distributed in the lower and middle elevation areas. The slopes falling in elevation classes of <824 m have the largest landslide abundances, which are 0.6687 in frequency and 2.57 in LCI. Nearly 92.54% of the landslides occurred at elevations of <2378 m. Human activities (such as mining, highway construction, and hydraulic engineering) are mainly concentrated in the lower and middle elevations and have a significant impact on the aforementioned landslide distribution pattern. Less than 0.29% of all landslides occurred at elevations >3803 m. The minimum LCI value (0.001) appears at elevations >4341 m. The landslide concentration drops abruptly at higher elevations. The high-elevation region meets the terrain conditions for the occurrence of landslides, but slippery material is lacking.

The total area with a stream power index (SPI) <6 is ∼398,132.8 km², accounting for 81.96% of the entire area. Approximately 82.79% of the landslides occurred at SPI values of <6. As the SPI increases, the number of landslides decreases (Figure 8(j)). Assessment of SPI showed that the classes <6 and >30 have a high correlation with landslide occurrence. The minimum LCI value appears in the 12–18 range, indicating that landslides are less likely to occur on these slopes. The maximum LCI value (1.02) is a product of the low frequency of occurrence of the high SPI in the region coupled with the high frequency of occurrence of landslides in this area.

4.2. Correlation with lithology

The lithological unit exerts an important influence on landslide occurrence. There are nine lithological units on the geological map (Figure 6(l)). Among the different lithological classes, Class 5 covers the largest area: 257,377.3 km², amounting to ∼52.99% of the study area. However, only ∼19% of the landslides occurred in this area. A total of 65.15% of all landslides occurred in Classes 3 and 4, which covered only 26.51% of the study area (Figure 8(k)). Landslide concentration sharply fluctuates in all lithological classes because each class has distinct physical properties that assist or resist slope failures. The highest LCI with a value of 3.31 appeared in Class 3 lithologies, followed by Class 4 and 7 lithologies with LCI values of 2.05 and 1.73, respectively. The LCI values of the remaining lithology classes are <1, indicating that these classes are not prone to landslides.

4.3. Correlation with seismic factor

The curve in Figure 8(l) does not show a positive correlation between peak ground acceleration (PGA) values and the landslides, and the landslide frequency in classes of PGA ≤ 0.15 g is much higher than those in higher PGA classes. The total area of the PGA classes of ≤0.15 g covers 341,617.0 km², accounting for ∼70.32% of the total. The total number of landslides in this area is 29,776, amounting to 85.34% of the total landslides. The highest landslide frequency (52.07%) occurred in the PGA 0.05 g class, which covered only 21.46% of the study area. Overall, the values of frequency and LCI decrease with increasing PGA values. The minimum frequency (0.13%) occurs in the PGA 0.4 g class. In fact, the area of PGA 0.4 g only covers 1370.6 km² or ∼0.28% of the total study area. The area is too small to obtain objective statistical results. Sichuan Province is one of the regions with the highest frequency of earthquakes and has suffered some of the deadliest earthquakes (such as the 2008 CE Wenchuan earthquake, the 2017 CE Jiuzhaigou earthquake, and the 2022 CE Luding earthquake). The 2008 Wenchuan earthquake (Ms 8.0), which occurred in Sichuan Province, China, generated more than 60,000 landslides and caused many casualties. Seismic and its additional stress can result in the group occurrence of landslides (Figure 9), especially on slopes composed of sandstones and mudstones. The above results show that the occurrence of landslides does not increase with increasing PGA values, as expected. This indicates that the impact of PGA on landslide occurrence is perhaps controlled by other conditioning factors, such as elevation, lithology, and rainfall.

Figure 9. Examples of historical seismic landslides in the study area.

4.4. Correlation with rainfall

Figure 8(m) shows the landslide distribution in different rainfall classes. Landslide frequency is not a simple function of rainfall; for example, the 24-h rainfall class of >160 mm is not characterized by the highest landslide frequency. The area of the rainfall class of 30–40 mm is ∼145,979.5 km², accounting for 30.05% of the entire area. However, only ∼5.98% of the landslides occurred in this area. Landslides are preferentially distributed in the rainfall class of 90–140 mm, occupying 71.75% of the total landslide number. The number of landslides begins to increase with the increase in rainfall and gradually decreases after the rainfall exceeds 110 mm. This indicates that 110 mm of rainfall may be critical for the development of landslides in the region. Areas with rainfall of more than 110 mm have good physical and mechanical properties of slope components and high vegetation coverage, resulting in fewer landslides. Some studies have elucidated the general relationship between landslides and rainfall, and various mechanisms of rainfall-induced landslides have been proposed (Alvioli et al. 2014). Long-term rains will enhance surface seepage and lower the suction of the soil on the slope, which may lead to a considerable reduction in the strength of the slope and structural surfaces, thereby causing landslides (Figure 10). Nevertheless, copious rainwater will create favourable conditions for vegetation development, reducing the likelihood of landslides. It should be noted that short-term heavy rainfall generates transient surface currents and hinders surface seepage. Therefore, there is a threshold for the incidence of rainfall-induced landslides, which is uncertain. The LCI values generally fluctuate throughout the entire area. When the rainfall class is 130–140 mm, the LCI value reaches the maximum (3.72), followed by >160 (3.53), 140–150 (3.15), 120–130 (2.69), 150–160 (2.46), 110–120 (2.34), 100–110 (2.23), and 90–100 mm (1.82). Other rainfall ranges result in LCI values <1 and thus a lower probability of landslides occurring.

Figure 10. Photos of the rainfall-triggered landslides in the study area.

4.5. Correlation with anthropogenic factors

Landslides can be easily induced by unreasonable excavations on the toes of slopes. The engineering construction associated with road networks has greatly changed local landforms and the hydrogeological conditions of slopes, which have weakened slope stability. The relationship between landslide distribution and road density is shown in Figure 8(n). A large number of landslides are observed in the road density of ≤0.6 km/km², accounting for 77.81% of the total. The classes 0.6–0.8 and >0.8 km/km² generated 10.39 and 11.80% of the landslides, which covered 5.48 and 4.80% of the area, respectively. The landslide frequency abruptly decreases at higher road densities. This finding can be attributed to the small area of these higher road densities. The classes 0.6–0.8 and >0.8 km/km² exhibit the lowest landslide frequencies, but higher LCI values (1.90 and 2.46, respectively) can be observed. Overall, the LCI values increase with increasing road density and increase rapidly for road density >0.2 km/km². The minimum LCI value (0.32) occurs in the class <0.2 km/km², although this class covers the largest area.

5. Discussion

Identifying the factors controlling landslide occurrence is difficult because the relationship between landslides and the influencing factors varies spatially and temporally. However, a full understanding of these factors is crucial for elucidating potential risks in landslide regions. Therefore, GDM is used to detect the relative importance of control factors. Table 3 shows the q statistic of the possible influencing factors, which is based on the GDM. For the role of a single influencing factor, the q statistics for each factor range from 0.000 to 0.347. The q statistical values range from 0. 000 to 0.022 for five factors, including plan curvature, slope aspect, profile curvature, SPI, and surface roughness. This value range implied that these factors alone are related to the spatial distribution of landslides, but their effects are limited. Except for the above five factors, the q statistics of other single factors are >0.026. In Sichuan Province, lithology is identified as the dominant single factor (0.347), followed by rainfall (0.345), road density (0.204), elevation (0.198), and PGA (0.188).

Table 3. The q statistic of possible influencing factors.

	1	2	3	4	5	6	7	8	9	10	11	12	13	14
Q	0.347	0.000	0.345	0.026	0.000	0.063	0.085	0.014	0.188	0.204	0.022	0.000	0.071	0.198

1: lithology, 2: plan curvature, 3: rainfall, 4: slope angle, 5: slope aspect, 6: relief amplitude, 7: river density, 8: SPI, 9: PGA, 10: road density, 11: surface roughness, 12: profile curvature, 13: cutting depth, 14: elevation, Q: q statistic.

The occurrence of landslides is controlled by the interaction of internal factors and external factors. Different combinations of these geoenvironmental factors lead to differences in the spatial distribution of landslides. Table 4 shows the factor interaction effects on the landslide distribution. The explanatory powers of the factor interactions are much higher than those of individual factors. The interactions between rainfall and PGA have the greatest impact on the spatial distribution of landslides (q = 0.502), which can explain ∼50.2% of the spatial variation in landslide density. This indicates that in specific rainfall ranges (such as 100–110 and 110–120 mm), PGA is more closely related to landslide occurrence. For example, the 1933 Ms 7.5 Diexi earthquake caused tensile cracks in the mountain where the Xinmo landslide was located (Fan et al. 2017). Several major earthquakes (such as the 1976 Ms 7.2 Pingwu earthquake and the 2008 Ms 8.0 Wenchuan earthquake) since then, and the long-term effects of gravity and rainfall, finally led to the occurrence of catastrophic landslides (Figure 11). The interaction between rainfall and river density and road density was also significant, and the q statistic value was more than 0.47. Although the q values of an individual factor are small, the interaction between the above two factors and rainfall accounts for ∼47% of the spatial variation of the landslides. The interaction between lithology and rainfall, river density, PGA, road density, and elevation ranges from 0.417 to 0.499, which suggests that these interactions also play an important role in the spatial distribution of landslides in Sichuan Province. In addition, the interaction between river density and PGA can also explain over 40.7% of the variations in landslide densities, which is higher than all individual factors.

Figure 11. Overview of the 24 June 2017 Xinmo landslide.

Table 4. The interaction relationships between pairs of possible influencing factors.

	1	2	3	4	5	6	7	8	9	10	11	12	13	14
1	0.347
2	0.348	0.000
3	0.499	0.345	0.345
4	0.360	0.026	0.352	0.026
5	0.350	0.002	0.348	0.028	0.000
6	0.371	0.064	0.368	0.068	0.065	0.063
7	0.426	0.086	0.497	0.125	0.087	0.170	0.085
8	0.351	0.016	0.348	0.038	0.017	0.072	0.106	0.014
9	0.438	0.188	0.502	0.217	0.191	0.244	0.407	0.191	0.188
10	0.449	0.205	0.473	0.215	0.207	0.228	0.327	0.212	0.376	0.204
11	0.355	0.023	0.348	0.026	0.024	0.066	0.120	0.034	0.207	0.210	0.022
12	0.348	0.000	0.345	0.026	0.001	0.064	0.086	0.015	0.188	0.205	0.023	0.000
13	0.374	0.071	0.374	0.075	0.073	0.076	0.180	0.080	0.253	0.232	0.073	0.071	0.071
14	0.417	0.198	0.365	0.208	0.200	0.219	0.338	0.203	0.253	0.275	0.203	0.199	0.226	0.198

Bold font means relatively high enhancement.

Many landslides in Sichuan Province occurred in the Jurassic and Cretaceous strata that are composed of a hard rock lithology (sandstone) intercalated with a high percentage of soft and weak layers (e.g. mudstone). Mudstone belongs to the impermeable layer, and it is generally difficult for groundwater to enter. However, under the influence of tectonic stress or gravity, shear stress concentration is likely to occur at the contact surface of mudstone and sandstone, which leads to the destruction of the mudstone structure and the increase in pores at the contact surface. As a result, groundwater enters it and forms a muddy interlayer. The interlayer tends to swell and weaken under saturated conditions due to the high content of clay minerals, including montmorillonite, illite, and kaolinite, which show expansion, disintegration, and softening under periodic dry and wet conditions. The muddy interlayer has low shear strength and is heavily fractured, as observed in the field, leading to higher concentrated landslide activity on slopes consisting of these weak and fragmented rocks (Figure 12).

Figure 12. Longitudinal profiles of slope characteristics for the (a) Zhaoerping and (b) Kualiangzi landslides.

The continuous uplift of the Qinghai-Tibet Plateau led to crustal uplift and deep cutting of rivers in the study area. This process shaped the current alpine gorge geomorphology in the study area. During the process of downcutting by rivers, with the lateral constraint of the slope being removed, the slope is subjected to rebound deformation, and stress redistribution occurs at a certain depth of the slope (Huang 2008). As shown in Figure 13, the stress redistribution of the slope has similar characteristics to the stress distribution of the surrounding rock of the tunnel, including the stress-reduction zone (σ < σ₀), the stress-increase zone (σ > σ₀), and the initial rock stress zone (σ = σ₀) (Huang 2008). The stress redistribution of the slope is referred to as the ‘hump stress distribution’, and the stress-reduction and stress-increase zones (i.e. ‘hump zone’) correspond to the range of the unloading deformation of the slope. The stress-reduction zone is close to the bank slope of the river valley. Due to the stress release of the valley slope, the stress (σ) in this zone is less than the initial rock stress (σ₀). The stress in the stress-increase zone is higher than the initial rock stress due to the adjustment of the stress field. The initial rock stress region remains in its original state due to being unaffected by river downcutting. The stress at the valley bottom is mainly characterized by shear relaxation unloading and high geostress zones.

Figure 13. Stress distribution of valley slopes (modified after Huang 2008).

Compared with the stress-increase zone, the influence of the stress-reduction zone on the slope rock mass close to the river is more significant because both the strong unloading zone and the weak unloading zone often appear in this zone. The strong unloading zone (Figures 14(a,b)) is commonly related to the tension fracture zone and pressure-induced tension fracture zone, while the weak unloading zone (Figures 14(c,d)) is generally related to the shear fracture zone. In general, the damage evolution process of the slope close to the river can be divided into three stages: first epigenetic deformation, then time-dependent deformation, and finally failure (Figure 15). Epigenetic deformation is the deformation and rupture of the slope rock mass driven by lateral unloading and stress release during the process of river incision. The structural planes formed by epigenetic deformation include two types: vertical fractures (i.e. Type 1 unloading cracks) and horizontal fractures (i.e. Type 2 unloading cracks). Epigenetic deformation can release the slope stress, form the unloading relaxation zone, and provide favourable boundary conditions for subsequent deformation of the slope. When the epigenetic deformation is completed, the slope stress field will be transformed into a state dominated by the self-weight stress field. At this time, the slope will form a new stable structure and be in equilibrium. Of course, if there are adverse geological conditions, the slope would also undergo creep deformation (i.e. ‘time-dependent deformation’) under the action of the self-weight stress field. The continuous development of the time-dependent deformation is also the gradual formation of the potential slide surface. When the slide surface develops to a certain stage, the slope will enter the final failure stage.

Figure 14. Unloading zone exposed by exploratory adits.

Figure 15. Sketch of three-stage mode of high slope evolution (modified after Huang 2008).

Since the implementation of the Western Development Policy, especially in the past two decades, the number of lifeline projects (such as railways, highways, and oil and gas pipelines) has increased rapidly, and a large number of factories and energy mining industries have been developed. With the acceleration of urbanization, more people are settling in cities, and people are urged to improve their living conditions. Therefore, to retain arable land and obtain more usable land in the surrounding mountains and deeply cut gullies, human beings have to alter the land through excavation, landfills, and irrigation, which further damages the environment of the study area. The improper selection of location and construction methods can easily lead to landslides (Chen et al. 2022). A large number of railways, highways, and other infrastructures in the study area are being planned or have been built, such as the Sichuan–Tibet Railway, Baoji-Chengdu Railway, Chengdu-Kunming Railway, Sichuan-Gansu Highway, and Sichuan-Shaanxi Highway. During the construction of these projects, many artificial slopes (i.e. cut slopes) were formed in this area. The improper excavation of the slope toe disturbs the equilibrium of natural slopes, changes the local landforms and the hydrogeological structure of slopes, and accelerates the formation of landslides (Figure 16). Mining activities produce a large amount of waste slag, and a considerable concentration of overloading is acting on the slope. At the same time, the structure of these waste slags is relatively loose, and rainwater tends to infiltrate along the existing cracks to form a weak layer, which markedly destabilizes the slope and increases the frequency of landslides (Figure 17). Overall, the area of mining-induced landslides is larger than that of road-related landslides, and their regional distribution is slightly consistent with the characteristics of deep-seated landslides, which is logical because mining activities are related to the lithologic characteristics of the landscape.

Figure 16. Photos of the road-related landslides in the study area.

Figure 17. Photos of the mining-induced landslides in the study area.

6. Conclusions

Sichuan Province is a fragile area prone to landslides, which are controlled by complex geo-environmental factors. In this study, a total of 34,893 landslides collected from the Department of Natural Resources of Sichuan Province are used to analyze the relationships between influencing factors and the spatial distribution of landslides. Topographic factors, such as plan curvature (concave and convex areas), profile curvature (the classes < −0.001 and >0.001), slope angle (<24°), slope aspect (southeast, south, southwest, and west), relief amplitude (<222 m), cutting depth (<275 m), river density (>0.64 km/km²), surface roughness (<1.23), elevation (<2378 m), and SPI (<6) play important roles in landslide occurrences. In addition, other factors also sensitive to landslides in the study area include lithology (Classes 3 and 4), PGA (≤0.15 g), rainfall (90–140 mm), and road density (≤0.6 km/km²).

Geographical detector analysis reveals that the role of individual influencing factors is related to the spatial distribution of landslides, but their effects are limited. The factor interaction between rainfall and PGA, lithology, river density, and road density has a greater impact on the spatial distribution of landslides in Sichuan Province. The completeness of the landslide inventory influences the reliability of the spatial distribution analysis of landslides. Landslide mapping is vulnerable to uncertainties and inaccuracies that are dependent on surveyors’ knowledge backgrounds and skills. The uncertainty produced during the landslide mapping process propagates to the study area’s landslide distribution analysis. As a result, the analysis results of this study could be enhanced if supported by a more comprehensive inventory in a further study. Overall, the results of the present study can assist future hazard zonation, prediction, and prevention of landslides in Sichuan Province.

Additional information

Funding

This study was funded by the National Natural Science Foundation of China (52208359 and 52109125), the China Postdoctoral Science Foundation (2020M680583), and the National Postdoctoral Program for Innovative Talent of China (BX20200191).