Assessing Debris Flow Risk at a Catchment Scale for an Economic Decision Based on the LiDAR DEM and Numerical Simulation

Various risk management measures have been applied to reduce risks associated with the debris flow; however, only a few studies have adopted the economic benefit to evaluate measure effectiveness. The present study sought to explore debris flow risks at a catchment scale and establish the appropriate risk-reducing measures. The Chengbei Gully debris flow in Shanxi province (China) was selected for the case study. High-resolution topographic data of the drainage basin were obtained using the airborne LiDAR technology. FLO-2D software was used to simulate the debris flow process to perform hazard zonation. Vulnerability was estimated based on the location of elements at risk within the hazard zones and the field survey. Several structural and non-structural measures for controlling risks were proposed based on the risk assessment results, and the benefit–cost ratio was used to analyze their effectiveness. The findings indicated that the rainfall event triggering the Chengbei Gully debris flow had an 80-year return period. The total risk under this rainfall condition was 2.3 × 105 $, which was an unacceptable level according to the criteria of tolerance risk. The findings showed that the engineering measure was the best mitigation approach for the Chengbei Gully debris flow with a benefit of 1.35 million $ and a benefit–cost ratio of 6.43.


INTRODUCTION
Debris flows are among the most dangerous mass movements worldwide, especially in the alpine environment (Glade, 2005;Hürlimann et al., 2014;Tiranti et al., 2018;Li et al., 2021). Most debris flows are characterized by large magnitudes, high velocities, and a mixture flow of sediment and water, thus posing potential significant threats to residents and infrastructures (Jakob, 2005;Gregoretti and Fontana, 2008;Gregoretti et al., 2016). Therefore, it is imperative to explore debris flow risks and propose professional mitigation measures based on risk mapping (van Westen et al., 2006;Qing et al., 2020).
The risk of natural hazards can be assessed through qualitative or quantitative methods based on the quantity and quality of datasets. Previously adopted approaches mainly involved expertbased and statistical methods (Guzzetti, 2000;Rossi et al., 2019) which mainly require a long and accurate landslide catalog. Results obtained using these procedures were relatively coarse and could not be compared to other areas in the form of numerical values (e.g., Guzzetti, 2000). Quantitative risk assessment is gradually being conducted at different scales (Dai et al., 2002;Remondo et al., 2005;Erener and Düzgün, 2013), owing to the requirement for risk mapping with a better quality and advances in development of GIS and RS techniques. To estimate the direct risks of natural hazards on elements at risk, studies have adopted statistical theories, such as the magnitude-frequency curve (Hungr et al., 2008;Corominas and Moya, 2010). In addition, previous studies explored quantitative risks based on the event tree analysis (Budetta, 2002;Mineo et al., 2017). In most cases, risk quantification was based on the product of hazard, vulnerability, and the value of elements at risk. Regarding the risk assessment of debris flows, two measurable components are mainly used including the following: 1) the temporal and spatial probability of one event occurrence and 2) the magnitude of harm that the event may result (Chen and Wu, 2016). Some attempts have been made for quantitative risk assessment of debris flows. For instance, Calvo and Savi (2009) proposed a method based on the Monte Carlo procedure for formal risk analysis in debris flow-prone areas. Jaiswal et al. (2011) presented a quantitative procedure for estimating the debris slide risk to life and property and applied it to estimate the risk of potential debris slides in a mountainous area in the Nilgiri hills of southern India. Lin et al. (2011) simulated the flow conditions of typhoon-triggered debris flows and generated risk mapping under different scenarios.
Notably, most case studies regarding the landslide risk assessment only analyzed the initiation area, whereas studies on the risk assessment by incorporating the runout distance are limited (Jaiswal et al., 2011;Guo et al., 2020). However, debris flows are a type of fast mass movements compared with slow-moving landslides; thus, there is a need to perform runout analysis to explore the potential inundation areas and establish actual risks (McDougall, 2017). A few solutions have been developed for this purpose, including constitutive methods and numerical models, based on fundamental assumption simulations (single-phase or multi-phase) (Baum and Godt, 2010;Gan and Zhang, 2019). Numerical modeling is one of the most efficient tools. Relevant approaches such as the DAN3D analysis (Salvatici et al., 2017), smoothed particle hydrodynamic (SPH) method Han et al., 2019), particle flow code (PFC) (Lo et al., 2018), and depth-averaged material point method (DAMPM) (Abe and Konagai, 2016) have received wide applications. The two-dimensional water flood and mudflow simulation program FLO-2D is a useful tool, and it has been successfully applied for the runout analysis of debris flows in several areas worldwide (e.g., Lin et al., 2011;Gomes et al., 2013). FLO-2D can predict important parameters during the movement processes (such as the runout distance, volume of accumulation, flow velocity, and depth), which are critical for estimating the debris flow intensity. Therefore, FLO-2D is an appropriate method for hazard assessment.
Vulnerability assessment is a frequent conceptual and operational challenge in the risk framework (Fuchs et al., 2011;Jaiswal and van Westen, 2013). Although multiple components (e.g., elements at risk and capacity) and various parameters (e.g., information on people and infrastructures) are required for accurate assessment, most cases lack sufficient data. This explains why the methods for vulnerability assessment are generally qualitative or semi-quantitative. Previous studies developed quantitative models; however, most are used for buildings exposed to natural hazards and are primarily based on physical tests or empirical theory. For example, Uzielli et al. (2008), Li et al. (2010), and Peduto et al. (2017) explored the vulnerability estimation function or the curve for the scenariobased landslide hazards. A unified quantitative method for vulnerability estimation is yet to be proposed; thus, a semiquantitative method was adopted in the present study.
The final result of risk assessment is mainly expressed by the expected losses resulting from a hazardous event of a given intensity. Studies thus use different indicators to show the risk distribution of hillslopes subjected to natural hazards, including the non-dimensional risk level (Bell and Glade, 2004;Arksey and VanDine, 2008), direct economic damage (Vranken et al., 2013), indirect damage (Jaiswal et al., 2010), and the total number of casualties (Guzzetti et al., 2005). Moreover, some studies applied the index of annual risk to assess the risk within 1 year; thus, the results for different cases and scales can be compared directly (Chen and Wu, 2016). In addition to mapping the risk zonation, the obtained risk results can be used for further objectives. For instance, parameters of existing models can be calibrated using the results of past events, to make them more accurate in predicting the risk under future scenarios (Bertolo and Wieczorek, 2005). Some studies used risk results as a reminder to the local authorities and residents to understand the potential threats (Graff, 2014). The quantitative risk assessment mainly seeks to provide the basis for the design and implementation of mitigating measures (Holcombe et al., 2012;Vranken et al., 2013). Although several measures can effectively reduce the risk, their time or economic cost is unacceptable for organizations and individuals faced with the risks. Studies on risk management show that risk mitigation measures are not primarily geared towards risk elimination but to maintain the risk to an acceptable level. Therefore, it is necessary to evaluate and compare the social or economic effects of the potential measures to select suitable measures (Peila and Guardini, 2008). Unfortunately, the present case studies on this topic (e.g., Chen and Wu, 2016;Guo et al., 2020) are limited and should be updated, especially regarding fast mass movements.
The current study sought to perform a test to provide a basis for making the decision of risk mitigation measures through cost-benefit analysis. The specific objectives include the following: 1) Simulation of the runout process of the debris flow using remote sensing DEM data and FLO-2D software; 2) Vulnerability estimation based on field survey data and the semiquantitative method; 3) Application of the procedure to the Chengbei Gully of Shanxi province (China) (approximately 1 km 2 ), to calculate the risk level in a given rainfall scenario; and 4) Proposing several measures as potential choices for risk reduction and comparing their effectiveness from an economic point of view.

