Characteristics of Debris Flow Activities at Different Scales after the Disturbance of Strong Earthquakes—A Case Study of the Wenchuan Earthquake-Affected Area

Abstract

Of the catastrophic earthquakes over the past few decades, the 2008 Wenchuan earthquake triggered the greatest number of landslides and deposited a large amount of loose material on steep terrains and deep gullies, which was highly conducive to the occurrence of post-earthquake debris flows. It is of great importance to clarify the evolution of debris flow activity for hazard evaluation, prediction, and prevention after a strong earthquake, especially in the face of large debris flow hazards. We established a long-time span database consisting of 1668 debris flow events before and after the earthquake, with information including the occurrence time, location, and scale (small, medium, and large). In order to analyze how the environmental background before and after the earthquake controlled the debris flow activity, we examined various controlling factors, including the material source, topography (relative relief and slope degree), rainfall, normalized vegetation index, and lithology. After completing the analysis of the spatial and temporal evolution of the debris flow events in the database, a 10 × 10 km grid was introduced to grade the controlling factors in ArcGIS. Based on the same grid, the density of debris flow events for each scale in different time periods was calculated and graded. We introduced the certainty factor to figure out the spatial–temporal relationships between debris flow activities at each scale and the controlling factors. The results can provide guidance on how to dynamically adjust our strategies for debris flow prevention after a strong earthquake. Lastly, Spearman rank correlation analysis was performed to clarify the variation in the magnitude of the influence of controlling factors on the debris flow activities of different scales with time. This can provide a reference for the dynamic evaluation of debris flow hazards in the Wenchuan earthquake-affected area.

1. Introduction

Debris flows are water-laden masses of soil and debris that rush down mountainsides into channels, entraining objects along their paths and spilling onto valley floors to form lobed deposits [1]. This phenomenon is one of the most dangerous types of mass movement because of its characteristics of high velocity, high force of impact, and long runoff distance [2,3,4,5,6]. It is widely spread on relevant slope inclinations where grained materials are present [7]. A large amount of loose material, surface deformation, and cracks resulting from a strong earthquake can significantly increase surface material erosion and transport activities. Furthermore, more frequent and destructive debris flows occur in earthquake-affected areas, and rainfall thresholds for triggering debris flows are significantly lowered [8,9,10,11]. In recent decades, a number of strong earthquakes have occurred in different regions around the world. All of the post-earthquake secondary disasters have had serious impacts, such as the 1923 Kanto earthquake in Japan [12], the 1999 Chi-Chi earthquake in Taiwan [13,14], the 2005 Kashmir earthquake in Pakistan [15,16], the 2008 Wenchuan earthquake in China [17,18,19], the 2015 Gorkha earthquake in Nepal [20], and the 2016 Kaikoura earthquake in New Zealand [21]. These catastrophic earthquakes all generated a large number of loose deposits from co-seismic landslides, significantly increasing the activity of debris flows.

The 2008 Wenchuan earthquake (Mw 7.9) injured approximately 459,000 people, killed over 87,000 [22], and left at least 15 million people fleeing their homes and over 5 million homeless [23]. A large number of houses and infrastructures were damaged, with economic losses amounting to USD 150 billion [24]. Through high-resolution remote sensing interpretation, Xu, et al. [25] identified about 200,000 co-seismic landslides with a total area of approximately 1,160 km², the largest number of co-seismic landslides triggered by any global earthquake in the past few decades [11]. These co-seismic deposits were prone to being reactivated or even converted to debris flows by rainfall after the earthquake [26]. Several studies tracked the evolution of debris flows after the Wenchuan earthquake [11,27,28,29], but none of them reached the earthquake-affected area for more than 10 years of follow-up studies. Heavy rainfall on 20 August 2019 triggered another cluster of debris flows in Wenchuan, resulting in 26 people being lost, 12 deaths, the disruption of several major transportation routes, and economic losses of nearly USD 3.6 billion [30,31]. This incident shattered some of our previous understanding of post-earthquake debris flows and prompted the necessity to study the evolutionary activity over a longer horizon time [32].

