Yongji county is located in central eastern Jilin Province, China (Figure 1), which covers a total area of 2,620 km². The number of debris flows in Yongji county has increased from 71 in 2007 to 129 in 2021, causing several deaths, destroying hundreds of houses and thousands of acres of farmland. The debris flows scoured the roadbed and piled up on the road, resulting in traffic paralysis. It is necessary and urgent to map the susceptibility of debris flows in Yongji county.

The study area lies between 125°48′09″E to 126°40′01″E longitude and 43°18′07″N to 43°35′00″N latitude. There are four landforms in the entire area: middle mountains, low mountains, platform, and river valley. From southeast to northwest, the landforms of the study area are middle mountains, low mountains and platform with the altitudes ranging from 1,386 to 182 m. In addition to several andesites and metamorphic rocks, the main rock type is Yanshan Early Granite. The study area lies in the Tianshan–Xingan geosyncline fold area of the Jilin and Heilongjiang fold system (Qin et al., 2019). Folds and faults are relatively developed in Yongji county, which provides conditions for the occurrence of geological disasters (Figure 2A). Yongji county is in the mid-latitude subtemperate continental climate zone with an annual average precipitation of 722.75 mm. There are 39 rivers covering an area of more than 20 km². The main rivers include the Yinma River, Wende River, Chalu River and Aolong River.

A debris flow inventory map is a prerequisite for DFSM(Xu et al., 2012; Arabameri et al., 2020; Dash et al., 2022). A total of 129 debris flows were collected based on field surveys and historical materials. Figure 2A shows that debris flows are mainly distributed across the southeast mountain area. Statistics show that among 129 debris flows, only 7 are medium in size and 122 are small. In recent years, the increase in debris flow frequency in Yongji county has been closely related to deforestation and reclamation. With the destruction of forest vegetation, rainfall is more likely to cause soil erosion, which gradually forms a series of gullies. These gullies provide circulation conditions for debris flows. Figures 2B–E shows some images of occurred debris flows in the study area.

In this study, the extraction of watershed units was completed in ArcGIS 10.2 software (Tien Bui et al., 2015; Cao et al., 2019). The most commonly used watershed units (HWUs) are classified by the hydrological principles (Fausto Guzzetti et al., 1999). HWUs are derived based on an 8-direction flow algorithm (Horton et al., 2013). Establishing the HWUs consists of the following six steps: 1) filling the original DEM, 2) extracting the flow direction, 3) calculating the flow accumulation, 4) extracting river networks based on a threshold, 5) stream linking, 6) dividing HWUs based on flow direction and stream linking. The detailed classification process is shown in Figure 3A.

In addition, watershed units can be generated based on the mean curvature method (CWUs). The mean curvature is a simple combination of profile curvature and plan curvature. Its maximum and minimum values can indicate the changes in aspect and slope positions at the same time. Therefore, the mean curvature can reflect the ridge line, valley line, platform edge and wide valley edge (Romstad and Etzelmüller, 2012). Establishing the CWUs consists of the following five steps: 1) smoothing the original DEM, 2) calculating the mean curvature, 3) extracting the flow direction, 4) filling depressions based on flow direction data, and 5) dividing CWUs based on flow direction and depressions. The detailed classification process is shown in Figure 3B.

For HWUs, the number and size are closely related to DEM resolution and flow threshold, but for CWUs, the control factor is only DEM resolution. For HWUs, flow threshold values of 500, 1,000, 2000, 5,000, and 10,000 were chosen based on a DEM with a resolution of 30 m. For CWUs, we resampled the DEM with resolutions of 50, 100, 200, 300, 500, and 1,000. To ensure that the number and size of the two types of watershed units were not much different and consistent with the actual watersheds, a flow threshold of 1,000 and a DEM resolution of 300 were selected to divide the watershed units. For the HWUs, the study area was divided into 1,092 watershed units. The smallest unit was 0.10 km², the largest unit was 13.63 km², and the mean size was 2.40 km² (Figure 4A). For CWUs, the study area was divided into 1,211 watershed units. The smallest unit was 0.11 km², the largest unit was 8.87 km², and the mean size was 2.17 km² (Figure 4B).

The occurrence of debris flows is affected by many factors including topographic, geomorphologic, geological, ecological and meteorological factors (Zhang et al., 2012; Bregoli et al., 2014; Hu et al., 2014). Based on field observations, available literature and expert experience, fourteen influencing factors were considered, such as watershed area, relative height difference, watershed average elevation, watershed slope, mean curvature, fault density, river density, stream power index (SPI), topographic wetness index (TWI), plan normalized difference vegetation index (NDVI), landforms, precipitation, land use and lithology.