Study Area and the Debris Flow Event
The Chengbei Gully is located in the Lvliang mountainous area to the east of Loess plateau (China), a well-known geohazard-prone region (Juang et al., 2019). The Chengbei Gully administratively belongs to Ji County in Shanxi province, and its coordinates are E110°40′25″, N36°6′14″ (Figure 1). A heavy rainfall event with an accumulated value above 100 mm and a maximum hourly rainfall of 54.6 mm was recorded in this area on July 3rd of 2013 . After the rains, rain water mixed with soil particles, rubbles, and stones flowed along the gully downstream and destroyed courtyards and residents' houses (Figures 2A,B). Additionally, the surface land cover was also affected much. As seen in Figure 2C, the remote sensing image showed the impact of the event on the land cover in the channel up to April 2014.
As defined by Hungr et al. (2014), the debris flow is very rapid to the extremely rapid surging flow of saturated debris in a steep channel and strong entrainment of material and water from the flow path. Meanwhile, the rainfall and intensity in the study area are similar with some other debris flow events reported in China (260 mm for total rainfall and 60 mm for maximum daily rainfall in Wei et al., 2018). Hence, the event observed on the Chengbei Gully on July 3rd of 2013 was considered as a debris flow triggered by an extreme rainfall.
The drainage area at the upper stream of the basin covers approximately 3.39 km 2 and is 6.76 km long. The average slope of the initiation area at the channel is approximately 30°with an elevation difference of 289 m. Geomorphologically, the entire gully can be categorized into three sections ( Figure 3). The upper part has two branch channels with the transverse sections exhibiting a typical "V" shape, with steep slopes at both sides (~45°). A channel with a gentle valley bottom and relatively narrow width occurs at the middle section. The middle section is 1,600 m long. The lower section is the widest and gentlest part of the gully with a relatively high discharge, and water flows into the Qingshui River which is a branch of the Yellow River. The geological setting of the Chengbei Gully mainly comprises Quaternary loess deposits with an average thickness of 100~200 m. Notably, silty clay and rubbles loosely accumulate at the toe of slopes. These materials provide good sources for debris flows. A field investigation reported that the total volume of these materials was approximately 1 × 10 5 m 3 . Studies show that no humans have settled in the upper and middle stream areas; however, there is a small village with a few residential houses in the lower stream.