In this study, we established a database of debris flow events for a long time before and after the earthquake, with a total of 1668 cases. The information includes occurrence time, location, and scale (small, medium, and large). After completing the statistical analysis of the database, a 10 × 10 km grid was introduced to grade the material source, topography, rainfall, normalized vegetation index, and lithology factors in ArcGIS. Based on the same grid, the density of debris flow events for each scale in different time periods was calculated and graded. Finally, the certainty factor and Spearman correlation coefficients were calculated to analyze the characteristics of debris flow activities for each scale under different controlling factors. We hope this study helps to elucidate the evolution of debris flow activity and the evolving impact of controlling factors in the Wenchuan earthquake-affected area, thus providing referable factual data for further studies.

2. Study Area

The study area is located in southwestern China, northeastern Sichuan Province. It extends along the Longmenshan fault zone, mainly involving Wenchuan County, Dujiangyan City, Mianzhu City, Beichuan County, Mao County, An County, Pingwu County, Qingchuan County, and other areas along the Yingxiu–Beichuan fault. The total area is about 31,630 km², and the specific geographical location is shown in Figure 1. The topography of the study area rises gradually from the Sichuan Basin in the southeast to the Hengduan Mountains in the northwest, with the lowest point in the Chengdu Plain at 477 m above msl, which is an alluvial plain and a flat dam in the intermountain valley. The highest point is in the Hengduan Mountains in the west at 6197 m above msl, which is a tectonic erosion alpine and extremely alpine landform. The southern mountainous region has an elevation difference of up to 4000 m over a distance of less than 50 km, making it one of the most treacherous regions in the world [33,34].

The study area is under the combined influence of a humid subtropical monsoon climate and a warm temperate continental semi-arid monsoon climate, superimposed on huge topographic differences, which make rainfall vary spatially, showing a gradual decrease from southeast to northwest with an annual average rainfall of about 800 mm for the entire area [35]. Under the influence of the monsoon, rainfall is mainly concentrated in May–September, accounting for 78% of the annual rainfall, with July and August being the most concentrated. There are many rivers in the study area that belong to three major water systems from north to south, namely the Jialing River system, the Tuo River system, and the Min River system, all flowing from the northwest to the southeast and finally merging into the Yangtze River. The Min River is the largest river in the study area, with a maximum runoff of 2340 m³/s and a maximum flow velocity of 6.9 m/s [36].

A variety of lithological types, mainly metamorphic sandstones, slates, magmatic rocks, and carbonates, with others being mudstones, sandstones, shales, phyllite, and quaternary deposits, were exposed in the study area [37]. The conflict between the Indian Plate and Eurasian Plate is the origin of the shaping of the topographic and geomorphic features of the region, and along with the uplift of the Qinghai–Tibet Plateau, the Longmenshan fracture zone was formed. It is generally believed that the Longmenshan fault originated in the Mesozoic and is located on a steeply northwest-dipping plane with an inclination of 60~70 degrees, which is a retrograde fault [10]. On 12 May 2008, the Wenchuan 7.9-magnitude earthquake occurred near Yingxiu in Wenchuan County. The epicenter was about 14 km from the surface and extended in a unilateral rupture toward the northeast to the vicinity of Qingchuan, forming a surface rupture zone of about 240 km. The Wenchuan earthquake had a mainshock that lasted about 120 s, and it was the worst natural disaster in China following the 1976 Tangshan earthquake. This earthquake caused severe damage not only in Sichuan but also in Gansu, Shanxi, and Ningxia provinces [38].

3. Data and Method

3.1. Debris Flow Events

In order to study the long-term characteristics of the debris flow activity in the earthquake-affected area, we need to collect as many debris flow events as possible. Usually, debris flow events are collected through the following methods: (i) field surveys; (ii) literature searches; (iii) government department records; (iv) existing geological hazard databases; (v) remote sensing images; and (vi) news reports. Through a long period of field surveys and the literature collection, supplemented by disaster and danger records from the Natural Resources Department and other debris flow event databases combined with remote sensing images and news reports as verification, a comprehensive database of debris flow events before and after the Wenchuan earthquake in the study area was established, including the time of occurrence of the event, Global Positioning System (GPS) coordinates of the gully, and the event scale level. As shown in

Table 1. Data sources of the debris flow events.

Source	Time Span	Quantity	Total Number of Events
			Time of Occurrence	GPS Coordinates	Scale Level
Field surveys and the literature	2008–2020	568	1668
Natural Resources Department	2000–2020	691
Geological Disaster Database of the Southwest Topographic Abrupt Change Zone	1952–2009	409

