Variation in Debris-Flow-Prone Areas with Ecosystem Stability: A Case Study of the Qipan Catchment in the Wenchuan Earthquake Region

Abstract

The spatial distribution of vegetation in a basin has a far-reaching influence on the potential for sediment separation and transport capacity. However, many landslides induced by strong earthquakes have greatly changed the existing pattern, which further increases the probability of debris flow in a basin during heavy rainfall and has a significant impact on the stability of the basin. Thus, this study selected the debris flow basin in the Qipan catchment of the Wenchuan earthquake area as the research object. Multisource and high-precision remote sensing images were used to analyze the land use changes in the basin, and the index of connectivity (IC) was introduced to analyze the evolution of sediment transport capacity. An ecosystem stability assessment method suitable for post-earthquake debris flow basins was proposed. Through quantitative assessment of the ecosystem stability of the basin after the Wenchuan earthquake in 2008 and the two debris flow events after the earthquake, the dynamic relationship between the debris-flow-prone area and the ecosystem stability of the basin was revealed. The results showed that the stability of the ecosystem in the Qipan catchment increased annually, indicating a stable and substable state. The spatial distribution characteristics were lower in the north and south and greater in the middle. By comparing the evaluation results with the actual terrain change trend, the accuracy and feasibility of the evaluation method are verified. The results of this study provide a scientific basis for the formulation of regional disaster prevention strategies and help to accelerate the improvement of regional stability in debris-flow-prone areas.

1. Introduction

Ecosystems play a vital role in maintaining regional ecological security and health management and are also a core factor in regulating sediment transport, soil stability, and groundwater flow [1,2]. However, the intensification of global environmental pressure, including abnormal climate changes (such as global warming, extreme drought, and rainstorms) and enhanced human activities, has caused severe external disturbances to the global ecosystem and has had a profound impact on the function and structure of the ecosystem, which seriously threatens its stability [3,4]. In the past decade, global warming has significantly increased the frequency of extreme weather events, such as rainstorm events caused by extreme precipitation, resulting in the frequent occurrence of debris flow disasters in China [5,6]. Therefore, systematically monitoring and evaluating the stability of regional ecosystems is crucial for developing effective strategies to protect ecological security and resist the adverse effects of climate change.

In the field of ecology, ecosystem stability is a multidimensional concept defined as the ability of ecosystems to resist environmental disturbances in their structure and function [7,8,9]. Donohue et al. [10] proposed that environmental disturbances not only increase the complexity of stability but also change the interaction between the components of stability in many ways. Therefore, the assessment of ecosystem stability cannot be simply summarized as a single measurement standard; it must be analyzed as a multidimensional structure, and for different evaluation objects, the selected indicators need to be adjusted according to the characteristics of their respective systems [11,12]. In addition, the vast majority of ecosystem stability assessment objects are forests, lakes, wetlands, arid deserts, or alpine regions [13,14]. In contrast, relatively few studies have investigated the ecological stability of postdisaster areas, especially debris-flow-prone areas.

The stability of ecosystems in montane landscapes is increasingly threatened by natural hazards (e.g., debris flows and landslides), especially in regions with strong earthquakes [15]. Strong earthquakes generally induce large-scale landslides that remove vegetation from the surface and transport it to sedimentary zones [16,17]; under heavy rainfall conditions, these landslides are prone to transform into debris flows, resulting in the destruction of ecosystems [18,19]. Examples include the 1999 Chi-Chi earthquake [20], 2008 Wenchuan earthquake (M_w 7.9) [21], and 2017 Jiuzhaigou earthquake (M_w 6.5) [22], which significantly accelerated the decline in regional ecosystem stability. In addition, a decrease in ecosystem stability increases the probability of debris flow, which leads to a further decline in ecosystem stability and the formation of a vicious cycle [23,24]. Therefore, based on the close relationship between debris flow activity and ecosystem stability, it is very important to evaluate the stability of debris-flow-prone areas.

With the deepening of research on ecosystem stability assessments, research methods have changed from qualitative to quantitative evaluation [25,26]. The development of 3S technology has promoted the wide application of remote sensing images in the evaluation of regional ecosystem stability. Data such as the normalized difference vegetation index (NDVI) and enhanced vegetation index (EVI) provide new means for evaluating region stability [27]. In debris-flow-prone areas, vegetation can stabilize slopes and effectively reduce sediment accumulation, which plays an important role in restricting the development of debris flows [28]. Therefore, many scholars have studied the regional stability of vegetation restoration after earthquakes and secondary disasters. For example, the vegetation restoration process in the Wenchuan earthquake area was predicted using the long-term sequence change in the NDVI as an indicator [29,30], but these studies ignored the internal factors influencing debris flow occurrence. If the characteristic factors of debris flow areas are included in regional stability assessments, the judgment and evaluation of debris flow susceptibility can be significantly improved. Among them, sediment connectivity, as a representation of sediment transport capacity, has been used by many scholars to evaluate and understand the transport efficiency of sediment from sources to sedimentary areas in basin systems [31,32,33]. This change is directly related to the probability of a debris flow, which is helpful for reflecting the impact of a debris flow on an ecosystem.

By comprehensively considering the factors related to vegetation and debris flows, we can understand the ecological stability of debris-flow-prone areas more deeply and identify potential risks and vulnerabilities more accurately. Therefore, we take the Qipan catchment, a typical debris-flow-prone area, as an example and use sediment connectivity theory to construct an ecological stability assessment method suitable for watershed-scale, debris-flow-prone areas, explore the ecological effects of vegetation and sediment connectivity, and reveal the relationships between the internal dynamics and ecological decay laws of landslides after earthquakes to aid in disaster prevention and mitigation.