Remote Sensing Data
Light detection and ranging (LiDAR) sensors can generate highresolution three-dimensional (3D) coordinates of objects. Height information (that is Z values) can be used to characterize the topography of an area (Wu et al., 2021). The Chengbei Gully is not a large-scale area; therefore, a high-resolution digital elevation model (DEM) is important for subsequent simulations. In the current study, the manned airborne LiDAR ( Figure 4A) was used to obtain the DEM data.
LiDAR mapping was conducted on July 2017 using the biplane DJI F550 purchased from Da-Jiang Innovations Science and Technology Co., Ltd. (DJI), Shenzhen, China. The biplane comprised laser-scan systems in near-infrared sensors. The load capacity of the biplane was 3 kg, and the flight lasted for approximately 40 min. The overflight was planned and operated by a professional employee of the China Geological Survey (CGS). The path was roughly rectangular in shape with a length of 4 km and a width of 2 km. The flight altitude and velocity were fixed to ensure a relatively consistent overlap between images. Retrieved LiDAR data were processed by Terrasolid software (Terrasolid Ltd., Helsinki, Finland). Point clouds were generated in a LAS format and subsequently converted to DEM data through ArcGIS software with a resolution of 1 m ( Figures 4B,C). As seen in Figure 1C, two control points were used to evaluate the data quality. The accurate elevations of the two control points obtained from the contour line were 908 m (point #1) and 842 m (point #2) asl, respectively. Their elevations revealed in the LiDAR DEM were 908.2 and 842.1 m asl, respectively.

In Situ Data and Laboratory Tests
The field survey was conducted between July and October 2017. Targets of the survey mainly included 1) investigating and recording deformation and movement of materials during the debris flow by engaging with residents; 2) obtaining information on the vulnerability of the elements at risk (population, Frontiers in Earth Science | www.frontiersin.org February 2022 | Volume 10 | Article 821735 constructions, and local price level); 3) exploring information on geological settings (e.g., topography, materials, and geological structures), to have a deeper insight on risk reduction; and 4) collecting soil samples to perform laboratory tests to determine the physical parameters for numerical simulation.
The following laboratory tests of soil samples were performed in the Xi'an Center of China Geological Survey: 1) the soil unit weight test; 2) the strain-controlled test to determine soil viscosity; and 3) the soil bulk concentration test. All these followed the national criteria in China (MOHURD, 2019).

METHODOLOGIES
The general procedure of this study includes the following: 1) risk calculation, 2) analysis of the risk level based on the tolerance criteria, and 3) specific designs of risk mitigation measures and effectiveness analysis. Figure 5 outlines the research design. The principles for the methods used in the current study are described in the following sections.

FLO-2D Simulation
The movement process of the Chengbei Gully debris flow was simulated by FLO-2D software, a two-dimension numerical model developed by O'Brien et al. (1993). The software adopts the quadratic rheological model which allows simulations for various conditions ranging from clear water to sediment flows. The FLO-2D modeling is controlled by the terrain and is based on the continuity equation and dynamic wave equation (Lin et al., 2011;Gomes et al., 2013).
In this study, DEM data were used as the input for the FLO-2D BASE environment to obtain topographic attributes and to develop the geological model. Rheological parameters defined by the FLO-2D manual (O'Brien, 2007) were determined by laboratory tests (section 2.2) and empirical values. In order to reveal the inherent uncertainties in these parameters, a potential range for these parameters was also provided (O'Brien, 2007;Gomes et al., 2013;Haas and Woerkom, 2016). The range can also be used to ensure that the future structural countermeasure for risk reduction would work as planned.
As seen in Figure 1C, the junction of two branches in the upper stream area was used as the inflow point for the catchment model. The grid with the lowest position in the watershed served as the outflow point, where the water flows confluence with the water from another gully. Data on rainfall events, including cumulative rainfall and maximum hourly rainfall, were also required. The grid size was set to 2 m because the simulation duration under this condition was about 3 h, which was similar to the duration of the event. The parameters used are presented in Table 1. These parameters are static, and the software uses them to analyze the flowing process. Some key indicators in the flowing process are the output results of it, including the flowing depth and velocity, which we will describe in following parts.
The Manning's roughness coefficient was set at 0.2 because forest was the main land use in this area (see Figure 3). According to Brater and King (1976) and Montes (1998), this type of land use has a significantly higher value of n. The value used in this study was higher than that used in most cases [e.g., Liu et al. (2013); Haas and Woerkom (2016)] but close to the value reported by Cui et al. (2011).

Hazard Mapping Based on Runout Analysis
A hazard can be defined as the degree of debris flow damage within a given period in a specific area. It is mainly associated with various parameters, including the probability of occurrence  and intensity (; Corominas et al., 2014). In this study, the probability of occurrence of the debris flow under the given rainfall condition was 1; therefore, the main task was to determine the intensity. In practical application, qualitative methods are mainly used to assess the intensity, based on two crucial parameters (flow depth and flow velocity) in the movement process of debris flows. A three-degree hazard matrix was defined for the Chengbei Gully debris flow based on a few examples of reclassification of the intensity proposed in previous studies (Hürlimann et al., 2006;Tang and Okimura, 2006;Chen et al., 2010) (Figure 6). Most buildings in the area had a floor height of~2.5 m; therefore, this value was used as a threshold of flow depth that can significantly damage the buildings. At a flow depth lower than 0.5 m, people can breathe even if they are caught by the deposits; thus, this was considered as a relatively safe depth. However, it should be noted that some studies consider that breathability cannot be used as a measure in a debris flow disaster, so the depth of 0.5 m may also be fatal (e.g., Raetzo et al., 2002;Hürlimann et al., 2006). Most residents living on the slope were adults aged between 40 and 60; thus, they can escape from debris flow with a velocity of 0.5 m/s. However, if the flow velocity is faster than 1 m/s, ordinary people can hardly escape.