, a total of 568 post-earthquake events were accumulated from field surveys and the literature collection up to 2020. Then, they were supplemented with 691 debris flow events collected and screened by the Natural Resources Department from 2000 to 2020. It should be noted that since the reporting mechanism of geological hazards in Sichuan Province was gradually established and improved after the Wenchuan earthquake, some debris flow events were missing before 2010. Therefore, 409 debris flow events during the period of 1952–2009 were also supplemented and screened from the geological disaster database of the Southwest Topographic Abrupt Change Zone. Finally, we established a database of debris flow events in the study area from 1952 to 2020, with a total of 1668 events. All debris flow events were triggered by rainfall.

The scale level corresponding to all events was classified by the investigators after the field investigation of debris flow events or visits to historical debris flow events according to the criteria (

Table 2. Classification of debris flow scale level (extracted from C.4-DZ/T 0261-2014).

Classification Indicator	Level	Criteria (×104 m3)
Volume of event washout (V)	Small	V < 2
	Medium	2 ≤ V < 20
	Large	V ≥ 20

) in the “DZ/T 0261-2014 Landslide Collapse Debris Flow Disaster Investigation Specification” issued by the Ministry of Land and Resources. Cases of different levels of events are shown in Figure 2.

3.2. Controlling Factors

We examined four factors controlling debris flow occurrences in the study area (

Table 3. Information of controlling factors.

Controlling Factors	Resolution	Time Coverage	Source	Factor Grading Criteria Based on a 10 × 10 km Grid
Landslide area percentage (LAP)	1 × 1 km	/	Xu et al., 2014	1. 0 2. 0–0.05 3. 0.05–0.1 4. 0.1–0.15 5. >0.15
Relative relief (RR)	30 × 30 m	/	ALOS	1. <300 m 2. 300–400 m 3. 400–500 m 4. 500–600 m 5. >600 m
Slope degree (SD)	30 × 30 m	/		1. <5° 2. 5–10° 3. 10–15° 4. 15–20° 5. >20°
Rainfall (RAIN)	25 × 25 km	1952–2019	Global Precipitation Climatology Centre dataset (GPCC)	1. 500–600 mm 2. 600–700 mm 3. 700–800 mm 4. 800–900 mm 5. 900–1000 mm 6. 1000–1100 mm 7. 1100–1200 mm 8. >1200 mm
Normalized vegetation index(NDVI)	1 × 1 km	1998–2019	Research and Environment Science and Data Center	1. <0.7 2. 0.7–0.75 3. 0.75–0.8 4. 0.8–0.85 5. >0.85
Lithology (LITH)	1:500,000	/	China Geological Survey (CGS)	1. Quaternary Sediments (QS): gravel, sand, and clay. 2. Metamorphic Rocks (MR): fragmentized phyllite and metamorphosed slate. 3. Flysch Deposits (FD): marl, sandy shale, mud interbedded with greywacke, coarse and fine sandstone, and sandy conglomerates. 4. Carbonate Rocks (CR): limestone, dolomite, and marl 5. Volcanic Rocks (VR): granite, andesite, basalt, and pyroclastic deposits.

Note: Bold represents dynamic controlling factors.

). (1) Material source factor: based on the co-seismic landslide database established by Xu, et al. [25], we used the landslide area percentage (LAP) to depict the spatial distribution of loose material after the Wenchuan earthquake. LAP was defined as the percentage of the area affected by landslides. (2) Topographic factors: based on the DEM (30 × 30 m) data downloaded from the ALOS PALSAR RTC dataset, ArcGIS was used to calculate the relative relief (RR) and slope degree (SD) of the study area to portray the spatial distribution differences of topographic conditions. (3) Water source factor: based on the monthly rainfall data (25 × 25 km) collected and extracted from the Global Precipitation Climatology Centre dataset (GPCC) from 1952 to 2019, the precipitation in each annual monsoon season (May–October) was calculated to portray the spatial and temporal variation in rainfall in the study area. (4) Other factors: normalized vegetation index (NDVI) and lithology (LITH). The Resource and Environment Science and Data Center of the Chinese Academy of Sciences (CAS) only provided annual NDVI data (1 × 1 km) from 1998 to 2019. These data were used to characterize the dynamic changes in the regional vegetation cover. We also collected a geological map on a scale of 1:500,000 from the Geological Survey of China (GSC) to provide information on the lithology of the study area. It is widely recognized that lithology plays an important role in determining landslide hazards because different geological units have different susceptibilities to landslide occurrences, even when the slope failure is not triggered by an earthquake [25].