2. Study Area

The Qipan catchment, located in the Aba Tibetan Autonomous Region, Sichuan Province, approximately 5 km away from Wenchuan County, is a debris-flow-prone area. It is the first tributary of the Minjiang River [33]. The Minjiang River has 15 branches, with a mainstream length of 15.8 km and a drainage area of 54.2 km². The topography is mainly composed of rugged mountains and deeply incised valleys. Most of the slopes are steep, ranging from 40° to 45°, and the altitudes are approximately 1320–4372 m, as shown in Figure 1. Precipitation is abundant, mostly concentrated from May to September every year, and heavy rain or heavy rainstorm events occur in some areas. This abundant precipitation provides favorable conditions for the reactivation of earthquake-induced geological disasters, the erosion of loose accumulations, and the restoration of vegetation growth.

The vegetation coverage in this area is as high as 60%, and the diversity is high. It can be divided into five main vertical vegetation zones with increasing altitude: subtropical, evergreen, broad-leaved forest; warm–temperate mixed coniferous and broad-leaved forest; cold–temperate, dark, coniferous forest; subcold shrub meadow; and cold meadow. Due to geological erosion and strong mountain seismic activity, many coseismic landslide events occurred, which further led to the frequent occurrence of debris flows in this basin and caused serious damage to local vegetation coverage and growth.

The earliest recorded debris flow occurred after the Diexi earthquake (1933), with a maximum discharge of 150 m³/s, and it washed away Xuehuaping village. From 1961 to 1978, 11 debris flows resulted in many casualties and severe economic losses. In 1978, the local government established a drainage channel network in the lower part of the gully that effectively reduced the frequency of debris flows in the basin. However, the “5.12” Wenchuan earthquake (M_w 7.5) in 2008 triggered landslides and collapses of many slopes, resulting in the evolution of the Qipan catchment into a high-frequency debris flow gully. From 8 to 12 July 2013, Wenchuan County suffered a heavy rainstorm (cumulative precipitation > 111.6 mm in 4 days) and triggered a catastrophic debris flow event [34]. After 2013, the government constructed a barrage across the Qipan catchment to prevent debris flows from running into the Minjiang River. However, three low-magnitude debris flows occurred on 5 July 2017, 22 August 2018, and 20 August 2019 (

Table 1. Historical debris flow events in the Qipan catchment, Wenchuan, China [35].

Data	Rainfall Intensity(mm)				Debris Flow Type	Peak Discharge (m3/s)	Duration (min)	Debris Flow Volume (104m3)
	3 Days	Daily	Hourly	10 min
1933	–	–	–	–	Viscous	150	–	–
6 July 1961	99.5	79.9	–	–		75	60	13.5
23 July 1964	48.3	41.7	–	1.2	Diluted	65	50	9.1
16 July 1965	69.5	41.2	–	–		65	50	9.9
28 July 1970	56.5	33.0	–	–		60	60	5.8
24 July 1971	79.4	53.4	–	–		62	45	8.4
29 July 1975	–	32.5	9.6	3.8		81	40	9.8
7 July 1977	–	39.4	7.6	1.6		65	30	5.8
15 July 1978	79.5	66.7	36.4	17.0	Viscous	90	50	13.5
15 August 1979	48.0	30.8	–	6.1	Diluted	42	30	3.8
26 July 1980	–	–	–	4.4		65	20	5.4
12 August 1981	–	53.8	9.5	2.1		90	25	6.7
19 July 1983	–	31.3	8.1	1.7		50	15	2.3
11 July 2013	109.6	54.3	6.4	–	Viscous	1745	30	78.2
5 July 2017	–	18.6	–	–		–	–	18.5
22 August 2018	–	33.4	–	–		–	–	11.5
20 August 2019	–	28.1	–	–		–	–	15

Note: “–” denotes missing data.

), and the basins upstream of the dams were nearly filled with debris flow material.

3. Data and Methods

3.1. Data Collection and Preprocessing

To account for factors such as resolution, cloud cover, and shadow, multitemporal, high-resolution remote sensing images of clear and cloudless periods in the growing season were obtained from Sentinel-2A, Gaofen-1, RapidEye, and Google Earth images (Landsat series), covering the time range from 2007 to 2019 (

Table 2. Data sources and descriptions.

Data Source, Accuracy	Type	Time	Acquisition Data
ALOS PALSAR, 12.5 m	DEM	2011	https://search.asf.alaska.edu/ accessed on July 2022.
Landsat 4, 30 m		2008/5–7	https://earthexplorer.nasa.gov/ https://earthexplorer.usgs.gov/ http://www.gscloud.cn/ accessed on June 2023.
Spot-5, 2.5–3.2 m		2008/6
Landsat 4, 30 m		2010/8
Rapideye-3A, 5 m		2012/7
Spot-5, 2.5–3.2 m	Remote sensing image	2012/7
Sentinel-2A, 10 m		2015/7
Sentinel-2A, 10 m		2018/8
Sentinel-2A, 10 m		2019/9
Gaofen-1, 1 m		2019/8	Beijing Digital Space Co., Ltd. (Beijing, China) accessed on June 2023.
LAADS DAAC, 250 m	MODIS Data	2007–2020	https://ladsweb.modaps.eosdis.nasa.gov/search/ accessed on September 2023.

), and the latest remote sensing images were obtained using unmanned aerial vehicles (UAVs) (Figure 2). After preprocessing the data, ENVI 5.3 and ArcMap 10.7 were used to process the image dataset.