Calculation of the Return Period
Every risk scenario can be attributed to triggering events with a specific return period (Zêzere et al., 2008). Therefore, it is imperative to determine the return period when assessing the temporal probability of the risk or predicting and modeling different future scenarios (e.g., Guo et al., 2020). An extreme value type I (i.e., Gumbel distribution) method was used to fit the rainfall events to explore the rainfall frequency. This statistical approach is widely applied in analysis and prediction of hydrological events, including flooding events (Onen and Bagatur, 2017). The cumulative distribution function of Gumbel distribution is expressed as shown below (Gumbel, 1941;Matti et al., 2016): where X is a random variable, x is a possible value of X, ξ and α are the location parameter and scale parameter, respectively, which are calculated as follows (Gumbel, 1941;Matti et al., 2016): where μ X and σ 2 X are the mean and variance of the dataset, respectively. To explore changes in rainfall events over the entire recording period, the Gumbel distribution, fitting to daily rainfall peaks for every rainy season in 50 years (from 1964 to 2013), was performed using the MATLAB tool.

Vulnerability Estimation
Two types of elements at risk, including property and life, are considered when determining the debris flow . Vulnerability of these elements varies with the hazard intensity and resistance to the debris flow and requires several indicators for evaluation. However, total quantitative vulnerability estimation is not performed in most cases due to lack of sufficient damage data (Fell et al., 2008;Jaiswal et al., 2011). In this study, a semi-quantitative estimation based on hazard mapping and information from the fieldwork was used to determine the vulnerability. The hazard intensity was evaluated using the acquired information mainly based on the depth and the velocity of the debris flow reported in section 3.1. In addition, information from the fieldwork was used to explore the resistance of elements at risk. For constructions, the impact force was the direct indicator associated with the hazard intensity. High values of flow depth and velocity imply a larger impact force of the debris flow, resulting in higher vulnerability of the building. The main factors contributing to resistance of constructions include height and materials. Hence, the relationship between the building height and flow depth is important for vulnerability estimation. For estimation of vulnerability of people, the model must account for indicators such as the age, healthy state, and the implementation of early warning systems within the community (Du et al., 2016;Guo et al., 2020).

Analysis of the Risk Level Based on Risk Tolerance Criteria
Risk assessment determines whether the present risk level is tolerable to the society or existing mitigation measures are adequate. There is a need to evaluate whether alternative control measures can be implemented if the risk is not tolerable and if mitigation measure are not effective. Life risk is the main indicator of the social risk tolerance criteria, for analysis of the number of people threatened in a disaster and the probability of death. Although some risk tolerance criteria have been proposed and implemented (e.g., Bottelberghs, 2000;AGS, 2007), the findings are controversial. Therefore, the criterion developed by Hong Kong Geotechnical Engineering Office (1998) was adopted in the present study. The criterion determines the tolerable social risk based on the F-N curve, where N represents the number of fatalities and F represents the frequency of N or more fatalities per year. The curve ( Figure 7) comprises four parts with different risk levels, including the acceptable risk, unacceptable risk, intense scrutiny region, and the "as low as reasonably practicable" (ALARP) region. Hence, if the risk level of a given scenario is unacceptable or is at the ALARP level, measures to reduce and mitigate risks must be considered. This implies that a measure is considered reasonable when it reduces the risk from the unacceptable or ALARP level to an acceptable level.

Risk Calculation and Mapping
The risk can be calculated using the formula below (Fell et al., 2008;Corominas et al., 2014): reaching a specific point. If an element at risk is located at the inundation area of the debris flow, this value is 1. P (S:T) is the temporal-spatial probability of the elements being at a certain point when the debris flow occurs. The P (S:T) for the constructions at risk in the inundation area obtained from the runout analysis is 1. Analysis of the people residing at the area showed that their time at home was approximately 10 h per day; therefore, the P (S:T) was set at 0.42. E is the economic value of the elements at risk or the number of people. Values for E were assigned to every property based on their type (constructions or indoor possessions), local price level (The people's Government of Shanxi Province, 2020), and field survey. The simple classification method used and the specific values obtained are shown in Table 2. This study was carried out for economic decision-making; thus, the values for life (E i ) were estimated using the method described by Shang (2012) as follows: where g represents the average gross domestic product per person per year, e denotes the average expected lifetime of people, and w denotes the ratio of working time to leisure time. The National Bureau of Statistics of China (2020) reports that g is 10,280 $/year, e is 77 years old, and w is approximately 0.2 in villages. Therefore, the total economic risk posed by the Chengbei Gully can be expressed as follows: where R T is the total economic risk of the Chengbei Gully and R p and R L represent the property risk and population risk, respectively.

Principles for Implementation of Measures
Potential measures considered for risk mitigation of the debris flow generally can be categorized into structural and nonstructural approaches. The structural approach mainly decreases the hazard intensity (or failure probability), whereas the non-structural approach improves vulnerability (or the resilience capacity) of elements at risk. Structural measures comprise two aspects: 1) Reducing the probability of the debris flow occurrence. This can be achieved by reducing the amount of sediment source to relieve the occurrence of debris flows. 2) Reducing the probability of the debris flow reaching the elements at risk. Considering the aforementioned two aspects , structural measures may include constructing blocking dams in  the upper or middle stream, drainage channel in the lower stream, and fence walls at the back of houses. Non-structural measures include the following options: 1) Implementation of early warning systems, which involves professional monitoring, real-time warning, and a public broadcasting system. These measures can help residents to timely evacuate from dangerous places, thus decreasing the spatial-temporal probability of elements at risk. 2) Community-based public education can improve disaster awareness and the response ability of people to reduce vulnerability of the elements at risk. 3) Relocation of residents in new settlements can reduce the number of the elements at risk.
In this study, considering the available information before the debris flow occurrence and the potential space that can be used for structures, building blocking dams in the mainstream of the channel is a suitable structural measure. For the non-structural measures, early warning systems, community-based public education, and resident relocation are feasible.
Hence, four single categories of measures were selected for the risk mitigation of the Chengbei Gully debris flow. Moreover, community-based public education, which is a non-engineering measure, is often implemented together with engineering measures or early warning systems. Therefore, the two measures, namely, community-based public education + early warning system and community-based public education + dam were also included.
A detailed design on these measures is described in the section below.
where Q S is the total volume of solid materials, Q is the volume of all materials, γ c is the unit weight of the debris flow, γ S is the unit weight of solid materials, and γ w is the unit weight of water. Q can be obtained as follows: where T is the continual time of the debris flow, and K is an empirical coefficient which is proposed as 0.202 in the guidebook. Q c is the peak quantity discharge of the debris flow and is expressed by where φ is the amendment coefficient for the unit weight of the debris flow, Q p is the peak quantity of the water flow with a given return period, and D c is the blockage coefficient of the debris flow which is set to 1.5. φ is calculated as follows: Geometrical characteristics of the dams designed based on Eqs 8-10 are presented in Figure 8. The lengths of dam 1# and dam 2# were 25 and 45 m, respectively. The dam type was selected as the rubble gravity dam, owing to the low cost and easily acquired materials.