As a compromise between preserving the spatial fidelity of the debris flow catchment and detecting the large-scale spatial trends of the controlling factors in the entire study area, we generated 371 10 × 10 km grid cells in the study area using ArcGIS. Then, we calculated the mean values of each factor, including LAP, RAIN, RR, SD, and NDVI, in each grid cell and graded the factors by the criteria in Table 3. For LITH, we divided the rock types in the study area into five major categories (Table 3) and gridded them using the same 10 × 10 km grid. Finally, Figure 3 shows the grading results of the factors after gridding. As RAIN and NDVI were used as dynamic factors, the grading results obtained from the average values of 2011–2013 were used as representative displays.

3.3. Certainty Factor

The certainty factor (CF) is a probability analysis method proposed and improved by Shortliffe and Buchanan [39] and Heckerman [40]. It can quantitatively calculate the sensitivity of factors affecting the occurrence of events. Many studies have used certainty factor models to establish quantitative relationships between landslide activities and controlling factors, allowing them to study the long-term evolution characteristics of landslides [41,42]. Based on the quantitative density of debris flow events under each grade of factor, we used the following formulas to calculate the certainty factor (CF) of debris flow activities controlled by the factors:

$$P_s=N/A$$

(1)

$$P_i=N_i/A_i\left(i=1, 2\dots 5\right)$$

(2)