Digital elevation models (DEMs) were constructed from Advanced Land Observation Satellite (ALOS) phased array L-band synthetic aperture radar (PALSAR) data with a resolution of 12.5 m for 2009, and a 1:10,000 isoline map produced by the Sichuan Surveying and Mapping Geographic Information Bureau in 2014 was used to extract the elevation, slope, aspect, flow direction, slope length, and other factors.

The land use data are based on Landsat 4-5TM, Landsat 8, and Sentinel-2A image data. Referring to the artificial classification criteria of land use types on the resource data cloud platform of the Chinese Academy of Sciences, combined with field investigations and UAV images [33], the land use type map of the Qipan catchment was obtained by using the ArcGIS 10.7 platform through computer and artificial visual interpretation.

The normalized difference vegetation index (NDVI), a measure of photosynthetic capacity, has been widely used to evaluate vegetation coverage [27]. Therefore, NDVI products are often used to assess vegetation destruction and its subsequent restoration process after earthquakes. The NDVI data used in this study were derived from the vegetation index product data (MOD13Q1) of the National Aeronautics and Space Administration (NASA) with spatial and temporal resolutions of 250 m and 16 days, respectively. The original data were preprocessed using the MODIS Reprojection Tool (MRT). The maximum value composite (MVC) method was used to further obtain the NDVI data at monthly, seasonal, and annual scales. The MVC method can reduce or eliminate the influence of clouds, atmosphere, and the solar altitude angle on images [36]. To eliminate the influence of seasonal changes in some vegetation from winter to summer, data from the growth period from May to September each year were selected as the representative data. In view of the high incidence of debris flow events during this period, NDVI data during the growth period can be used to accurately determine the vegetation status during the debris flow period in the basin. The main reason for using MODIS data instead of Landsat series data was that the latter have low temporal resolution, so there were not enough cloudless images covering the whole study area each year. To intuitively analyze the spatial and temporal changes in the study area, all spatial data (such as DEM, land use, and NDVI) were resampled using the nearest neighbor method with the same size and resolution as the study area.

3.2. Landslide Interpretation

Landslides directly or indirectly promote the occurrence of subsequent debris flows by changing the material supply conditions [37]. The transformation methods include dynamic transformation (direct transformation into debris flow during landslide movement) and static transformation (the evolution of landslide deposits is restarted again to form a debris flow) [38,39]. After investigation, was determined that the study area belongs to the static transformation mode. In this transformation process, there is a certain positive proportional relationship between the area of a landslide and the amount of recharged debris flow [33]. The scale of a debris flow can be estimated by using the area of the landslide and the connectivity of the landslide. A decrease in the area of the landslide body means that the supply to the debris flow source is reduced, resulting in a decrease in the probability of occurrence in the high-incidence area of the debris flow, thus effectively reducing the susceptibility of the debris flow and significantly enhancing the overall stability of the region [40,41]. Therefore, we conducted a field survey of the landslide area in the basin and extracted the landslide boundary from high-resolution images and drone images through computer and artificial visual interpretation. At the same time, the landslide range in the basin was determined by referring to the landslide data published by Fan et al. [42].

The 2008, the Wenchuan earthquake triggered large-scale collapses and landslides, resulting in significant damage to the ecosystem of the Qipan catchment. Over time, the ecosystem in the study area showed a gradual recovery trend [15,30]. In March 2013, a 4.5 magnitude earthquake further triggered many new landslides, providing a rich material source for debris flows. The heavy rainfall in the rainy season in July of the same year transformed a large amount of loose material on the slope into debris flow, which was the largest debris flow since 1933. Although there were several small-scale mudslides in the following years, the impact was relatively small [35]. By 2019, two medium-sized debris flows occurred again in the basin. A comparison of the remote sensing images revealed that vegetation restoration in the basin in 2019 was good, and most of the landslides disappeared [33]. Therefore, the selection of 2008, 2013, and 2019 as the key time nodes for the interpretation of landslide changes is based on the relatively uniform time intervals between these years, and relatively strong natural disasters occurred within three years. It is possible to accurately assess and understand the dynamic changes in landslides and the natural recovery process of a basin after major geological disasters, which is helpful for studying the impact of vegetation restoration on ecosystem stability.

According to the results of the field investigation and visual interpretation, the spatial distribution map of the induced landslides in the Qipan catchment after the earthquake was created (Figure 3). After the earthquake in 2008, 142 collapses and landslides were induced in the Qipan catchment (Figure 3a). In 2013, there were 134 landslides in the study area (Figure 3b). In 2019, the number of landslides decreased to 61 (Figure 3c). The interpretation results show that the restoration of vegetation significantly reduces the landslide area and the source of the landslide recharge debris flow, which, in turn, reduces the susceptibility to debris flow and gradually enhances the stability of the basin.

3.3. Sediment Connectivity

Borselli et al. [43] first proposed the index of connectivity (IC) to quantify sediment connectivity in 2008. The index comprehensively considers the characteristics of landscape topography and the impedance factors of water flow and sediment movement (such as terrain roughness and vegetation coverage). By calculating the potential migration probability of sediment from each position of the basin to the river channel, the spatial distribution of the potential sediment transport capacity in the basin is analyzed, and field-measured data are not needed [44]. After many improvements and developments, the connectivity index has become a flexible tool that can be adjusted according to different environmental conditions and can be used to effectively evaluate sediment transport and deposition processes [45,46]. The potential capacity for sediment transport affects the susceptibility of debris flows, thus affecting the stability of basin ecosystems. The numerical variation in this index is closely related to the environmental factors affecting runoff and deposition processes [47]. According to the probability of runoff from the upslope area and the probability of sediment transport to the downslope area [43], the corresponding formula is given, as follows:

$$IC={\mathit{log}}_{10}\left(\frac{D_{up}}{D_{down}}\right)={\mathit{log}}_{10}\left(\frac{\overline{W}\overline{S}\sqrt{A}}{{\displaystyle \sum \frac{d_i}{w_i{s}_i}}}\right)$$

(1)

where D_up is the potential for sediment detachment on the upslope, including soil erosion indices and topographic factors, and D_down represents the potential for sediment transport along the path to sinks, considering the flow path, flow distance, slope gradient, and land cover. $\overline{W}$ is the average weight factor (dimensionless); $\overline{S}$ and A represent the average slope angle (m/m) and upslope contribution area (m²) of the catchment, respectively; d_i is the length of the flow path along the ith cell (m); and w_i and s_i are the weight factors and the slope gradient of the ith cell, respectively. The IC ranges from −∞ to + ∞, with sediment connectivity increasing for larger IC values.

In the IC model, the weight factor $\overline{W}$ represents the sediment separation and transport resistance related to surface characteristics [47]. In previous studies, the vegetation cover and management factor (C value) in the universal soil loss equation (USLE) model directly replaced $\overline{W}$ in the IC model and provided a fixed experience value according to the land use type or tillage management system [48,49]. Considering that the land use in the study area has large spatial and temporal changes under the joint action of post-earthquake geological disasters and vegetation restoration, based on previous research results and land use changes, reference C factor assignments ranging from 0 (rock) to 1 (bare soil) were used to calculate the C factor values of different land use types [32,50,51], as shown in

Table 3. The C values of different land use types.

Land Use Type	C Value
Grassland	0.150
Alluvial fan	0.651
Farmland	0.250
Channel	0.001
Landslide	0.690
Built-up land	0.020
Forest	0.060
Barren land	0.651

.

3.4. Assessment of the Ecosystem Stability of Debris-Flow-Prone Areas

Basin ecosystem stability in debris-flow-prone areas is a prerequisite for ensuring the normal functioning of ecosystems and is closely related to sustainable livelihoods and development [2,52]. Debris flow, which is linked to topography, rainfall runoff paths, and land use or source materials, is the major factor driving regional differences in the stability of basin ecosystems. There are many methods for analyzing the potential danger of debris flows to ecosystem stability [53,54]. However, for debris-flow-prone areas, the degree of debris flow susceptibility and various factors affecting its susceptibility should be considered. These factors can be subdivided into dozens, but few can quantify the abovementioned factors uniformly.

Sediment connectivity is a key parameter reflecting the runoff–deposition process and indicates the sediment transport potential of a debris flow [43]. It represents the comprehensive effect of the runoff path, terrain elevation change, contribution area (which refers to the area of sediment supply in the basin and includes the effective supply area of landslides), slope change, land use type change, and other factors on the occurrence of debris flows [44,55,56]. As a key indicator of basin erosion and sediment movement characteristics, it has a direct impact on the occurrence probability and scale of debris flows. High sediment connectivity means a greater risk of debris flow, which threatens the stability of the basin [50,55,57]. Moreover, ecosystem stability cannot be separated from vegetation [58]. In an earthquake-disturbed basin, the revegetation zones in the landslide area provide much less material for debris flows than in the previous period when the vegetation was destroyed [59,60]. In short, in debris-flow-prone areas, it is meaningful to introduce vegetation evaluation into the variation in ecosystem stability. In this paper, we use the normalized difference vegetation index (NDVI), which is the most commonly used vegetation index and can better reveal the vegetation growth state and spatiotemporal changes in a basin, to evaluate vegetation [27,61]. Therefore, this study defines the ecosystem stability index based on the NDVI and IC, which is then used to comprehensively evaluate the ecological stability of debris-flow-prone areas after earthquakes. The evaluation process is provided in Figure 4.

The ecosystem stability index (ESI) of debris-flow-prone areas is clearly defined as a comprehensive calculation of the IC and NDVI. The index is used to represent the effect of vegetation interception and sediment transport potential on the supply of debris flows from watershed sources, to determine the dynamic changes in debris flow susceptibility after earthquakes and to characterize the stability of watershed ecosystems. The formula can be written as follows:

$$ESI={\left[\left(1-\frac{IC-{IC}_{\min }}{{IC}_{\max }-{\mathrm{IC}}_{\min }}\right)\times NDVI\right]}^{1/2}$$

(2)

Previous studies have shown that the NDVI has a negative correlation with the IC [31,62] and have noted that an NDVI value within 0.6–1 indicates good vegetation conditions, a value within 0.2–0.6 indicates medium vegetation conditions, and a value within 0–0.2 indicates poor vegetation conditions. A normalized IC within 0.8–1 indicates a strong sediment transport capacity in the basin, a value within 0.4–0.8 indicates a medium sediment transport capacity, and a value within 0–0.4 indicates a weak sediment transport capacity. Therefore, this paper comprehensively analyses the vegetation characteristics, debris flow characteristics, and sediment transport capacity in the basin and divides the ecosystem stability of the region into five grades. The specifics are shown in

Table 4. The threshold of the comprehensive ecosystem stability grade.