Early Warning System
The main components in the early warning system include the following: 1) Monitoring of rainfall. Three rain gauges were installed in the channel, two in branches upstream, and one downstream. The gauges provide an alert in case of heavy rainfall exceeding the given threshold (Wei et al., 2017). 2) Monitoring of infrasonic signals. Low-frequency infrasonic waves emitted by movement of the debris flow can be detected over a kilometer radius (Liu et al., 2015). Hence, installation of a monitoring device in the upper stream can be used for remote monitoring. 3) Monitoring of maximum depth. A set of wires stretched and fixed across the channel can continuously measure the distance between the detector and the flow surface. If the maximum depth of the debris flow increases to a detectable level, then the alert is made. 4) Video monitoring. Three cameras were installed in the upper, middle, and lower stream separately to record and monitor the movement of flows 24 h per day.

Community-Based Public Education
This measure has been validated and is effective in some natural hazard-prone areas of China (Liu et al., 2016;Guo et al., 2019;Guo et al., 2020). Hence, it was also considered alone in this study. This measure mainly includes 1) posters which present and advertise the benefit of risk mitigation, 2) brochures for residents, 3) classes on disaster prevention to young students in school, and 4) professional lectures to adults given by experts in this field.

Resident Relocation
This measure involves relocation of residents from the inundation area of the debris flow. The local authority is responsible for the relocation and construction of new settlements. The new settlement area should accommodate the relocated families as per the principle of average distribution (Guo et al., 2020).
Specific implementation and locations of all potential measures are shown in Figure 9.

Economic Cost Analysis
The total cost of each measure (T C1 ) must be evaluated in the first step as shown below: where T 1 , T 2 , and T n represent the total cost of the 1 st , 2 nd , and n th item, respectively. P 1 , P 2 , and P n represent the unit price of the 1 st , 2 nd , and n th item, respectively, whereas N 1 , N 2 , and N n represent the number needed for the 1 st , 2 nd , and n th item, respectively. All prices are estimated according to the price level as of 2015~2020 in Shanxi province (National Bureau of Statistics of China, 2020).   All the items and their unit prices included in every measure are presented in Table 3. Although most of these costs are paid for at present, a few items will be acquired in the future (e.g., operation cost of the early warning system and public education); therefore, it is necessary to convert them into present economic values (T PC ) as following: where C 1 , C 2 , . . . C n are the single-year cost for future 1, 2, . . . , n years, n is the life cycle of the measure, and i is the annual interest rate. Hence, the total present economic cost (T C ) for a specific measure can be denoted as follows:

Effectiveness Analysis and Comparison
Effectiveness of measures is a crucial decision-making index because it can help determine the measures that reduce the risk at the largest level at a relatively small cost. Reduced total risk achieved by mitigation measures can be considered beneficial since effectiveness of the measure is expressed as the benefit-cost ratio (Fuchs et al., 2007;Jaiswal and van Westen, 2013). A higher ratio indicates a more effective measure. The structural measures have a life cycle, and their discounting is applied to convert the future benefit to present economic values. Therefore, different effectiveness of measures can be compared. The total present value of the benefit of a specific measure (T B ) can be calculated as follows (Chen and Wu, 2016): where B 1 , B 2 , . . . B n are the single-year benefit for future 1, 2, . . . , n years, n is the life cycle of the measure, and i is the annual interest rate. Hence, the effectiveness of the measure (Ef) is denoted as following:

Hazard Mapping and Vulnerability Estimation
Simulation results obtained from FLO-2D, including the flow depth and flow velocity maps, are shown in Figure 10. The results showed that the area inundated by the debris flow under the rainfall condition in this study was approximately 1.6 × 10 5 m 2 . The area with the flow depth less than 1 m accounted for 40% of the inundation area, with the area of depth from 1 to 3 m accounting for 48%, and the area of depth more than 3 m accounted for only 12% of the area inundated by the debris flow. Two places with relatively large flow depths were identified. The first place was at the boundary between middle and down streams with 6.73 m as the maximum flow depth. This area was characterized by a wider and gentler valley which provided a favorable terrain for the deposition. The second place was at the mouth of the gulley with a flow depth of exceeding 4 m. Moreover, an accumulated alluvial fan took this place as the center covering an area of 4 × 10 4 m 2 . Analysis of the flow velocity showed that its distribution presented a similar characteristic to the flow depth. As a whole, the middle of the channel had larger velocities than that of its edge whereas the overall velocity in the upstream and middle stream was higher than the velocity downstream. This finding was associated with the terrain. For instance, the gully channel was narrow and relatively steep in the upper part but wider at the lower part. The average flow velocity of the whole gully was approximately 1.2 m/s, whereas the maximum velocity was approximately 4.0 m/s and was close to the boundary of middle and down streams. The hazard map ( Figure 10C) of the debris flow was obtained based on the criteria outlined in Figure 6. High-level hazard areas were mostly located in the middle of the channel in the upstream and middle stream. Moreover, the mouth of the gully was another important highlevel hazard area. Areas of high, moderate, and low hazard levels accounted for 30, 40,, and 30% of the total area, respectively. It should be noted that there may be a certain relationship between the flow depth and velocity, but it should depend on the actual terrain that the debris flow occurs (Shu et al., 2018). In this study, the overall distribution of the flow depth and velocity in the channel is similar, but some specific characteristics vary. For instance, at the bottom of the gully, the depth in the middle of the channel is larger than surrounding places, but the velocity values of these places do not have much difference. Additionally, the upper section of the channel has continuous relative large depths while the velocity of this part has values with different levels.
Settlement of residents was distributed on the lower part of the gully; thus, the current study only described the vulnerability of this area (Figure 11). Previous studies (Fuchs et al., 2007;Jakob et al., 2012;Mavrouli et al., 2014) mainly expressed vulnerability as a function of either flow velocity or deposit depth. However, vulnerability significantly increases at a deposit thickness over 2 m or the velocity above 4 m/s. The vulnerability value of a property in the high-vulnerability area was set as 0.9, whereas the values for middle-and low-vulnerability areas were 0.6 and 0.1, respectively, based on the processes of hazard zonation. The findings on the vulnerability of people showed that most of the people in the settlements were adults between 40 and 60 years (the middle level on physical function). In addition, some aspects which can reduce vulnerability were lacking, such as implementation of early warning systems, training on emergency responses, and public education. Hence, the vulnerability for the people in the area with a high-level hazard was set to 0.6, whereas the vulnerability for the other two areas was set to 0.3 and 0.1, respectively. Vulnerability of the current study was slightly higher than those reported by Fuchs et al. (2007) and Guo et al. (2020); however, it was lower than the value reported by Solari et al. (2020). Analysis of the distribution of the elements at risk showed that the total area of settlements in high-, medium-, and low-hazard intensities was 3,775 m 2 , 10,460 m 2 , and 11,070 m 2 , respectively. The total number of residents living in the three areas was 24 (high-hazard intensity), 45 (middle-hazard intensity), and 60 (low-hazard intensity), respectively.

Risk Assessment Results
The rainfall data in the study area are shown in Figure 12. The results showed that the triggering rainfall of the Chengbei Gully  (Table 4) was obtained according to Eq. 4. The total economic risk for the property and population was 8.9 × 10 4 $ and 1.4 × 10 5 $, respectively. Analysis of spatial distribution showed that the risk distribution was similar to that of the hazard distribution. This is because the elements at risk with the same hazard intensity were assigned the same vulnerability values. Although the middle-hazard intensity area had the largest total property risk, the average risk for each house was medium due to a larger total area of elements at risk. Notably, some areas with a high hazard had low-risk levels owing to very few elements at risk located mainly at the boundary between middle and downstreams. Hazard and risk levels at the mouth of the gully were high because it was close to the community settlements. The risk levels in the areas with three hazard intensities (low, moderate, and high) were then calculated according to Eq. 4. If 5% of the total population was considered as the number of fatalities (Petley, 2012;Guo et al., 2020), the corresponding frequency was calculated as historical data on the number of fatalities caused by the debris flow were not available. As seen in Figure 13, the risk levels of all three areas were unacceptable when evaluated by the risk tolerance criteria, thus indicating the risk mitigation measures must be considered and implemented to reduce the total risk of the Chengbei Gully.

