Debris Flow Risk Assessment Based on a Water–Soil Process Model at the Watershed Scale Under Climate Change: A Case Study in a Debris-Flow-Prone Area of Southwest China

Risk assessment lays a foundation for disaster risk reduction management, especially in relation to climate change. Intensified extreme weather and climate events driven by climate change may increase related disaster susceptibility. This may interact with exposed and vulnerable socioeconomic systems to aggravate the impacts and impede progress towards regional development. In this study, debris flow risk under climate change was assessed by an integrated debris flow mechanism model and an inclusive socioeconomic status evaluation. We implemented the method for a debris flow-prone area in the eastern part of the Qinghai-Tibet Plateau, China. Based on the analysis of three general circulation models (GCMs)—Beijing Climate Center Climate System Model version 1 (BCC_CSM), model for Interdisciplinary Research on ClimateEarth System, version 5 (MIROC5, and the Community Climate System Model version 4 (CCSM4)—the water–soil process model was applied to assess debris flow susceptibility. For the vulnerability evaluation, an index system established from the categories of bearing elements was analyzed by principle component analysis (PCA) methods. Our results showed that 432 to 1106 watersheds (accounting for 23% to 52% of the study area) were identified as debris-flow watersheds, although extreme rainfall would occur in most of the area from 2007 to 2060. The distributions of debris flow watersheds were concentrated in the north and transition zones of the study area. Additionally, the result of the index and PCA suggested that most areas had relatively low socioeconomic scores and such areas were considered as high-vulnerability human systems (accounts for 91%). Further analysis found that population density, road density, and gross domestic production made great contributions to vulnerability reduction. For practical mitigation strategies, we suggested that the enhancement of road density may be the most efficient risk reduction strategy.


Introduction
Among the greatest continuing threats to mountainous areas, rainfall-induced debris flows are particularly hazardous for the direct losses of life and property and ready mobilization into chain-styled catastrophes [1][2][3]. The risk of rainfall-induced debris flows will be increased by the underlined drivers of climate change and variability [4][5][6][7][8][9]. According to the IPCC Fifth Assessment Report (AR5) [10], the impacts of climate change may be profound in the future with an associated increase flow dynamics model, and on this basis, the related future risk was assessed through a comprehensive vulnerability evaluation of a debris-flow-prone area of southwest China, with an area of approximately 112,710 km 2 . Different from previous research, each watershed in the study area was simulated by a debris flow trigger mechanism model to assess whether it is a potential debris flow watershed. Though more precise debris flow trigger mechanism models at the gully scale exist, this study employed a water-soil coupling mechanism model for the following reasons: The ability to synthesize the genetic conditions of the debris flow at a watershed scale [46,47]; interaction with the vulnerability of elements at the township scale; and the development of reasonable risk reduction measures [37]. Meanwhile, we aimed to analyze the impacts of debris flows on societal and economic developments in mountain areas. The major models, procedures, and evaluation included: (1) Downscaling climate data and assessing extreme rainfall events in the future; (2) predicting the distribution of debris flow susceptibility by a water-soil coupling mechanism model; and (3) analyzing the vulnerability of disaster-prone watersheds and assessing the ultimate risk.

Study Area
The debris flow-prone area is a part of the Hengduan Mountains, located in the west of Sichuan Province, China ( Figure 1). Covering 112,710 km 2 , this area is mainly comprised of hills, mountains, and a plateau while more than 95% is occupied by mountains. It is in one of the most tectonically active regions in the world. Among the historical recorded debris flow gullies, this area occupies about 48% (1575) of Sichuan province (3200). It has a subtropical monsoon climate with humid summers and dry winters. The spatial patterns of rainfall are remarkably different within the debris flow-prone area. A combination of the complex topography and monsoon climate creates a higher frequency of extreme rainfall. Meanwhile, this area is sensitive to climate change and the impacts of climate change may be profound due to the complex topography. Moreover, numerous loose deposits have been produced recently because of the 5.12 Wenchuan earthquake and the 4.20 Lushan earthquake, which have drastically increased the risk of rainfall-triggered debris flows in the study area [48,49]. For these reasons, a risk assessment of debris flow under climate change is imperative in this debris flow-prone area. area of approximately 112,710 km 2 . Different from previous research, each watershed in the study area was simulated by a debris flow trigger mechanism model to assess whether it is a potential debris flow watershed. Though more precise debris flow trigger mechanism models at the gully scale exist, this study employed a water-soil coupling mechanism model for the following reasons: The ability to synthesize the genetic conditions of the debris flow at a watershed scale [46,47]; interaction with the vulnerability of elements at the township scale; and the development of reasonable risk reduction measures [37]. Meanwhile, we aimed to analyze the impacts of debris flows on societal and economic developments in mountain areas. The major models, procedures, and evaluation included: (1) Downscaling climate data and assessing extreme rainfall events in the future; (2) predicting the distribution of debris flow susceptibility by a water-soil coupling mechanism model; and (3) analyzing the vulnerability of disaster-prone watersheds and assessing the ultimate risk.

Study Area
The debris flow-prone area is a part of the Hengduan Mountains, located in the west of Sichuan Province, China ( Figure 1). Covering 112,710 km 2 , this area is mainly comprised of hills, mountains, and a plateau while more than 95% is occupied by mountains. It is in one of the most tectonically active regions in the world. Among the historical recorded debris flow gullies, this area occupies about 48% (1575) of Sichuan province (3200). It has a subtropical monsoon climate with humid summers and dry winters. The spatial patterns of rainfall are remarkably different within the debris flow-prone area. A combination of the complex topography and monsoon climate creates a higher frequency of extreme rainfall. Meanwhile, this area is sensitive to climate change and the impacts of climate change may be profound due to the complex topography. Moreover, numerous loose deposits have been produced recently because of the 5.12 Wenchuan earthquake and the 4.20 Lushan earthquake, which have drastically increased the risk of rainfall-triggered debris flows in the study area [48,49]. For these reasons, a risk assessment of debris flow under climate change is imperative in this debris flow-prone area.

Assessment Framework
The assessment of debris flow risk was comprised of four steps, which are illustrated in Figure 2. The assessment of debris flow risk was comprised of four steps, which are illustrated in Figure  2.

Debris Flow Susceptibility
A water-soil coupling model was used for the debris flow susceptibility assessment [50]. Of all types of debris flow, shallow soil slopes sliding along bedrock is the most common in mountain areas. Furthermore, such a type is described as thick over soil failure caused by shear breakage under major exogenic forces, like precipitation, and is substantially smaller in scale, yet larger in number and greater in devastation. The formations of this type of debris flow depend on the energy, material, and triggering conditions. The relative height and terrain slope provide energy conversion gradients, while extreme rainfalls work as a triggering condition. However, it is difficult to estimate accurately the loose material of each plot for a large area. This problem can be solved by evaluating the loose material generation in each plot. In this way, the material condition of debris flow formation is influenced by geological structure, land cover, and human activity. Thus, the water-soil coupling model was established for such types of debris flow based on a geomorphology-based distributed hydrological model and the slope instability mechanism model. The water-soil coupling process of debris flow formation can be described as having three basic steps. The first step is to simulate the surface runoff formation process. Due to the large surface slope of the mountain area and the weak river recharge of groundwater, the 1-D Richards equation set up by the finite element method was applied to present water flow in unsaturated soils for each raster, as shown in Equations (1)-(3). Surface runoff yield was based on the excess infiltration mechanism and calculated by Manning's formula [51,52]. Secondly, based on the limit equilibrium method of infinite slopes, the slope stability alongside the instability depth under rainfall infiltration was analyzed by the safety factors calculation model as expressed in Equation (5) for each raster. If Fs < 1, the soil mass is unstable and about to transform into sediments. Third, the density of the water-soil mixture was quantified by the volume of surface runoff and the gross amount of unstable soil mass at the watershed scale. With this, a qualitative relationship between the density and probability of debris flow was presented. The density of the water-soil mixture was calculated using Equation (6):

Debris Flow Susceptibility
A water-soil coupling model was used for the debris flow susceptibility assessment [50]. Of all types of debris flow, shallow soil slopes sliding along bedrock is the most common in mountain areas. Furthermore, such a type is described as thick over soil failure caused by shear breakage under major exogenic forces, like precipitation, and is substantially smaller in scale, yet larger in number and greater in devastation. The formations of this type of debris flow depend on the energy, material, and triggering conditions. The relative height and terrain slope provide energy conversion gradients, while extreme rainfalls work as a triggering condition. However, it is difficult to estimate accurately the loose material of each plot for a large area. This problem can be solved by evaluating the loose material generation in each plot. In this way, the material condition of debris flow formation is influenced by geological structure, land cover, and human activity. Thus, the water-soil coupling model was established for such types of debris flow based on a geomorphology-based distributed hydrological model and the slope instability mechanism model. The water-soil coupling process of debris flow formation can be described as having three basic steps. The first step is to simulate the surface runoff formation process. Due to the large surface slope of the mountain area and the weak river recharge of groundwater, the 1-D Richards equation set up by the finite element method was applied to present water flow in unsaturated soils for each raster, as shown in Equations (1)-(3). Surface runoff yield was based on the excess infiltration mechanism and calculated by Manning's formula [51,52]. Secondly, based on the limit equilibrium method of infinite slopes, the slope stability alongside the instability depth under rainfall infiltration was analyzed by the safety factors calculation model as expressed in Equation (5) for each raster. If F s < 1, the soil mass is unstable and about to transform into sediments. Third, the density of the water-soil mixture was quantified by the volume of surface runoff and the gross amount of unstable soil mass at the watershed scale. With this, a qualitative relationship between the density and probability of debris flow was presented. The density of the water-soil mixture was calculated using Equation (6): In Equation (1), k(θ,z) is the unsaturated soil hydraulic conductivity, and its concrete expression is shown in Equation (2). k 0 is the saturated water conductivity of the surface soil layer. S e is the saturated degree as described by Equation (3). q v (θ,z) is the soil water flux. ψ = (u a − u w ) is the matric suction of the soil from the Van Genuchten model of Equation (4) (u a represents the atmospheric pressure and was set to zero) [31]. n and m are the parameters of the curve, in accordance with n = 1 − 1/m. θ s is the saturated water content of the soil, θ r is the residual water content of the soil, and θ ∆t is the soil water content of the current hour.
To calculate the safety factors in Equation (5), ϕ is the internal friction angle, β is the gradient, c is the soil cohesion force, ϕ b is related to the matrix suction, γ t is the wet soil density at time t, H s is the soil depth, and S e is the same as above. In Equation (6), w is the density of the water ( w = 1.0 g/cm 3 ), and s is the density of solid particles ( s = 2.7 g/cm 3 ). v w is the volume of the water and v s is the volume of the failure soil mass.
According to these equations, the probability of debris flow can be predicted for each day in watershed units. In order to reduce uncertainties, the threshold value of the probability of debris flow was set to 0-20%, 20-40%, 40-60%, 60-80%, and 80-100%. Then, the debris flow susceptibility was expressed by the number of events for the future.

Vulnerability
The vulnerability of a local debris flow disaster system varies between research targets and across scales, depending on the environmental governance, demographic features, and development status [53]. In this study, vulnerability was expressed as a comprehensive index estimated through principal component analysis (PCA) and set as a range from 1 to 100. the construction of an evaluation index system is a commonly used method to integrate different categories of affected elements, to synthesize the vulnerability index by statistical methods or other mathematical methods, and to represent the relative degree of vulnerability of the assessment unit [54][55][56][57]. Considering the availability and consistency of socio-economic data in debris flow watersheds which normally correspond to the township scale, 11 variables (Table 1) were included and were mostly inferred from previous representative studies [58][59][60]. The procedures of PCA are described by Jolliffe [61]. In order to facilitate local disaster mitigation and guide decision-making, vulnerability assessment was established in the current state. Table 1. Summary of the variables used in the comprehensive vulnerability index.

Variables
Description Methods and Source Population density D pop = P/S

Non-agricultural work force
The number of people engage in work other than agriculture, i.e., industrial and service employees Census data (2014)

Town-constructed area
Area of regions with complete infrastructure facility and public utilities in town Census data (2014)

population proportion
The proportion of population of a certain town to that of the superior county B = P/P 0

Population concentration
The relative importance of population distribution for a certain town R = B/(S/S 0 ) Traffic net density D ti = L ti /S

Traffic trunk influence
Influence of different types of traffic trunk lines

Primacy index
Primacy index means the transportation radiation force of a primate city and represents the social and economic conditions of towns been affected by the primate cities.
,1} parameters D pop -density of population, P-the permanent population of a town, P 0 -the permanent population of a county D ip -density of industrial production, V-annual industrial production, S-Town area, S 0 -county area D ti -density of traffic network, L ti -line length or traffic nodes of traffic network C im -influence of m traffic trunk lines, I-traffic type (Roads, railways, airports, etc.), m-traffic trunk level u-minimum line length from e to f. If l ef ≤ 100 km, traffic type is set to roads, otherwise railways. Then, the shortest path function, f e (x), is determined by comparing the minimum value of the path from the primate city to region e. According to the function, an important node, f, is selected as a junction node for region e. H i is a Boolean value. When the f e (x) function is true, H i is 1 or vice versa. Using the function in a GIS software, the hinterland range of the important node, f, H m (i) is defined through search and comparison.
In order to tailor targeted disaster adaptation, the variables retaining a large proportion of the principal components were analyzed separately to avoid any possible unaffordable adverse impacts of debris flow. Lower vulnerability scores mean a more resilient watershed to debris flow risk.