Level	Eco-Stability Index	Degree of Eco-Stability	Description
Ⅰ	0–0.2	Extremely unstable	Vegetation is sparse and sediment connectivity is high in the evaluation unit, NDVI: 0–0.2, and IC: 0.8–1.
Ⅱ	0.2–0.4	Unstable	NDVI: 0.2–0.4, and IC: 0.6–0.8.
Ⅲ	0.4–0.6	Substable	NDVI: 0.4–0.6, and IC: 0.4–0.6.
Ⅳ	0.6–0.8	Stable	NDVI: 0.6–0.8, and IC: 0.2–0.4.
Ⅴ	0.8–1	Extremely stable	NDVI: 0.8–1, and IC: 0–0.2.

.

4. Results

4.1. Spatial–Temporal Variations in Land Use

Figure 5 provides the land use map of the Qipan catchment in three different years. There are eight types, i.e., forest, grassland, farmland, gully, barren land, alluvial fan, landslide, and built-up land, in this basin.

From 2008 to 2019, forestland was dominant in the Qipan catchment, accounting for more than 52.14% of the total area of the basin. Alluvial fans are mainly distributed in the southeastern corner of the basin, accounting for approximately 13.18%; grasslands are scattered in ditches, accounting for approximately 10.73%; construction land and farmland are mainly concentrated in small areas at the outlets of the basin; and landslides are scattered on both sides of the channel. In terms of the spatial distribution pattern, the vegetation restoration activities in the first six years were mainly concentrated in the middle and lower reaches of the basin, after which the restoration activities gradually expanded to the upper reaches. At the same time, the vegetation coverage along the east–west boundary in the middle reaches of the basin remained stable and relatively high. As of 2019, the vegetation cover of the basin has been basically restored, and a significant reduction in the number of landslides has been observed.

Table 5. Land use transition matrix in the Qipan catchment from 2008 to 2013 (km²).

Year	Type	2013
		Grassland	Alluvial Fan	Farmland	Channel	Landslide	Built-Up Land	Forest	Barren Land	Total
2008	Grassland	1.82	1.04	0.28	0.12	0.20	0.00	2.07	0.08	5.61
	Alluvial fan	1.45	5.15	0.00	0.01	0.05	0.00	0.06	0.18	6.90
	Farmland	0.28	0.00	0.61	0.07	0.01	0.08	0.22	0.05	1.31
	Channel	0.00	0.00	0.00	0.25	0.04	0.01	0.18	0.00	0.48
	Landslide	0.20	0.00	0.00	0.40	1.30	0.01	3.24	0.10	5.24
	Built-up land	0.00	0.00	0.01	0.20	0.00	0.09	0.00	0.00	0.31
	Forest	1.12	0.07	0.04	0.99	2.07	0.00	23.18	0.16	27.63
	Barren land	1.18	0.29	0.02	0.25	0.54	0.01	2.45	0.13	4.87
	Total	6.05	6.57	0.96	2.29	4.21	0.20	31.40	0.70	52.36

,

Table 6. Land use transition matrix in the Qipan catchment from 2013 to 2019 (km²).

Year	Type	2019
		Grassland	Alluvial Fan	Farmland	Channel	Landslide	Built-Up Land	Forest	Barren Land	Total
2013	Grassland	2.03	1.31	0.17	0.04	0.03	0.01	1.93	0.53	6.05
	Alluvial fan	0.36	5.51	0.00	0.00	0.04	0.00	0.20	0.44	6.57
	Farmland	0.02	0.00	0.68	0.01	0.00	0.01	0.24	0.00	0.96
	Channel	0.42	0.01	0.01	0.58	0.08	0.08	1.03	0.09	2.29
	Landslide	0.63	0.06	0.00	0.02	0.89	0.00	1.98	0.62	4.21
	Built-up land	0.02	0.00	0.00	0.02	0.00	0.13	0.02	0.01	0.20
	Forest	3.34	0.10	0.04	0.13	0.23	0.01	27.04	0.53	31.40
	Barren land	0.14	0.18	0.02	0.00	0.01	0.01	0.23	0.11	0.70
	Total	6.97	7.17	0.92	0.80	1.28	0.25	32.66	2.33	52.36

and

Table 7. Land use transition matrix in the Qipan catchment from 2008 to 2019 (km²).

Year	Type	2019
		Grassland	Alluvial Fan	Farmland	Channel	Landslide	Built-Up Land	Forest	Barren Land	Total
2008	Grassland	1.96	0.77	0.24	0.03	0.07	0.00	2.16	0.38	5.62
	Alluvial fan	0.41	6.10	0.00	0.00	0.00	0.00	0.01	0.37	6.90
	Farmland	0.20	0.00	0.60	0.02	0.00	0.08	0.38	0.04	1.31
	Channel	0.04	0.00	0.00	0.21	0.02	0.01	0.20	0.00	0.48
	Landslide	0.73	0.04	0.00	0.20	0.43	0.01	3.52	0.31	5.24
	Built-up land	0.05	0.00	0.00	0.12	0.00	0.13	0.01	0.00	0.31
	Forest	2.44	0.01	0.05	0.17	0.56	0.00	23.62	0.78	27.63
	Barren land	1.14	0.25	0.03	0.05	0.20	0.01	2.75	0.44	4.87
	Total	6.97	7.17	0.92	0.80	1.28	0.25	32.66	2.33	52.36

show the land use transition matrix of the Qipan catchment over three years, which clearly shows the change trends among the different land use types. From 2008 to 2013, the land cover transformation was severe, and the areas of alluvial fans, landslides, bare land, and farmland decreased. Among them, the bare land area decreased by 4.17 km², up to 7.97%, and the decrease was the most obvious, while the forestland area increased most significantly, with an increase of 7.2% and an area of 3.77 km². The main reason for the increase in forestland and the decrease in bare land was the restoration of vegetation in the landslide zone and bare land after the earthquake. From 2013 to 2019, the scale of land use transformation slowed; the areas of farmland, gullies, and landslides decreased by 0.08%, 2.85%, and 5.59%, respectively; and the areas of forestland and grassland increased by 2.4% and 1.75%, respectively. In general, the land use types in the basin changed from bare land and landslides with low vegetation coverage to forestland and grassland with medium–high and high vegetation coverage.