Effectiveness Analysis and Comparison
Analysis of the cost measure showed that the cost on dam construction and the overall relocation were incurred only once whereas other measures (early warning system and public education) comprised the future investment. The life cycle of public education was approximately 5 years (Chen et al., 2010). The life cycle of the early warning system was estimated to 20 years owing to its association with the service life of the piece of equipment. The annual interest rate i ranges from 0 to 0.12 (Staats, 1969) and was estimated as 0.05 for the Chengbei Gully. This value was higher than the value reported in Europe (e.g., i = 0.016 in Ferlisi et al., 2021) and slightly smaller than the value used in Taiwan, China (e.g., i = 0.06 in Chen and Wu, 2016).
The single-year benefit (B in Eq. 14) can be estimated by the annual benefit rate. The annual benefit rate of each measure was determined by the semi-quantitative estimation method. Previous studies (e.g., 40% in Salbego, 2014;30% in Salbego et al., 2015;40% in Guo et al., 2020) report that the benefit of structural measures accounts for approximately 30~40% of the  total cost; thus, the annual benefit rate of the dam construction was determined as 30%. The early warning system and public education are under the scope of community-based disaster reduction, so both were considered to have similar effectiveness for risk mitigation. China began implementing these two measures from 2005. Based on the saved lives by these measures and number of deaths caused by geohazards since 2005, the annual benefit rate of these measures has been estimated as 2% (Guo et al., 2020). Because the life cycle of the early warning system was designed as 20 years, its annual benefit rate from the 21st year was set as 1%. Relocation removes the elements at risk from the submerged area of the debris flow, thus eliminating the risk (the benefit rate is 100%). Therefore, the effectiveness of each measure can be obtained, as shown in Table 5, according to Eqs 11-15, Table 2, and Table 3. Analysis of the four single measures showed that the total cost ranged from 7.37 × 10 4 $ to 2.53 × 10 6 $. Relocation had the largest budget, followed by construction of the dam, the early warning system, and public education. Regarding the total benefit, relocation is the highest, higher than dam + public education, followed by the dam. The other measures have lower benefits than these three measures. Analysis of the effectiveness showed that the benefit-cost ratios of all measures were greater than 1. The total benefit was higher than the total cost, thus indicating that all these measures were cost-effective. Specifically, the dam construction had the best effectiveness with the benefit-cost ratio of 6.43, followed by dam + public education (the benefit-cost ratio was 5.07). The benefit-cost ratio of these two measures was evidently higher than the other measures. The relocation had the highest total cost and total benefit at the same time, so its effectiveness was not outstanding. To sum up all the aspects, the dam and dam + public education had good benefits (ranking second and third, respectively) and the highest benefit-cost ratio. Hence, these two measures can be the best options for risk reduction in this study.

DISCUSSION
Definition of the debris flow intensity is rather controversial owing to the several factors that have to be considered. However, FIGURE 13 | Calculation of the risk level for the Chengbei Gully debris flow using the risk tolerance criteria. sufficient data are not available in the case of the Chengbei Gully. Flow velocity is perhaps an appropriate choice because it can significantly represent the impact forces caused by debris flows.
The flow depth was also considered in the current study in addition to flow velocity. Hence, the three-degree intensity matrix was defined by combining these two indicators. Given that there is no agreement on this topic, we determined the thresholds between every two degrees according to the literature data, characteristics of elements at risk, and our experiences obtained from the mountainous areas of China. The final intensity matrix was similar to matrices used in engineering practices of other mountainous areas in China and can be considered for comparison with other criteria developed in other countries (e.g., Hürlimann et al., 2006). However, the intensity matrix should be combined with the probability of occurrence to obtain the final hazard matrix. In the current study, the triggering condition has been fixed, and the probability of occurrence for the debris flow is 1 under this condition. The matrix can be replicated in other debris flow catchments, owing to its simplicity and relatively low time cost of the intensity matrix. Furthermore, the thresholds in the matrix can be changed to adapt to new situations. Flow velocity and flow depth within a specific point were determined using the FLO-2D algorithm, which could also simulate initiation of the debris flow and identify the inundation areas. This enables analysis of the risk in these areas. However, simulation depends significantly on the quality of input data. In this study, the digital terrain model used the data with a high spatial resolution obtained from airborne LiDAR whereas rheological parameters were mainly obtained from the FLO-2D manual. Hence, inaccuracies of results would be mainly attributed to uncertainties of the rheological parameters owing to lack of the parameter calibration. However, according to Table 1, the potential range of parameters can be roughly determined from the laboratory tests and the literature, which can be used to ensure that the structural countermeasure would work as planned. On the other hand, the robustness of the mitigation measure with respect to the uncertainty linked with the debris flow rheology can be estimated from the parameter range. Although the estimation for the robustness is not quantitative, it can provide an insight for the reliability of the measures.
Additionally, poor validation of the FLO-2D model was a limitation of this study because calculation of flow velocity and flow depth was the basis for assessing debris flow risk. Although some examinations on the morphological aspect in the field after the event revealed some similarities between actual and simulation values, it does not effectively prove the model accuracy. However, it has been clarified in the literature that uncertainties introduced by the topography may be much larger than uncertainties in other information (e.g., Luppichini et al., 2019). The uncertainty level from the modeling process is, thus, acceptable in this analysis. Notably, the past event can provide useful data and experience for future studies in the parameter calibration and accuracy validation, which will benefit the risk prediction of a given scenario (Cai et al., 2021).
Uncertainty analysis is an important issue in quantitative risk estimation, given its inherent existence in every input factor (components). A preliminary analysis on uncertainties in this study included the uncertainty level, the reason for uncertainties, and potential measures to improve the issue ( Table 6). The highest uncertainty in the procedure exists in the vulnerability estimation, as reported in many previous studies (e.g., Jaiswal and van Westen, 2013;Guo et al., 2020). A semi-quantitative model was used to analyze the resistance of elements at risk against the disaster which was close to an empirical method because a high subjectivity was included. A potential solution for this issue is to develop physical models based on vulnerability curves. Some efforts have been made on this topic (e.g., Peduto et al., 2017), but most of them are with regard to vulnerability of constructions and should be validated by a series of laboratory tests. Other studies performed vulnerability calculation using various indicators such as the structure type, foundation depth, and maintenance condition. However, this type of model requires detailed field survey and sufficient information on elements at risk. A high uncertainty level on the risk results implies that risk mapping should be treated as a relative result and not absolute. Similarly, in this study, the results are not refuted but highlighted the importance of uncertainty analysis. As Bell and Glade (2004) reported, although some improvements can be performed, uncertainty in risk analyses will always occur. Therefore, honestly informing the user of the associated limitations is important, instead of ignoring them. The present procedure provides a relatively simple approach to handle risk assessment, especially in the case of insufficient data. This modified procedure presents a relative risk level of elements at risk to users and the local authorities. Furthermore, resultant mapping is a fundamental task for economic evaluation for all potential measures of risk mitigation. The main objective of the current study was to compare and determine risk mitigation measures; thus, this present procedure of risk assessment is appropriate from the economic viewpoint.
Uncertainty also exists in the analysis of the risk reduction measures. For example, the dam was set to intercept all the solid material; however, the total solid material is hard to be accurately determined. Besides, if the dam was used to decline the peak discharge once the flow would not threat citizens in the downstream, the dam could be much lower, and the cost could also decrease. However, it should be noted that this kind of uncertainty is related with engineering practices, which is acceptable for users. Moreover, the cost and benefit will change with the dam design at the same time, so the final effectiveness will not be affected much.
The variation of risk values of a debris flow risk map within a given area under a particular rainfall condition is influenced by intensity, spatial probability reaching a specific point, vulnerability of elements at risk, and temporal probability. The findings from the current study indicate that if mitigation measures are not considered from the beginning, the risks posed by the debris flow to the constructions and lives may be significantly high. Notably, the risk of the Chengbei Gully debris flow was at an unacceptable level, as shown by the tolerance criteria. Risk mitigation measures can, thus, be implemented Frontiers in Earth Science | www.frontiersin.org February 2022 | Volume 10 | Article 821735 through reduction of the values of the aforementioned four aspects. In the current study, potential measures were mainly taken based on three components in the risk, namely, temporal and spatial probability of the debris flow reaching a specific point and vulnerability of elements at risk. The subsequent analysis showed effectiveness of these measures; Here, the proposed measures decreased the risk by 1/3 of its original risk.
Although the absolute value of the risk reduction was not ascertained, the risk posed by the investigated event was significantly decreased. The findings imply that elements at risk can be well protected. The aforementioned results showed that it is difficult to eliminate the risk, but mitigation measures can reduce the unacceptable risk to an acceptable level to avoid significant losses. However, local authorities must be aware that new risks still may be caused by future extreme scenarios; therefore, continuous and systematic solutions must be included in the long-term plan of risk management and control. Overall, the risk assessment of the debris flow is a multidisciplinary procedure where quantitative analysis and qualitative analysis inherently coexist (Glade, 2005;Hürlimann et al., 2006). However, the focus of this study is on the global understanding of the present approach and the experiences acquired during the case of the Chengbei Gully debris flow. Hence, the abovementioned exhaustive discussion already presented the uncertainties in this study and can help readers to better understand the rigor of the article.