Debris Flow Risk Assessment
According to the risk assessment framework [62], the debris flow risk is equal to the product of debris flow susceptibility and vulnerability as shown in Equation (7):

Data
The debris flow risk assessment is closely related to both natural conditions and socio-economic development level, which mainly involve datasets about precipitation, lithology, elevation, soil, land use, road, population, and GDP.

Daily Precipitation Data
As input data of the water-soil coupling model, daily rainfall projections were obtained from general circulation models (GCMs) of the NASA Earth Exchange Global Daily Downscaled Projections (NEX-GDDP) dataset (https://cds.nccs.nasa.gov/nex-gddp/) in RCP4.5 emissions scenarios during a time period of 2007-2060.
Since uncertainties always exist in climate projection, different GCMs can provide a more useful reference for debris flow susceptibility prediction so that the uncertainties in climate projection can be revealed and reduced. Three GCMs from CMIP5 were chosen to project future changes in precipitation, including the Beijing Climate Center Climate System Model (BCC_CCSM), the Model for Interdisciplinary Research on Climate, version 5 (MIROC5), and the Community Climate System Model version 4 (CCSM4).
Additionally, the results from NEX-GDDP could not be directly used in the water-soil coupling model due to the coarse spatial resolution (25 km × 25 km). Downscaled GCMs' projection was capable of representing the long-term ground observations patterns. The previous literatures also indicates that downscaled GCMs have great similarity to local extreme events in the flood season [63,64]. Thus, downscaling techniques were needed to deal with raw data. In general, the local climate is closely related to land surface, elevation, and other factors. In most previous studies, linking of the statistical downscaling and spatial interpolation method was an effective approach to obtain high-resolution rainfall data [65]. In this study, the statistical downscaling approach proposed by Holland [66] was used, which used a linear correction algorithm to correct the bias of the mean climate in GCMs data. Moreover, the ANUSPLIN software was applied to interpolate the rainfall surface at a 1 km scale in terms of the bias-corrected data. Considering that an extreme rainfall event is a major driver of debris flow susceptibility, the indices of intensity, duration, and frequency were selected to analyze the differences of extreme rainfall and debris flow susceptibility.
In detail, historical data collected by the local meteorological stations in and around the debris flow-prone area (Figure 1) from 1981 to 2010 were provided by the National Meteorological Information Center of China (http://www.cma.gov.cn).

Environment Data
The water-soil coupling model employed multiple data, including land use, elevation (relative height and slope), surface soil, and a lithology map. Land use data were interpreted from remote sensing data of Landsat 7, 2015 (https://earthdata.nasa.gov/). A digital elevation model (DEM) of a 30-m resolution was obtained from the United States Geological Survey (https://earthexplorer.usgs.gov/). Soil datasets including soil type were provided by the institute of soil science, Chinese Academy of Sciences, and were converted into a 1 km resolution. The lithology distributions were obtained from a regional geological map while the soil cohesion, c, internal friction angle, ϕ, and other mechanical parameters were in accordance with a rock mechanics handbook. The relationship between the matrix suction and the soil water content was calculated by the Van Genuchten model [67]. The parameters were also determined.
All data were converted to a 100 m resolution raster in order to calculate in the assessment model.

Projection of Extreme Rainfall
As shown in Figure 3, the spatial variability of extreme rainfall was rather large, as indicated by the different extreme rainfall indices. This was probably a reflection of the impacts of the complex topography and monsoons on extreme rainfall events. High frequencies of extreme rainfall were found in the transition zone from the basin to the plateau. The extreme precipitation index from BCC_CSM had larger differences, while the CCSM4 model had a higher value than the other two. Spatially, extreme rainfall was mainly concentrated in the southwest of this area. In detail, almost the whole area would experience more than 10 extreme rainfall events of Rd25 (Figure 3a-c), and about 92% of the area would contain more than 100 extreme rainfall events of R5D50 during the time period of 2007 to 2060 (Figure 3d-f). The maximum one-day precipitation ranged from less than 40 mm to more than 400 mm, while the maximum intensity was about 406 mm. About 96% of the area would experience extreme rainfall events with a total accumulated precipitation of more than 40 mm. Results from the extreme rainfall showed that areas with a maximum daily precipitation more than 50 mm are concentrated in the western part of the study area (Figure 3g-f). This suggests that the increase of extreme rainfall events would increase the possibility of debris flow events in the study areas under climate change, while the number of historical debris flow events was more than 10 in most watersheds.

The Distribution of Debris Flow Susceptibility
Applying the water-soil coupling mechanism model to each watershed obtained results for different levels of debris flow susceptibility. Results from the BBC_CSM model showed that watersheds with a probability of debris flow above 80% had the largest range and contained those with a lower probability. Thus, watersheds with an 80% probability were used for further analysis in order to reduce uncertainties. In total, 432 of all watersheds were identified as debris flow watersheds and occupied 23% of the whole study area. High-risk areas (the number of events with debris flows >30 t) were mainly distributed in the counties of Mianzhu, Wenchuan, Leshan, and Mianning (accounting for 29% of all debris flow watersheds). As extreme rainfall events were inherently random, moderate and low-risk watersheds had greater quantities (accounting for 87% of all debris flow watersheds). Spatially, watersheds with a moderate and low-risk were more broadly distributed than high-risk watersheds, while the former were concentrated in the northern part of the study area. Additionally, the spatial distributions of potential debris flow simulated from MIROC5 and CCSM4 data were basically the same. In total, 1106 watersheds (accounting for 52%) were forecasted as debris flow watersheds. Although the hazard areas were larger, either probabilities or frequencies of future debris flow simulated by the MIROC5 and CCSM4 models were lower than the BCC_CSM model. Comparison of the distributions of extreme precipitation with the debris flow susceptibility generally showed agreement in the debris flow-prone area. However, the mechanism model helped to identify

The Distribution of Debris Flow Susceptibility
Applying the water-soil coupling mechanism model to each watershed obtained results for different levels of debris flow susceptibility. Results from the BBC_CSM model showed that watersheds with a probability of debris flow above 80% had the largest range and contained those with a lower probability. Thus, watersheds with an 80% probability were used for further analysis in order to reduce uncertainties. In total, 432 of all watersheds were identified as debris flow watersheds and occupied 23% of the whole study area. High-risk areas (the number of events with debris flows >30 t) were mainly distributed in the counties of Mianzhu, Wenchuan, Leshan, and Mianning (accounting for 29% of all debris flow watersheds). As extreme rainfall events were inherently random, moderate and low-risk watersheds had greater quantities (accounting for 87% of all debris flow watersheds). Spatially, watersheds with a moderate and low-risk were more broadly distributed than high-risk watersheds, while the former were concentrated in the northern part of the study area. Additionally, the spatial distributions of potential debris flow simulated from MIROC5 and CCSM4 data were basically the same. In total, 1106 watersheds (accounting for 52%) were forecasted as debris flow watersheds. Although the hazard areas were larger, either probabilities or frequencies of future debris flow simulated by the MIROC5 and CCSM4 models were lower than the BCC_CSM model. Comparison of the distributions of extreme precipitation with the debris flow susceptibility generally showed agreement in the debris flow-prone area. However, the mechanism model helped to identify a low debris flow susceptibility under high frequencies of extreme rainfall events, such as the northwestern part of this area, by comparison between Figures 3 and 4. Further analysis suggested that the rainfall threshold of debris flow varied from 30 mm to 50 mm in Figure 4d, indicating that the threshold value varies from different watersheds.

Vulnerability and Debris Flow Risk Assessment.
Three principal components were retained in PCA for the whole debris-prone area regional analysis. Together these first three principal components accounted for 73.03% of the variation in the original 11 variables included in the analysis. The factor loading was dominated by the population and workforce. The loadings of each variable for the retained principal components are detailed in Table 2. As shown in Table 2, the first component was heavily loaded on population and workforce. The second PCA was mainly loaded on transportation related variables. The third component was loaded on production-related variables. Population density, road density, and gross domestic production with the heaviest loads of each component were thus selected for further analysis.
A measurement of vulnerability was calculated by integrating different data and was standardized to a range of 1 to 100 (Figure 5a), where higher values indicated a higher vulnerability towards debris flow. Consistent with relevant studies, the high vulnerability area towards debris flow (vulnerability value >90) comprised the most area because of the relatively low social-economic conditions in mountain areas, which accounted for over 91% of the study area. However, the level of vulnerability could also determine the debris flow risk and management strategies. Moreover, the number of debris flow events had an approximate number by comparing Figure 5c with Figure 5b.

Vulnerability and Debris Flow Risk Assessment.
Three principal components were retained in PCA for the whole debris-prone area regional analysis. Together these first three principal components accounted for 73.03% of the variation in the original 11 variables included in the analysis. The factor loading was dominated by the population and workforce. The loadings of each variable for the retained principal components are detailed in Table 2. As shown in Table 2, the first component was heavily loaded on population and workforce. The second PCA was mainly loaded on transportation related variables. The third component was loaded on production-related variables. Population density, road density, and gross domestic production with the heaviest loads of each component were thus selected for further analysis.
A measurement of vulnerability was calculated by integrating different data and was standardized to a range of 1 to 100 (Figure 5a), where higher values indicated a higher vulnerability towards debris flow. Consistent with relevant studies, the high vulnerability area towards debris flow (vulnerability value >90) comprised the most area because of the relatively low social-economic conditions in mountain areas, which accounted for over 91% of the study area. However, the level of vulnerability could also determine the debris flow risk and management strategies. Moreover, the number of debris flow events had an approximate number by comparing Figure 5c with Figure 5b. Meanwhile, the number of days with debris flows showed a significant difference between the debris flow susceptibility (Figure 4a-c) and debris flow risk (Figure 5b,c) by a paired Wilcoxon signed rank test (p < 0.05), indicating that our indices were suitable for debris flow risk assessment. As shown in Figure 5b-d, the distribution of debris flow risk showed that 53%, 30%, and 98% of debris flow watersheds would result in more than 20 times, indicating that a plan for risk management should take climate change adaptation into consideration in debris flow-prone areas. Meanwhile, the number of days with debris flows showed a significant difference between the debris flow susceptibility (Figure 4a-c) and debris flow risk (Figure 5b-c) by a paired Wilcoxon signed rank test (p < 0.05), indicating that our indices were suitable for debris flow risk assessment. As shown in Figure 5b-d, the distribution of debris flow risk showed that 53%, 30%, and 98% of debris flow watersheds would result in more than 20 times, indicating that a plan for risk management should take climate change adaptation into consideration in debris flow-prone areas. In Figure 6, different vulnerability elements were mainly explained by the GDP, population, and road density, indicating the effects of socio-economic development on debris risk at a watershed scale. Results from the impacts of different vulnerability elements on debris flow risk showed that each selected factor could change the risk level in the study area. Meanwhile, further analysis showed that changes in one of elements were sufficient to cause the reduction in debris flow risk (Figure 5ac). These results suggested that the improvement of socio-economic conditions was pivotal for risk management practices in the development stage of the study area.

Debris Flow Process Models Enhance Susceptibility Evaluation
Compared with the approaches of debris flow risk assessment using statistical models, results from the mechanism model were more accurate in describing higher debris flow susceptibility in potential watersheds. The analysis of climate change potentially increasing the debris flow susceptibility. The differences between the mechanism and statistic model are due to the influence of environmental conditions. Although there was a simple relationship between rainfall intensity and debris flow event [68], other environmental factors, including topography, soil, and land use, are also significant drivers of debris flow [69][70][71]. Different watersheds might offer sufficient conditions to In Figure 6, different vulnerability elements were mainly explained by the GDP, population, and road density, indicating the effects of socio-economic development on debris risk at a watershed scale. Results from the impacts of different vulnerability elements on debris flow risk showed that each selected factor could change the risk level in the study area. Meanwhile, further analysis showed that changes in one of elements were sufficient to cause the reduction in debris flow risk (Figure 5a-c). These results suggested that the improvement of socio-economic conditions was pivotal for risk management practices in the development stage of the study area. Meanwhile, the number of days with debris flows showed a significant difference between the debris flow susceptibility (Figure 4a-c) and debris flow risk (Figure 5b-c) by a paired Wilcoxon signed rank test (p < 0.05), indicating that our indices were suitable for debris flow risk assessment. As shown in Figure 5b-d, the distribution of debris flow risk showed that 53%, 30%, and 98% of debris flow watersheds would result in more than 20 times, indicating that a plan for risk management should take climate change adaptation into consideration in debris flow-prone areas. In Figure 6, different vulnerability elements were mainly explained by the GDP, population, and road density, indicating the effects of socio-economic development on debris risk at a watershed scale. Results from the impacts of different vulnerability elements on debris flow risk showed that each selected factor could change the risk level in the study area. Meanwhile, further analysis showed that changes in one of elements were sufficient to cause the reduction in debris flow risk (Figure 5ac). These results suggested that the improvement of socio-economic conditions was pivotal for risk management practices in the development stage of the study area.

Debris Flow Process Models Enhance Susceptibility Evaluation
Compared with the approaches of debris flow risk assessment using statistical models, results from the mechanism model were more accurate in describing higher debris flow susceptibility in potential watersheds. The analysis of climate change potentially increasing the debris flow susceptibility. The differences between the mechanism and statistic model are due to the influence of environmental conditions. Although there was a simple relationship between rainfall intensity and debris flow event [68], other environmental factors, including topography, soil, and land use, are also significant drivers of debris flow [69][70][71]. Different watersheds might offer sufficient conditions to

Debris Flow Process Models Enhance Susceptibility Evaluation
Compared with the approaches of debris flow risk assessment using statistical models, results from the mechanism model were more accurate in describing higher debris flow susceptibility in potential watersheds. The analysis of climate change potentially increasing the debris flow susceptibility. The differences between the mechanism and statistic model are due to the influence of environmental conditions. Although there was a simple relationship between rainfall intensity and debris flow event [68], other environmental factors, including topography, soil, and land use, are also significant drivers of debris flow [69][70][71]. Different watersheds might offer sufficient conditions to decrease the rainfall threshold, or otherwise increase debris flow susceptibility. Our results showed that the threshold ranged from 20 mm to 110 mm. According to previous studies, a rainfall intensity greater than 30 mm could trigger debris flow in the most debris flow-prone areas [72], which results in a drastic underestimate of the climate change-affected debris flow occurrences to the selection of a higher threshold value.
Rather, our results also suggested that the area of debris flow susceptibility was smaller in the whole debris-flow-prone area. Although many extreme rainfalls would occur in the future, few debris flow watersheds could be identified by the mechanism model. In consequence, a single precipitation threshold for all regions is insufficient to indicate debris flow susceptibility under climate change in risk management implementation. Applying a mechanism model may integrate different influence factors to assess debris flow susceptibility, and could reduce uncertainties arising from climate change projection and the process simulation as well.

Debris Flow Risk Reduction Management in the Future
For our study area, the rainfall threshold of rainfall-triggered debris flows had a wide range. Meanwhile, the variation in the rainfall intensity that triggered debris flow was large enough to take precautions so that differentiated risk management should be implemented in debris flow watersheds to address climate change.
Although a greater than 25 mm rainfall intensity would occur in the whole debris flow-prone area, the distribution of debris flow watersheds was mainly concentrated in the eastern part of the disaster-prone area under climate change, where the number of debris flow has increased continually in recent years [49]. Rainfall intensity, steep slope, and loose materials are pivotal factors for rainfall-triggered debris flows, which are sufficient in our debris flow watersheds. Meanwhile, previous studies suggested that the rainfall thresholds for triggering may decrease after an earthquake [73,74]. For risk reduction management under climate change, these watersheds need more attention, regardless of whether debris flows have been triggered in recent years.
Alongside the continual influence of climate change, the identification of factors of vulnerability in debris flow watersheds lays the foundation for risk assessment and management. Our results found that most of the study area showed a relatively high vulnerability to debris flow and the impacts of elements varied with debris flow watersheds. The debris flow susceptibility was assessed by the mechanism model and key factors of vulnerability were identified in each of the statistically derived watersheds, thus we were able to take different measurements for different watersheds. However, the vulnerability assessment mainly aims to serve risk reduction strategy deployment. In our study, enhancing road network density was found to be the most efficient strategy among the vulnerability factors. Moreover, the effects of the road network density could decrease the risk when the GDP decreased. Even though population density was also decreased in these watersheds, the risk was not increased rapidly. Because road density is beneficial to the circulation of regional development factors, it is regarded as the source of disaster resilience. Additionally, regional GDP also had a significant influence on the debris flow risk of watersheds at which the road density and population were at comparatively low levels. Meanwhile, rural areas have undergone rapid urbanization in China during recent decades, which has also provided opportunities for resilient community construction plans to address climate change and reduce natural hazards [75].
The integration of a debris flow mechanism model and analysis of vulnerability may be more effective for risk assessment at a watershed scale. Compared with others based on threshold values, our results highlighted the higher risk areas. Some studies suggested that moderate debris flows could also damage structures and endanger human lives [76]. Therefore, the selection result of higher risk areas would provide a more conservative strategy for risk reduction management.

Conclusions
In this study, it was found that integrated models which involve a debris flow mechanism model and an analysis of vulnerability may be more effective for risk assessment at a watershed scale. Results from the water-soil process model suggested that the occurrence of debris flows were more than 10 t, indicating higher debris flow susceptibility under climate change in the study area. Spatially, the risk area was mainly concentrated in the eastern part. Meanwhile, many non-debris flow watersheds were identified under climate change. Although debris flow susceptibility is the dominant factor of the increasing risk, strategies to reduce vulnerability cannot be ignored. To mitigate the impacts of debris flows under climate change, any improvement of socio-economic conditions can reduce risk on the mountainous lagging development background, among which the enhancement of road network density may be the most efficient strategy in practice.