After the 2008 Wenchuan earthquake (M_w 7.9), many geological disasters, such as collapses and landslides, were induced, and many collapse deposits were deposited in the main channel and on the slope surface. A total of 142 collapses and landslides were induced, accounting for 10.01% of the whole basin. From 2008 to 2019, after the Chinese government repaired the Qipan catchment sediment-retaining dam and restored the vegetation, the landslide area decreased significantly. Until 2019, the proportion of landslides decreased to 2.44%, while the total proportion of grassland and forestland increased from 63.50% in 2008 to 75.69%, and the vegetation coverage rate in the basin increased significantly.

4.2. Spatial–Temporal Variations in Vegetation Coverage

The NDVI satellite remote sensing inversion data from 2007 to 2020 were processed using the maximum value composite (MVC) method, and the annual average NDVI value of the Qipan catchment each year was obtained (Figure 6). Furthermore, the interannual variation in the NDVI in the five periods of spring, summer, autumn, winter, and the growing period was obtained. In this study, the growing period (May to September) on the western Sichuan Plateau was used as the research period.

The change in the normalized vegetation index (NDVI) reveals the change in vegetation area in the basin. The natural restoration of vegetation reduces the area and volume of landslides. Figure 6 shows that the vegetation coverage characterized by the NDVI in different periods of the basin was the highest in the growing period (0.72–0.83), followed by summer (0.71–0.83) and winter (0.53–0.65). The NDVI of each period decreased significantly in 2008, and the vegetation recovered rapidly 2 years after the earthquake. It began to increase slowly in 2010, and the NDVI decreased significantly during the two debris flow events of 2013 and 2019.

In the three growing seasons, the spatial distribution pattern in the basin is shown in Figure 7, which reveals the characteristics of higher values in the central area and lower values on both sides, where the two sides correspond to the gully mouth area and the accumulation fan area, respectively. There are several scattered low-value areas in the center of the basin, which are mainly caused by surface fragmentation caused by natural disaster events, resulting in land bareness and extensive destruction of vegetation. During the analysis period, the average NDVIs of the three years were 0.72, 0.77, and 0.79, respectively, indicating an interannual increase in vegetation coverage in the basin, which reflects the positive progress of regional vegetation restoration and growth. In addition, the change in the distribution of the NDVI frequency further emphasized the improvement in vegetation coverage each year.

Figure 8 shows that the area with NDVI values > 0.8 before the earthquake (representing high vegetation coverage) was dominant in the catchment, accounting for 76.06% of the total area, while the severe surface disturbance caused by the Wenchuan earthquake in 2008 reduced this area to 38.13%. The area with an NDVI between 0.4 and 0.8 increased by 34.67%, and the area with an NDVI < 0.4 increased by 3.24%. The surface vegetation in these areas has been seriously damaged, resulting in a significant decrease in the NDVI. Over time, the succession and restoration of damaged vegetation occurred, the soil stability and soil water holding capacity increased, and the growth rate of vegetation accelerated significantly. By 2013, the area with an NDVI > 0.8 was restored to 65.56%, and the area with an NDVI < 0.4 was reduced to 2.3%. In the later period, the vegetation recovery was relatively slow. By 2019, the area with an NDVI > 0.8 had recovered to 66.84%. Until 2020, the vegetation in the study area recovered to approximately the pre-earthquake level.

4.3. Catchment Connectivity

The IC values for the whole basin (Figure 9) in 2008, 2013, and 2019 were calculated using Formula 1, and the sediment connectivity of the basin showed a downward trend overall.

In 2008, the IC values ranged from −7.34 to 3.57, with an average value of −4.66, which was the highest of the three years (Figure 9a and Figure 10). Areas with high sediment connectivity were distributed in the basin outlet area and the river runoff path. After the occurrence of a debris flow in 2013, the IC values ranged from −8.30 to 3.57, with an average value of −5.67, which was lower than that in 2008. This was due to the landslide and collapse caused by the earthquake, which destroyed the surface, caused enormous amounts of gravel to be stored on the slope, accumulated rich material that blocked the channel, increased the impedance of sediment movement, and reduced sediment connectivity. In 2019, the IC value continued to decline, ranging from −9.61 to 3.57, with an average value of −6.93. This change may have occurred because after the debris flow in 2013, the government-built, sand-retaining dam greatly reduced the transport of sediments, and the rapid restoration of slope vegetation increased the roughness of the slope and weakened the sediment transport capacity of the basin.

The IC values of 15 tributaries in the basin were unevenly distributed, as shown in Figure 11 and Figure 12 and

Table 8. IC values in 2008, 2013, and 2019.

Region	IC Values
	2008			2013			2019
	Max	Min	Mean	Max	Min	Mean	Max	Min	Mean