$$CF=\{\begin{array}{c}\frac{P_i-P_s}{P_i\left(1-P_s\right)},P_i\ge P_s\\ \frac{P_i-P_s}{P_s\left(1-P_i\right)},P_i<P_s\end{array}$$

(3)

where N is the total number of debris flow events in the study area; A is the area of the study area (km²); $P_s$ is the quantity density of debris flow in the study area (/km²); $N_i$ is the number of debris flow events in each factor grade; $A_i$ is the area of each grade of the controlling factor (km²); $P_i$ is the quantity density of debris flow in each grade of the controlling factor (/km²); and CF is the certainty factor of the grade corresponding to the controlling factors, and the value range is [−1,1]. A positive value indicates an increase in the certainty of debris flow activity. In contrast, a negative value corresponds to a decrease in the certainty of debris flow activity. A CF value close to 0 means that the prior probability is similar to the conditional probability. It is difficult to give any indication about the certainty of the debris flow activity.

3.4. Spearman’s Rank Correlation Coefficient

The Spearman rank correlation coefficient is a nonparametric statistical method that uses the rank magnitude of two variables for correlation analysis [43,44]. It operates on data ranks rather than the original data, has no requirement for the distribution of the original variables, and has a wider application than the Pearson correlation coefficient [45]. The idea of the rank correlation coefficient is simple. First, each variable is sorted and graded from low to high. Then, the grade difference of each data pair under two variables is recorded, and the calculation is performed. In this study, the debris flow densities of each grid cell in different time periods were graded equidistantly at 0.1/km² intervals and then correlated with the graded results of the controlling factors (Table 3). Lastly, Spearman rank correlation analysis was performed based on the following formulas:

$$D_i=\left(U_i-V_i\right)$$

(4)

$$R_s=1-\frac{6{{\displaystyle \sum }}_{i=1}^n{D}_i{}^2}{n\left(n^2-1\right)}$$

(5)

where $U_i$ and $V_i$ are the rankings of the two variables in order of size or superiority, respectively; $D_i$ is the difference in grading between the two variables; N is the number of data pairs; $R_s$ is the Spearman correlation coefficient of the two variables. The value interval of $R_s$ is [−1,1]. When $R_s$ > 0, there is a positive rank correlation between the two variables; when $R_s$ < 0, there is a negative rank correlation between the two variables; when $R_s$ = 1, this indicates that the two variables have exactly the same rank and there is a perfectly positive correlation; when $R_s$ = −1, this indicates that the two variables have exactly opposite ranks and there is a perfectly negative correlation; when $R_s$ = 0, this indicates that the two variables are not correlated.

4. Results

4.1. Evolution of Debris Flow Activity

Table 4. Statistics of the debris flow events over time.

Temporal Interval	Number of Debris Flow Events				Major Debris Flows
	Subtotal	Small	Medium	Large	Date	Affected Districts
1952–2007	98	40	47	11
2008	323	97	159	67	9.24	Beichuan
2009	68	34	29	5
2010	97	73	11	13	8.13	Mianzhu/Dujiangyan/Wenchuan
2011	167	137	23	7	7.3/7.28/8.20	Wenchuan/Maoxian/Anxian/Lixian
2012	207	146	52	9	8.17	Pengzhou/Shifang/Mianzhu/Anxian/Beichuan
2013	287	130	104	53	7.10	Pengzhou/Anxian/Mianzhu/Beichuan/Maoxian/Wenchuan
2014	34	29	4	1	7.9	Lixian/Wenchuan
2015	3	3	0	0
2016	2	1	0	1
2017	15	14	1	0
2018	133	90	39	4	6.25/7.10	Beichuan/Pingwu
2019	53	38	11	4	8.20	Wenchuan
2020	181	122	49	10	8.17	Beichuan/Pingwu/Mianzhu/Wenchuan
Total	1668	954	529	185

and Figure 4 show the statistics and presentation of the debris flow events in the database. There were 98 events recorded before the earthquake (1952–2007), and 85% were small and medium. After the earthquake (2008–2020), there were 1570 records, with small and medium events accounting for 89%. The debris flow activity in the study area increased significantly after the earthquake, with the strongest years being 2008 and 2013. The rainfall triggered several debris flows during this period, including the 9.24 Beichuan debris flow in 2008, the 8.13 Yingxiu, Longchi, and Mianzhu debris flows in 2010, and the 7.10 Wenchuan debris flow in 2013. We found that a large number of debris flow events occurred in 2008 when rainfall was lower compared to 2010 and 2013. The distribution of these debris flow events was similar to the distribution of co-seismic landslides (see Figure 3, LAP). It is observed that loose material induced by the earthquake can be promptly initiated by rainfall to form debris flows. From 2014 to 2017, the debris flow activity in the study area diminished, with only one small cluster of debris flows occurring in 2014. In 2018, the debris flow activity intensified again, and debris flows still occurred in areas such as Wenchuan and Beichuan under the excitation of rainfall. Still, the number and scale of debris flows were attenuated compared to 2008–2013. According to the database, more than 85% of the debris flows in the earthquake area occurred during major debris flow activities (a cluster of debris flow events triggered by a particular rainfall process), and this percentage reached 95% for large events.

According to Figure 5a, the larger the debris flow event, the more concentrated it is in areas with a high density of co-seismic landslides (see Figure 3, LAP), because such debris flows require a greater supply of material sources. Under the influence of rainfall, loose material was continuously consumed, and debris flow activity showed a fluctuating downward trend in both number and scale. However, it still had not returned to pre-earthquake levels by 2020 (Figure 5b).

4.2. Relationship between Debris Flow Activity and Controlling Factors

Since events in 1952–2007 and 2014–2017 were few, we portray the events in these two intervals as a single entity, and the debris flow activity with time clearly shows four different periods: 1952–2007, 2008–2013, 2014–2017, and 2018–2020 (Figure 5b). To better reveal the characteristics of the activity in the first 6 years after the earthquake, we equated 2008–2013, thus dividing the debris flow activity into five periods: period Ⅰ (1952–2007), period Ⅱ (2008–2010), period Ⅲ (2011–2013), period Ⅳ (2014–2017), and period Ⅴ (2018–2020). Then, we assigned grade values for each controlling factor to debris flow events in each period and calculated the certainty factor (CF) of the debris flow events on three scales over time (Figure 6). It should be noted that, due to the missing 2020 data for the RAIN factor and the missing 1952–1997 and 2020 data for the NDVI factor, the debris flow events during the missing time periods for the two types of controlling factors are not included in the calculation and analysis.

Small and medium debris flow activities in the study area increased with the relative relief and slope degree prior to the earthquake. On the other hand, large debris flows were more likely to occur in areas with a gentle relative relief and slope. This may be due to the fact that the material sources of the pre-earthquake debris flow were not abundant, and the greater slope degree and relative relief were not conducive to the accumulation of material sources for large events in the gully. After the earthquake, debris flow activities of different scales in the study area were mainly concentrated in the areas with LAP > 0.05. Among them, the relative relief and slope of the small and medium debris flow-prone areas decreased significantly, from RR > 500 m and SD > 20° to RR: 300–400 m and SD: 10–15°. After 2014, due to the decrease in landslide activity and sediment supply capacity, debris flow-prone areas started to move to areas with a greater relative relief and slope (RR: > 600 m, SD: > 20°), with the most significant trend being large debris flows (note the increase in the RR and SD curves in the above-mentioned areas). In 2018 (10 years after the strong earthquake), the active areas of small and medium debris flows began to spread to areas with a smaller LAP while, at the same time, the main activity of large debris flows was still dominated by the LAP > 0.1 area, where co-seismic landslides were abundant. It can be seen that debris flows were mainly concentrated in areas rich in co-seismic landslides after the earthquake, the requirements of debris flow for topography were reduced compared with those before the earthquake, and the larger debris flow events had higher requirements for topography (relative relief and slope degree) over time.

In terms of rainfall, the stronger the rainfall before the earthquake, the more active the debris flow. In the post-earthquake period of 2008–2010, the rainfall in the study area did not increase significantly compared with the pre-earthquake period, but the concentration of medium and large debris flow activities at this time was not in the area with greater rainfall but concentrated in the area with 700–800 mm of rainfall. At that time, less rainfall could initiate material sources to form larger debris flows. In 2011–2013, the rainfall in the study area increased significantly, and the debris flow activities were concentrated in the area with >800 mm of rainfall. From 2014 to 2017, the rainfall in the study area weakened, and the debris flow-prone area coincided with the most abundant rainfall area (800–900 mm). After entering 2018, the rainfall in the study area increased again, and the small debris flow activity in the >1200 mm region was significantly and substantially attenuated, but large debris flow events continued to occur. It can be seen that the requirement of rainfall for debris flow events in the study area decreased after the earthquake compared to the pre-earthquake period and then gradually increased, especially for large debris flows, where this trend was most pronounced.

Combining the results of the analysis of LAP, topography, and rainfall after 2018 with the massive depletion of material sources in the dense LAP area, there were no longer enough material sources in small and medium gullies to participate in debris flow activities. At this time, only catchment basins with extremely abundant material source reserves and steeper topographic conditions could initiate the remaining material sources to form large debris flows under the effect of heavy rainfall. Therefore, this stage of the large debris flow-prone area had higher requirements for the three major conditions of material sources, topography, and rainfall.

The distribution of debris flow activity under the NDVI before the earthquake did not show obvious characteristics, while the Wenchuan earthquake caused different degrees of damage to the vegetation and the NDVI in areas with high debris flow activity was low (0.7–0.8). After 2014, with the recovery of vegetation, the NDVI in small and medium debris flow-prone areas increased rapidly, but the recovery of vegetation in large debris flow-prone areas was slow. It is generally believed that the recovery of vegetation can play a role in the consolidation of loose materials. However, in large debris flow-prone areas with greater rainfall and topographic conditions, even nearly 10 years after the earthquake, the channel runoff formed by strong rainfall can still activate the deposits, which will make it difficult for the recovery rate of the NDVI to catch up with that of the other areas. In terms of lithology, the debris flows were prone in the flysch deposits (FD) before the earthquake. After the earthquake (2008–2013), small and medium events were still most active in the flysch deposits (FD), and large debris flows were more active in the carbonate rocks (CR) and volcanic rocks (VR). After 2014, the active area of small and medium debris flows gradually shifted to the metamorphic rocks (MR), while large debris flows were more concentrated in the volcanic rocks (VR). It is observed that large debris flows after the earthquake were more concentrated in areas with harder rock properties, while small and medium debris flows tended to be in softer rock areas. The Wenchuan earthquake caused more damage to the hard rock areas. On the other hand, the topography of hard rock areas is steeper, which is also conducive to the development of large debris flows after the earthquake.

4.3. Evolving Controls

To analyze the importance of LAP, topography (relative relief and slope degree), and rainfall in the debris flow activity of each scale in the study area over time, a Spearman rank correlation analysis was performed, and the results are shown in

Table 5. Changes in correlation coefficients between debris flow activity and the main controlling factors at different scales.

Level	Period (Ⅰ–Ⅴ)	Spearman Correlation Coefficient (RS)
		LAP	Rainfall	Relative Relief	Slope Degree
Small	Ⅰ: 1952–2007	/	0.20	0.12	0.10
	Ⅱ: 2008–2010	0.34	0.19	0.07	0.05
	Ⅲ: 2011–2013	0.34	0.30	0.07	0.07
	Ⅳ: 2014–2017	0.16	0.15	0.21	0.21
	Ⅴ: 2018–2019 (2020)	0.19	−0.05	0.00	0.00
Medium	Ⅰ: 1952–2007	/	0.15	0.17	0.14
	Ⅱ: 2008–2010	0.31	0.02	−0.16	−0.13
	Ⅲ: 2011–2013	0.34	0.27	−0.01	0.00
	Ⅳ: 2014–2017	0.11	0.11	0.08	0.08
	Ⅴ: 2018–2019 (2020)	0.07	−0.09	0.00	0.01
Large	Ⅰ: 1952–2007	/	0.06	−0.05	−0.06
	Ⅱ: 2008–2010	0.36	0.09	−0.17	−0.15
	Ⅲ: 2011–2013	0.37	0.22	0.04	0.05
	Ⅳ: 2014–2017	0.07	0.06	0.10	0.08
	Ⅴ: 2018–2019 (2020)	0.22	0.18	0.14	0.17

Note: Red: significant correlation at the 0.01 level; orange: significant correlation at the 0.05 level; and black: insignificant correlation.

. The significant correlation factors were extracted and ranked, and the results for each scale on the timeline are shown in Figure 7.

According to Table 5 and Figure 7, before the earthquake, small and medium debris flows in the study area were mainly controlled by hydrodynamic conditions, which were determined by a combination of rainfall and topography. However, as gentler topographies are more conducive to accumulating material resources, the control of hydrodynamic conditions in large debris flows was relatively weakened. During the post-earthquake period, i.e., 2008–2013, debris flow activity was most closely related to LAP, with correlation coefficients higher than 0.3 for different scales. Rainfall ranked after LAP in importance during this period, while topography did not influence debris flow activity for most of the time. Interestingly, topography even had a negative effect on debris flow activity in the first 3 years after the earthquake. This particular phenomenon should be attributed to the fact that huge amounts of loose materials could participate in debris flow activity, and many gullies with poor topographic conditions were able to form medium or even large debris flows. After 2014, the importance of LAP decreased significantly, and the correlation coefficient dropped to about 0.2 when the importance of topography began to emerge. We found that after 2018, 10 years after the earthquake, large debris flows were controlled by a combination of LAP, rainfall, and topography, and their requirements for hydrodynamic conditions increased significantly compared to the time just after the earthquake (note the increase in R_S values from Ⅱ to Ⅴ for large debris flows in Table 5).

5. Discussion

Based on the characteristics of debris flow activity, we delineated five typical subregions from the study area, but the alpine and extremely alpine mountain areas above the snow line were not included (Figure 1). Subregions 1, 2, and 4 are located in the south and belong to the Min River and Tuo River systems, while the rest of the subregions are located in the north and belong to the Jialing River system. The LAP, topography, and rainfall of each subregion as well as the debris flow activity statistics are shown in

Table 6. Statistical table of background factors and post-earthquake events in each subregion.

ID	Restricts Involved	LAP	Relative Relief (m)	Slope Degree (°)	Rainfall (mm)	Event Density (/km2)	Percentage of Large Events
1	Anxian/Mianzhu/Shifang/Pengzhou/Dujiangyan	13.81%	419.67	14.53	936.85	0.19	16%
2	Wenchuan/Maoxian	11.53%	582.47	20.32	1021.23	0.05	18%
3	Old Beichuan/East Pingwu	8.14%	394.13	13.86	751.23	0.24	8%
4	Lixian	0.99%	586.84	20.36	887.66	0.09	2%
5	West Pingwu	0.26%	391.77	13.93	753.3	0.05	1%

.

Subregions 1 and 3, which are distributed along surface ruptures with abundant sources of material, were extremely active in the first 6 years after the earthquake when debris flows were less demanding in terms of topography and rainfall. This phenomenon is similar to the changing conditions of post-earthquake debris flows in the Chi-Chi earthquake-disturbed area in Taiwan. In the Chenyulan Basin, Lin, et al. [4] found that the maximum hourly rainfall intensity and critical accumulated precipitation necessary to initiate debris flow were reduced to as low as 1/3 of the pre-earthquake values. Meanwhile, after the earthquake, debris flows were observed even in gullies with less topographic relief and an effective drainage area smaller than 0.03 km². A strong earthquake can transform many non-debris-flow gullies into debris-flow gullies. Back to the case of the Wenchuan earthquake, Huang and Li [46] collected geological hazard statistics released by the Ministry of Land and Resources and pointed out that a large number of newly formed debris flows were distributed along the surface ruptures. The Wenchuan earthquake also transformed many gullies into debris flows soon after. As time passed, the results of this study in the long-term observations show that the intensity of its activity decreased significantly after 2014, with small and medium events prevailing.

Compared with subregions 1 and 3, subregion 2 was also strongly affected by the earthquake. However, the topography is steeper, and the catchment area of the gullies on both sides of the river is generally larger, meaning that large debris flow events can still be triggered by heavy rainfall 10 years after the earthquake. Chen, et al. [47] used NDVI data to quantify the recovery rate of co-seismic landslides in the key areas in subregions 1, 2, and 3 as 0.033/year, 0.024/year, and 0.037/year, respectively. The recovery rate of co-seismic landslides within subregion 2 was slower, and the supply of material sources of debris flow was more sustainable. However, not all gullies within subregion 2 could experience debris flow continuously. After tracking the co-seismic landslides, Fan, et al. [26] and Xiong, et al. [48] found that co-seismic landslides located in larger catchment basins remained active for longer, and debris flow activities in the corresponding basins were sustained within subregion 2.

Subregion 4 is located in the mid-alpine valley, which is the steepest in the study area. Debris flows were more frequent than in other regions before the earthquake, but the density of co-seismic landslides after the earthquake was much lower than that of the subregions near the surface ruptures, and debris flow events were still mainly small and medium. Subregion 5 was the least affected by the Wenchuan earthquake and had the weakest intensity of post-earthquake debris flow activity.

In summary, the reasons subregion 2 became a hotspot with a frequent occurrence of large debris flows after the earthquake are as follows. Firstly, this subregion was strongly affected by the earthquake, resulting in a large amount of loose material, which is the prerequisite condition for the generation of large debris flows after strong earthquakes. Secondly, the steep topographic conditions were conducive to runoff convergence and erosion, which was key to the sustained occurrence of large debris flows long after the earthquake. Finally, rainfall was the final touch needed to pull the trigger.

6. Conclusions

Long-term monitoring of debris flows in the Wenchuan earthquake-affected area can not only support the study of the evolution of debris flow activity but also be of significant importance to the evaluation, prediction, and prevention of debris flow hazards in the disturbed area. In this study, we established a database of debris flow events over a long period before and after the earthquake, with a total of 1668 cases, where the information included occurrence time, location, and scale. Then, we statistically analyzed the database and used the certainty factor to characterize the evolution of different scale events over time based on the controlling factors (LAP, rainfall, relative relief, slope degree, NDVI, and lithology), which were graded using a 10×10 km grid. Based on the same grid, we graded the density of events on each scale and performed a Spearman correlation analysis with LAP, rainfall, and topography factors. The conclusions are as follows:

(1)

Debris flow activity intensified rapidly after the Wenchuan earthquake, followed by a fluctuating downward trend in the number and scale of events. As of 2020, debris flow activity in the earthquake-affected area had not returned to the pre-earthquake level;

(2)

Post-earthquake debris flows were mainly concentrated in areas with abundant co-seismic landslides, and the requirements for topography and rainfall were reduced compared to those before the earthquake. Over time, the larger the debris flow events, the higher their requirements for topography and rainfall. The large debris flows became more concentrated in volcanic rock areas with time, and vegetation recovery in these areas was slower;

(3)

Pre-earthquake debris flows were mainly controlled by rainfall and topography. From 2008 to 2013, the debris flow activity was most closely related to the co-seismic landslides, followed by rainfall, while topography was not important. After 2014, the importance of the co-seismic landslides diminished significantly, and the importance of topography became apparent at this time;

(4)

The Wenchuan district has been a hotspot for large debris flows that have continued to occur after the earthquake due to the abundance of co-seismic landslides, superior topographic conditions, and sufficient rainfall during the monsoon season in this region.