CONCLUSION
In this study, numerical modeling based on FLO-2D software and the semi-quantitative vulnerability estimation was performed to assess the debris flow risk in the Chengbei Gully of central China. Remote sensing and in situ field data were used to derive a high-resolution DEM and to obtain parameters for FLO-2D modeling. Key indicators during the movement process including the runout distance, flow velocity, and flow depth were simulated to determine the hazard map based on the debris flow event in July 2013 in this region. The findings showed that the rainfall event triggering the Chengbei Gully debris flow had a return period of 80 years. Under this rainfall condition, the total risk of the Chengbei Gully debris flow was 2.3 × 10 5 $. Notably, the risk level of the lower part of the gully was the highest due to the existing community settlement. The risk level in the study area was unacceptable since the risk level obtained from the risk tolerance criterion method was greater than the standard curve.
Four single mitigation measures and two mixed measures were considered and designed for risk reduction, including construction of the masonry gravity dam, implementation of the early warning system, community-based public education, and overall relocation of the residents. Cost-benefit analysis for these measures showed that the economic benefits varied from 7.37 × 10 4 to 4.51 × 10 6 $, whereas the benefit-cost ratio varied from 1.03 to 6.43. Benefit-cost ratios of all measures were above 1, implying that the current measures had larger benefits than the costs and were cost-effective. Although the total cost of the dam construction was not the smallest, it had the greatest economic benefit-cost ratio; thus, it was a better choice for risk reduction of the Chengbei Gully.
In summary, the present study indicates that effectiveness analysis for the risk mitigation measures should evaluate the absolute value of benefit or risk reduction and should also focus on the benefit-cost ratio obtained. In addition, it is difficult to eliminate the risk from debris flow events completely if the mitigation measures are not designed and performed on time. Therefore, measures to control the total risk to an acceptable level are important. The present modeling procedure allows users to implement the quantitative risk assessment during a short time as a basis for analysis of potential mitigation measures and economic decision. The quantity and quality of datasets in this study are at the normal level; thus, the entire process can be easily replicated and applied in other areas with similar settings. However, the uncertainty level of the obtained risk results may be significant, especially the estimation of vulnerability and damage.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
YT was responsible for data curation, methodology, software, and writing of the original draft. ZG contributed to the conceptualization, methodology, supervision, and writing of the original draft. LW carried out the methodology, software, visualization, and validation. BH assisted with software, visualization, and reviewing and editing of the manuscript. WF helped with the visualization and supervision. XS carried out the methodology and reviewing and editing of the manuscript. ZL was responsible for methodology and software. YZ helped with reviewing and editing of the manuscript.