Total region	−0.106	−7.345	−4.75	−2.492	−8.283	−5.748	−2.492	−9.612	−7.024
Tributary #1	−0.106	−5.849	−4.21	−2.492	−6.935	−5.280	−2.492	−8.145	−6.449
Tributary #2	−2.052	−6.011	−4.48	−3.215	−7.180	−5.652	−4.418	−8.363	−6.840
Tributary #3	−2.298	−6.721	−4.70	−3.358	−7.812	−5.791	−4.647	−9.048	−7.011
Tributary #4	−2.513	−6.334	−4.83	−3.552	−7.382	−5.871	−4.861	−8.595	−7.083
Tributary #5	−2.493	−6.366	−4.79	−3.480	−7.346	−5.756	−4.802	−8.654	−7.065
Tributary #6	−2.600	−5.941	−4.82	−3.548	−6.912	−5.761	−4.876	−8.246	−7.083
Tributary #7	−2.980	−6.115	−4.81	−3.922	−7.117	−5.758	−5.267	−8.411	−7.096
Tributary #8	−3.091	−7.229	−4.86	−4.024	−8.174	−5.797	−5.358	−9.508	−7.129
Tributary #9	−3.020	−7.345	−4.89	−3.948	−8.283	−5.814	−5.281	−9.612	−7.142
Tributary #10	−3.075	−6.751	−4.74	−4.013	−7.695	−5.674	−5.352	−9.030	−7.009
Tributary #11	−2.683	−6.011	−4.82	−3.621	−6.926	−5.758	−4.958	−8.259	−7.094
Tributary #12	−3.012	−6.181	−4.84	−3.950	−7.119	−5.768	−5.287	−8.450	−7.097
Tributary #13	−2.740	−6.218	−4.74	−3.783	−7.270	−5.793	−5.084	−8.504	−7.032
Tributary #14	−3.117	−5.858	−4.57	−4.198	−6.935	−5.645	−5.465	−8.187	−6.904
Tributary #15	−2.888	−5.512	−4.58	−3.995	−6.614	−5.684	−5.242	−7.855	−6.925

. In the three years, the IC values of tributaries near the outlet (downstream) of the basin were greater than those of other tributaries, and the IC values of tributary 1 closest to the downstream area were the highest, with annual average values of −4.21, −5.28 and −6.44. The IC value of tributary 9, which is the farthest from the outlet of the basin, was the smallest, with average annual values of −4.89, −5.81 and −7.14. These data indicate that in the upstream region the surface roughness is greater due to steep terrain, exposed rock, and strong erosion. This slows the flow of water and limits the migration and transport of sediment, thereby reducing the connectivity of the region. Relatively speaking, the terrain of the downstream area is flat, the sediment coverage is greater, and the surface roughness is lower, which makes it easier for the water flow to carry and transport sediment, thereby increasing the connectivity. Therefore, the sediment connectivity upstream is lower than that downstream.

4.4. Relationship between Vegetation and Connectivity in the Landslide Area

Pearson’s correlation was used to analyze the relationship between the IC and NDVI. As shown in Figure 13, there was a significant negative correlation between the NDVI and sediment connectivity (R < −0.3, p < 0.01), which indicates that in the study area, with increasing vegetation status, sediment connectivity tended to decrease. A negative correlation coefficient R indicates that there is a stability mechanism between an increase in vegetation and a decrease in sediment. The regeneration of vegetation increases the surface roughness, hinders the solid material transport process of the landslide, and reduces the potential sediment transport capacity of the channel. Low vegetation coverage, steep slopes, and high annual average rainfall provide favorable environmental factors and power sources for the transport of loose materials, thereby enhancing the potential sediment transport capacity of the channel. Therefore, the presence of vegetation can effectively reduce the connectivity of sediments and slow the development of debris flows. This result is highly important for understanding ecosystem health and stability in debris-flow-prone areas.

4.5. Comprehensive Evaluation of Ecosystem Stability in the Qipan Catchment

Figure 14 shows the results of the ecosystem stability assessment for the Qipan catchment in 2008, 2013, and 2019. The overall stability of the basin has gradually increased over the past 11 years, especially in 2019. The average ecosystem stability index values for the whole basin from 2008 to 2019 were 0.617, 0.642, and 0.685, respectively, showing an overall increasing trend. The study area was mainly stable and substable, while a few unstable areas were mainly concentrated in the outlet (downstream) and southeast corner (upstream) of the basin, and only 0.10% were extremely unstable areas.

From 2008 to 2013, the ecosystem stability state mostly recovered from a substable to stable state, and the ecologically substable area decreased from 33.99% to 27.70% (Figure 15). The ecosystem stability area increased from 63.11% to 70.47%, and the areas with increasing stability were mainly concentrated in the downstream gully. Compared with the results from the land use map (Figure 5), this area was mainly in the area of rapid vegetation recovery on the landslide, indicating that the ecological stability was closely related to vegetation restoration activities. From 2013 to 2019, the ecological stability increased more rapidly. The ecologically extremely stable area increased from 0.52% to 3.73%, the ecologically stable area increased by 6.74%, and the ecologically substable area decreased by 9.49%, indicating that the government’s construction of sand dams greatly reduced the transport of sediment, and the rapid restoration of slope vegetation had a greater effect on debris flow susceptibility; thus, the stability of the basin greatly improved.

Figure 16 shows the average ecosystem stability of each tributary. The areas with low stability are mostly located in the outlet (downstream) and southeast corner (upstream) of the basin. The reason for the low stability of the downstream area is that there are more human activities in the region, as well as low vegetation coverage and high sediment connectivity, while the upstream area is mostly comprised on debris flow accumulation with less vegetation coverage, resulting in the lowest stability of the region. The stability of D4 in the basin was the highest of the three years because of its high vegetation coverage and the small landslides in this area. From 2008 to 2013, the stability of only six branch ditches improved, but the stability of most branch ditches changed only slightly. From 2013 to 2019, except for the stability of D13, the stability of the other branch ditches increased, among which D5 and D9 exhibited the greatest increase. The stability recovery state of each tributary in the basin is different. Therefore, when implementing watershed stability management measures, it is also necessary to consider and coordinate the differences in ecological stability restoration among tributaries.