Because substantial collinearity will lead to model instability, collinearity analysis is essential before influencing factors are applied for DFSM(Qiu et al., 2022). Person’s correlation coefficient was calculated to test the collinear relationship among these factors, and the results are shown in Table 1. There is no correlation coefficient when the absolute value is less than 0.7 (Dormann et al., 2013; Yao et al., 2022). There was high collinearity between relative height difference and watershed average elevation, and the Person’s correlation coefficient was 0.86. In addition, the value of collinearity between watershed average elevation and precipitation was 0.69. Therefore, the watershed average elevation was eliminated.

The watershed area, relative height difference, watershed slope, mean curvature, SPI and TWI were extracted from the DEM with a resolution of 30 m. Fault, river, and lithology data were acquired from the geological map of Yongji county and field investigations. The Landsat 8 image taken on 11 August 2021, was used to produce the NDVI. Landforms, precipitation, and land use were provided by government reports. Thirteen influencing factors were converted to a grid cell with a resolution of 30 m in ArcGIS 10.2 (Chen et al., 2017). Table 2 shows date source and scale of influencing factors. When watershed units are applied to DFSM, grid patterns for each factor need to be transferred to the corresponding watershed units. For watershed area, geometric calculation in the attribute table was used to calculate the area of each watershed. The difference between the highest and the lowest points in each watershed was calculated as a relative height difference (Qin et al., 2019). For watershed slope, mean curvature, SPI, TWI, and NDVI, the zonal statistics tool in the spatial analysis was used and the statistical type was “mean.” The length of faults and rivers in each watershed was extracted by using the intersection tool, and then, the fault density and river density in each watershed were calculated using the field calculator. Precipitation for each watershed was determined based on the principle of majority, and this principle was also applied to factors of landforms, land use and lithology. The data types of precipitation, landforms, land use and lithology are discrete, while the data types of other factors are continuous. The influencing factor layers based on HWUs with a flow threshold of 1,000 are shown in Figure 6.

The FRs of the influencing factor subclass were used as the input variable of the DFSM models (Huang et al., 2020). Based on a series of previous studies (Xu et al., 2012; Aditian et al., 2018; Vakhshoori et al., 2019; Chang et al., 2020), we divided the continuous factor into eight levels using the natural fracture method. Taking HWUs with a flow threshold of 1,000 as an example, the FRs for each level of thirteen factors are shown in Table 3.

Logistic regression (LR) may be the most widely used statistical technique in susceptibility assessment (Colkesen et al., 2016). As a multivariate regression method, LR can find a model to describe the relationship between multiple independent variables and a dependent variable (Lee and Pradhan 2006; Lee 2007; Pourghasemi et al., 2013). For DFSM, the influencing factors are considered the independent variables and the occurrence and non-occurrence of debris flows are considered the dependent variables. For LR, variables may be continuous, discrete or arbitrary combinations of two types (Lee, 2007). LR can be expressed as follows (Ayalew and Yamagishi 2005; Yalcin et al., 2011; Schlögel et al., 2018): P = 1 1 + e − z (1) Z = α + β 1 x 1 + β 2 x 2 + … β n x n (2) where P denotes the probability of a debris flow occurrence in each watershed, varying between 0 and 1; Z represents the dependent variable including non-debris flows (0) and debris flows (1); α represents the intercept of the regression function, β 1 , β 2 , … β n are the regression coefficients; and x 1 , x 2 , … x n are the debris flow influencing factors.