5. Discussion

In this study, we found that in the Qipan catchment the NDVI immediately decreased sharply after the earthquake and gradually increased in the following years (Figure 6, Figure 7 and Figure 8). This finding is consistent with the research results of Zhang et al. [63], Shen et al. [64], Wu et al. [61], and others, indicating that coseismic landslides caused large areas of local vegetation to peel off, causing severe interference in the basin ecology, and these areas were prone to debris flows. Studies have confirmed that the decrease in landslide activity after an earthquake is related to the increase in vegetation restoration [65,66]. Moreover, we found that changes in sediment connectivity can reflect the influence of geological and human activities on sediment movement [47]. This change reveals the health status of the basin and the evolution of the landform after the debris flow disaster, which has a certain correlation with the stability of the basin. Some scholars’ research also proves that the stability of a disaster area is related to the connectivity [52,53,55]. Therefore, when evaluating the stability of basin ecosystems, it is helpful to comprehensively understand basin stability by considering various factors such as vegetation and sediment connectivity. Therefore, this study combines vegetation and sediment connectivity to construct a stability assessment model.

To ensure the reliability and accuracy of the evaluation method, a two-step verification strategy was adopted in this study. First, the consistency and accuracy between the remote sensing images and ecosystem stability assessment results for each year were verified via spatial comparison. A case analysis of the evaluation results of the Qipan catchment in 2013 was carried out (Figure 17). All the landslide-sensitive areas are clearly marked in the figure. The ESI range of all landslide areas in the basin are 0.403–0.734, and the average ESI is 0.643. The stability is still low. Three landslide points were selected for schematic comparison. The average ESI values of these three schematic areas are 0.602, 0.599 and 0.616 from left to right, which are in the substable stage and are highly consistent with the actual land use situation reflected in remote sensing images, thus verifying the rationality of the evaluation results.

Further verification involves an in-depth analysis of specific landslide areas, selecting the field-tested Laoyingyan landslide area (Figure 18a–f) and the largest landslide area (Figure 18g–l) in the basin for representative analysis. By comparing the evaluation results with the actual terrain dynamics, it was found that there is a high consistency. Although there is incomplete correspondence between the actual situation of the terrain and the evaluation results due to the influence of the abrupt change in the terrain, this consistency reveals the change in stability of the landslide area in the debris-flow-prone area and confirms the effectiveness of the evaluation model in capturing the dynamic changes in the ecosystem in the debris-flow-prone area. In addition, the results emphasize that an increase in vegetation can effectively slow the development speed of debris flows, thereby reducing their susceptibility [58,59]. Therefore, the research method in this paper is highly suitable for debris-flow-prone areas and can provide substantive guidance for ecological restoration and vegetation management in debris-flow-prone areas.

6. Conclusions

Based on the unique sediment connectivity of debris-flow-prone areas, a new method for evaluating the ecosystem stability of debris-flow-prone areas after earthquakes is proposed in this paper. This study comprehensively considered sediment connectivity and regional vegetation conditions, and the stability of the ecosystem in the basin was quantitatively evaluated by taking the Qipan catchment, a post-earthquake debris-flow-prone area in Wenchuan County, as an example. The following conclusions have been obtained from our study:

• The number of collapse and landslide events in the Qipan catchment in 2008, 2013, and 2019 was 142, 134, and 61, respectively, with a significant reduction in the landslide area. The significant reduction in this landslide event directly led to a decrease in the supply of debris flow sources, which effectively reduced the susceptibility to debris flows and enhanced the overall stability of the basin.
• From 2008 to 2019, the spatial pattern of the land use types in the Qipan catchment changed significantly: the vegetation coverage in the study area continued to increase, while the sediment connectivity showed a downward trend. In areas with low vegetation coverage, the IC value is relatively high, and there is a significant negative correlation between them. The restoration of vegetation effectively reduces the material source of debris flows through the interception of sediment, thus slowing the formation and development of debris flows.
• From 2008 to 2019, the ecosystem stability of the whole study area significantly improved, and its distribution was characterized by lower stability in the northern and southern regions and higher stability in the central region. The evaluation results of ecosystem stability are consistent with the surface evolution trend of debris-flow-prone areas, which confirms the applicability of this evaluation method in the evaluation of ecological stability in post-earthquake debris-flow-prone areas.

In this study, by considering the topographic characteristics and vegetation restoration speed of post-earthquake debris-flow-prone areas combined with traditional stability assessment methods, the NDVI and IC were used to assess the ecological stability of post-earthquake debris-flow-prone areas. Through field tests and comparisons, it was found that the evaluation method has high accuracy and strong applicability and can accurately reflect the stability change characteristics of post-earthquake debris-flow-prone areas. This assessment method not only considers the impact of the external representation of the ecosystem (vegetation coverage) on stability but also comprehensively considers the internal factors (sediment connectivity), thus providing a more comprehensive assessment perspective. In future research, other factors that may directly or indirectly affect the stability of the ecosystem could also be considered as evaluation factors to be added to the assessment of ecosystem stability in post-earthquake debris-flow-prone areas, thereby improving the depth and breadth of the assessment. The assessment of ecosystem stability in debris-flow-prone areas can provide more accurate and comprehensive information support for local governments and environmental protection agencies so that they can take targeted actions to help debris-flow-prone areas quickly improve regional stability.