Multilayer perceptron (MLP) is a kind of artificial neural network and has been widely used in classification (Tien Bui et al., 2015; Pham et al., 2017). The MLP generally consists of three main components, namely, input layers, hidden layers, and output layers (Kavzoglu and Mather 2003). For DFSM, the input layers are considered the influencing factors of debris flow, the output layers are considered the classification result of inferring debris flow or non-debris flow, and the hidden layers are considered the classification layers that convert input into output. The MLP model with only one hidden layer is the most basic three-tier structure model, which can fit and predict many non-linear problems (Li et al., 2019). In this study, a single-hidden-layer MLP model is used in DFSM. For example, n 0 , n 1 and n 2 represent the number of input, hidden and output layers, respectively, and the input variables are X = [ x 1 , x 2 … , x n 0 . Then, the input and output of the hidden layer are: (Li et al., 2019; Huang et al., 2020): z j = ∑ i = 1 n 0 w i j x i + b j (3) y j = f z j = 1 + e − z j − 1 (4) where z j , b j and y j represent the y th input variable, threshold value and output variable of the neuron in the hidden layer, respectively, w i j represents the weight value between the i th input neuron and the j th neuron in the hidden layer, and f z j represents the activation function. Then the input and output of neurons in the output layer are: z k = ∑ j = 1 n 1 w j k y i + b k (5) y k = z k (6) where z k , b k and y k represent the j th input variable, threshold value and output variable of the neuron in the output layer, respectively, w i k represents the weight value between the j th neuron in the hidden layer and the k th output neuron.

The decision tree model is a technique that uses a tree structure to discover and describe structural patterns in data. It does not require a preestablished relationship between all input variables and a target variable (Hitoshi Saito and Matsuyama, 2009). As an algorithm of the decision tree model, classification and regression tree (CART) was first proposed by Breiman et al. (1984) The CART consists of a root node, a set of internal nodes and a set of leaf nodes. The leaf nodes correspond to the classification result, and the other nodes correspond to the classification rules. CART was selected as the decision tree model in this study in view of its performance efficiency (Wang et al., 2015).

The Bayesian network (BN) is a graphical model for probabilistic relationships among a set of variables (Song et al., 2012). BN can be represented by directed acyclic graphs and conditional probabilities, reflecting the independent and interdependent relationship among various variables. The calculation formula is given as follows (Han et al., 2019): P L , M , N = P L × P M | L × P N | L , M (7) where P L is the prior probability, indicating the conditional probability without the parent node, P M | L is the conditional probability, indicating the occurrence probability of M under L conditions and P N | L , M is the conditional probability, indicating the occurrence probability of N under L and M conditions.

In this study, three commonly used criteria, including the predictive accuracy (ACC), the Kappa coefficient and the area under the receiver operating characteristic curve (AUC) were used to evaluate the prediction ability of DFSMs. The calculation of the three criteria is based on the confusion matrix (Ghasemian et al., 2022a). The confusion matrix, also known as the error matrix, is a standard format for accuracy evaluation. The confusion matrix can represent the difference between the model prediction results and the actual observation results (Xiong et al., 2020). In this study, the confusion matrix of the debris flow susceptibility predictive models is shown in Table 4. For example, a true positive (TP) suggests that the prediction result is ‘Debris-Flow’, and the actual observation result is ‘Debris-Flow’.

The predictive accuracy (ACC) represents the ratio of correctly predicted observations to total observations. This index shows how well the debris flow model works:

The Kappa index is used to assess the acceptability of debris flow models which can be calculated by: K = p 0 − p e 1 − p e (9) p e = A × a + B × b n 2 (10) where K is the Kappa coefficient; p 0 is overall classification accuracy, namely, ACC; A is the actual number of debris flows, and it is also the sum of TP and FN (Table 4); B is the actual number of non-debris flows, and it is also the sum of FP and TN; a is the predicted number of debris flows, and it is also the sum of TP and FP; b is the predicted number of non-debris flows, and it is also the sum of FN and TN; n is the total number of samples, and it is also the sum of A and B or a and b . The K varies from 0.0 to 1.0, the higher the K value is, the better the classification accuracy of the model.

The receiver operating characteristic (ROC) curve and area under the curve (AUC) can compare the prediction performance of different classifiers (Akgun et al., 2012). The abscissa and ordinate of the ROC are the false-positive rate (FPR) and true-positive rate (TPR) respectively. They can be obtained from the following equations (Pourghasemi et al., 2013): F P R = F P F P + T N (11) T P R = T P T P + F N (12)

AUC represents the quality of models that reliably predict the occurrence or non-occurrence of debris flows. The AUC varies from 0.5 to 1.0, and the higher the AUC value is, the better the prediction performance of the model.

The whole analysis process was implemented in IBM SPSS software (Sun et al., 2019; Sun et al., 2021). For LR, the forward step mode was adopted to screened variables. For BN model, the mechanism type was Tree Augmented naive Bayes (TAN), and Bayesian adjustment of small cell count was selected as a parameter learning method. For MLP, one hidden layer was selected, and the maximum training time was used as the termination rule. For CART, the maximum tree depth was set to 10, and percentage was used as the termination rule. Other parameters are default.

The Chi-Squared statistic was employed to identify the most important factors affecting the occurrence of debris flows in the study area (Ghasemian et al., 2022a). Figure 7 shows that landforms have the highest impact (128.5) on debris flows in the study area, followed by relative height difference and SPI (114.5), watershed slope (110.5), TWI (82.5), precipitation (74.5), lithology (60), land use (43.5), watershed area (27.5), river density (26.5), mean curvature (24), NDVI (23), and fault density (21.5).

According to field surveys and historical materials, a total of 129 debris flows were collected. Meanwhile, 129 non-debris flows were selected, which were at least 500 m away from the nearest debris flow (Figure 4) (Dou et al., 2019; Sun et al., 2020). Assigned 1 and 0 for debris flows and non-debris flows, respectively. The FRs of the thirteen influencing factors shown in Table 2 were taken as the input variables, and the debris flows and non-debris flows were taken as the output variables. For all 258 samples, 70% (n = 180) were selected randomly for training data, which were used to create the DFSM models. The remaining 30% (n = 78) were used as testing data, which were applied to validate the DFSM models. Based on two types of watershed units (HWUs and CWUs) and four models (LR, MLP, CART and BN), eight DFSMs of Yongji county were completed.

In this paper, IBM SPSS software was chosen to build the debris flow susceptibility predictive models. The model outputs are the debris flow susceptibility indices of all watershed units in the study area. Debris flow susceptibility indices are the probability of debris flow occurrence which varies from 0 to 1 (Xiong et al., 2020). Based on the ArcGIS software, the debris flow susceptibility indices were converted into raster format to produce the debris flow susceptibility map. Quantile classification was applied to divide the final maps into five classes, namely, very low susceptibility (VL), low susceptibility (L), moderate susceptibility (M), high susceptibility (H), and very high susceptibility (VH). (Martha et al., 2013; Hussin et al., 2016; Steger et al., 2017).

As shown in Figure 8, the susceptibility distributions of the eight models have common characteristics. Very high and high susceptibility areas are mainly distributed in the southeast, moderate susceptibility areas are mainly distributed in the middle, and very low and low susceptibility areas are mainly distributed in northwestern of Yongji county, which is consistent with previous research results (Qin et al., 2019). The landform in the southeast of the study area is mainly middle mountains, and the land use is mainly forest and farmland. The watershed units distributed in the southeast have large relative height differences and slopes, which leads to frequent debris flow disasters. The precipitation decreases from southeast to northwest, which is consistent with the susceptibility distribution. The lithology in southeastern Yongji county is hard massive rock, mainly granite. Weathered granite is a component of debris flows, which increases the density and destructive power of debris flows (Figures 2B–E).

For the eight DFSMs, the area ratios of the five susceptibility classes (very high, high, moderate, low, and very low) were 12.85–19.96, 13.82–21.51, 14.87–23.79, 17.70–28.33, and 16.56%–36.70%, respectively (Figure 9A); The debris flow ratios of the five susceptibility classes were 62.79–76.74, 12.40–24.81, 6.98–14.73,0.78%–3.10% and 0%–0.78%, respectively (Figure 9B). As shown in Figure 9C, the debris flow density was calculated to evaluate the performance of the DFSMs, that is, the ratio of debris flow percentage to area percentage on each susceptible class (Pham et al., 2016). The maximum values of the debris flow density of the eight models appear in the very high susceptibility class, varying from 3.15 to 5.61. The minimum values all appear in the very low susceptibility class, varying from 0.00 to 0.04. The debris flow density increases gradually from a very low class to a very high class, which provides a good visualization of the spatial predictions of debris flows (Pham et al., 2017; Asadi et al., 2022).

Model validation is a vital step in disaster susceptibility mapping (Wang et al., 2022). By considering the three commonly used performance metrics of ACC, AUC and Kappa, eight models were verified. The AUC, ACC and Kappa coefficient values of the eight models on the training and testing data are shown in Figure 10.

In the training phase, when HWUs were used as the mapping unit, the ACC stated that HWUs_CART model had the highest value (0.990), followed by HWUs_MLP (0.932), HWUs_BN (0.919) and HWUs_LR (0.881). It showed that the HWUs_CART model can correctly classify the debris flow and non-debris flow locations as debris flow and non-debris flow situations respectively. The highest and lowest Kappa values were 0.980 and 0.762, respectively for the HWUs_CART and HWUs_LR. Meanwhile, HWUs_MLP (0.864) and HWUs_BN (0.838) was ranked in other positions. In terms of AUC, results indicated that the HWUs_CART model with a value of 0.991 had higher performance than the HWUs_MLP (0.982), HWUs_BN (0.970) and HWUs_LR (0.966). When CWUs was used as the mapping unit, the ACC, Kappa and AUC values of the CWUs_CART model were 0.980, 0.960 and 0.985, which showed that the performance of the CWUs_CART model was the highest, followed by the CWUs_MLP (0.893, 0.787, 0.967), CWUs_BN (0.899, 0.778, and 0.946) and CWUs_LR (0.858, 0.716, and 0.946) (Figure 10). Although the results showed the excellent performance for all the four algorithms, the CART had the highest ability in debris flow classification and susceptibility mapping in the study area. In terms of watershed unit, ACC, Kappa and AUC values decreased when HWUs was replaced by CWUs, indicating that HWUs were more suitable for DFSM in the study area than CWUs.

Right side of Figure 10 showed the prediction capabilities of the eight models based on testing dataset. These results are very important for evaluating the applicability and robustness of the models. When HWUs were used as the mapping unit, the highest value of ACC was 0.834 for the HWUs_CART model, next for the HWUs_MLP (0.826), HWUs_BN (0.808) and HWUs_LR (0.802) models. The Kappa for the HWUs_CART model was 0.668 as the highest value, whereas this value was 0.652, 0.617, and 0.604 for HWUs_MLP, HWUs_BN, and HWUs_LR, respectively. The highest and lowest AUC values were 0.942 and 0.882, respectively for the HWUs_CART and HWUs_LR. Meanwhile, HWUs_MLP (0.902) and HWUs_BN (0.892) was ranked in other positions. Correspondingly, ACC, Kappa and AUC from CWUs were shown in Figure 9, which indicated a similar result with HWUs. CART model resulted in the highest ACC, Kappa and AUC values of 0.832, 0.663, and 0.909, which manifested it is the best model for the study area. At the same time, the HWU-based models had better performance than the CWU models for DFSM in the study area.

The results of the models are tested by one-way ANOVA in SPSS. For HWUs, there are significant differences between CART and each of the three methods (LR, MLP, and BN). There are no significant differences among LR, MLP, and BN. For CWUs, there are no significant differences between MLP and each of the two methods (LR and BN). There are significant differences between the other methods.

As shown in Figure 3, the extraction processes of HWUs are more complex than those of CWUs, because HWUs require six steps while CWUs require five steps. Model builder in ArcGIS is a workflow that connects a series of geoprocessing tools (Qin et al., 2019). It takes the output of one tool as the input of the other tool. Model builder can greatly reduce operation time and improve work efficiency. We had built two workflows for the processes of extracting HWUs and CWUs in the model builder. Experiments on two types of watershed units showed that HWUs extraction required 17 s, while CWUs extraction required only 3 s. In addition, for the division of HWUs, the influence of DEM resolution and flow threshold needs to be considered, while for CWUs, only DEM resolution needs to be considered. In summary, it takes more time and effort to extract HWUs than CWUs.

There are also significant differences between the two types of watershed units in the classification results. As shown in Figure 4, HWUs are mostly strip-shaped and widely different in size, while CWUs are nearly square and relatively uniform in size. The watershed unit boundaries extracted by the hydrological analysis method in areas with obvious topographic changes are relatively consistent with reality (Figure 11A). However, there are many small and narrow units in flat areas, because there is no clear flow direction in flat areas for hydrological analysis (Massimiliano et al., 2016) (Figure 11C). For the curvature method, the division of watershed units in flat areas is relatively satisfactory, and there is no parallel line problem similar to the hydrological analysis method (Figure 11D). However, in areas with obvious topographic changes, the boundaries of watershed units do not match well with the actual situation (Figure 11B).

Four models, LR, MLP, CART and BN, were used to complete the DFSMs of Yongji county in this study. Figure 10 shows the predictive ability of the eight models. When HWUs were used as mapping units, CART had the highest evaluation criteria with AUC, ACC and Kappa values of 0.991, 0.990, and 0.980 respectively, followed by MLP (0.982, 0.932, 0.864), BN (0.970, 0.919, 0.838) and LR (0.966, 0.881, 0.762) in the training stage. For the testing stage, the CART had the highest prediction accuracy with AUC, ACC and Kappa values of 0.942, 0.834 and 0.668 respectively, followed by MLP (0.902, 0.826, 0.652), BN (0.892, 0.808, 0.617) and LR (0.882, 0.802, 0.604). When CWUs were used as mapping units, the evaluation results showed the same trend as HWUs. The comparisons of the four evaluation models show that the CART had the best predictive ability over the other three models. The current research was in agreement with previous research results. Wang et al. (2015) analyzed landslide susceptibility based on five mathematical models (artificial neural network, frequency ratio, CART, LR and weights of evidence methods) and three sampling strategies. They indicated the results obtained from CART show steady prediction power with an AUC value larger than 0.7. Felicísimo et al. (2012) indicated that the CART is one of the most predictive models with the AUC value of 0.77. Using random forest (RF), boosted regression tree (BRT), classification and regression tree (CART), and general linear (GLM), Youssef et al. (2015) found the success rate for CART was 0.816 and for the prediction rate the CART was the highest with a value of 0.862. CART represents information in an intuitive and easy visual way, and is widely used in many fields (Bevilacqua et al., 2003; Malinowska 2014; Kim et al., 2015; Youssef et al., 2015; Yang et al., 2016).

Several studies have been conducted in similar areas. Selecting the frequency ratio (FR) model as the statistical method, Qin et al. (2019) explored the accuracy and practicability of HWUs and grid cell units (GCUs) in evaluating debris flow susceptibility in Yongji county. Qiao et al. (2021) built debris flow susceptibility models via the analytical hierarchy process (AHP) method and generated maps of Yongji county. The AUC values of the testing data in different studies are shown in Table 5. When the HWUs were chosen as mapping units, the AUC values of the DFSMs based on AHP, FR, LR, BN, MLP and CART were 0.812, 0.879, 0.882, 0.892, 0.902, and 0.942 respectively. The main difference among these DFSMs is the selection of different evaluation models, which indicates that machine learning models can improve the prediction accuracy of DFSMs. These results are consistent with previous studies, indicating that machine learning models are more suitable for DFSM than heuristic and general statistical models (Huang et al., 2020; Sun et al., 2021).

The selection of mapping units is one of the key issues for the rationality and correctness of disaster susceptibility mapping (Van Den Eeckhaut et al., 2009; Chen et al., 2019; Sun et al., 2020). The impact of different mapping units on disaster susceptibility mapping is greater than that of statistical methods (Zezere et al., 2017). Although many studies have tried to compare different evaluation models for disaster susceptibility mapping (Achour et al., 2018; Liang et al., 2020; Xiong et al., 2020; Dash et al., 2022; Qiu et al., 2022), very few studies have considered different mapping units. Qin et al. (Qin et al., 2019) explored the effect of grid cell unit and HWUs on the susceptibility mapping of debris flow, they found HWUs can reflect the geological and geomorphic environmental conditions of a debris flow accurately and perfectly. Li et al. (Li et al., 2021) discussed the influence of four different HWUs on debris flow susceptibility assessment results. The results show that the appropriate watershed division scheme can obtain more reasonable results. In this study, HWUs and CWUs were selected to map debris flow susceptibility. When the CART was selected as the machine learning model, the HWUs generated high AUC, ACC, and Kappa for training data (0.991, 0.990 and 0.980) compared to the CWUs (0.985, 0.980, and 0.960). For testing data, the AUC, ACC, and Kappa of HWUs were 0.942, 0.834, and 0.668, respectively. The AUC, ACC, and Kappa of CWUs were 0.909, 0.832, and 0.663, respectively. The results suggest that the HWU model has a higher debris flow prediction performance than the CWU model. The same trend can be observed in the LR, MLP, and BN models. Therefore, the HWU-based model is superior to the CWU-based model in debris flow susceptibility assessment due to higher training and testing accuracy.

As described in “6.1 Watershed unit classification processes and results comparison,” compared with CWUs, HWUs agree well with the actual watershed units in mountainous areas, but small and narrow units appear in plain areas. Since the frequency of debris flows in mountainous areas is much higher than that in plain areas, the division of watershed units in mountainous areas is more important than that in plain areas. Therefore, the HWU model is more practical than the CWU model. CWUs can also represent the distribution of watersheds and can be used as an alternative scheme.

Although this paper discussed the application of two types of watershed units in DFSM and obtained positive results, there are some limitations: 1) the number of debris flows is small, and 2) only HWUs with a threshold of 1,000 and CWUs with a resolution of 300 are selected for comparison. In future research, we will constantly update the debris flow database to improve the data quality. Moreover, it is necessary to explore the similarities and differences of multiscale watershed units in DFSM.