Debris flow susceptibility mapping using the Rock Engineering System (RES) method: a case study

The main purpose of the present study is to develop a debris flow susceptibility map of a mountain area (Susa Valley, Western Italian Alps) by using an upgraded version of the Bonetto et al. (Journal of Mountain Science 18, 2021) approach based on the Rock Engineering System (RES) method. In particular, the area under investigation was discretized in a 5 × 5-m grid on which GIS-based analyses were performed. Starting from available databases, several geological, geo-structural, morphological and hydrographical predisposing parameters were identified and codified into two interaction matrices (one for outcropping lithologies and one for Quaternary deposits), to evaluate their mutual interactions and their weight in the susceptibility estimation. The result for each grid point is the debris flow propensity index (DfPI), an index that estimates the susceptibility of the cell to be a potential debris flow source. The debris flow susceptibility map obtained was compared with those obtained from two expedited and universally recognized susceptibility methods, i.e. the Regional Qualitative Heuristic Susceptibility Mapping (RQHSM) and the Likelihood Ratio (LR). Each map was validated by using the Prediction Rate Curve method. The limitations and strong points of the approaches analysed are discussed, with a focus on the innovativeness and uniqueness of the RES. In fact, in the study site, the RES method was the most efficient for the detection of potential source areas. These results prove its robustness, cost-effectiveness and speed of application in the identification and mapping of sectors capable of triggering debris flow.


Introduction
Debris flows are rapid landslides of mixed and unconsolidated sediments and water which occur when soil and rock fragments become saturated and flow down into a steep channel driven by the force of gravity. Debris flow events can be extremely dangerous for humans and infrastructures due to their high velocity (up to 20 m/s), their large mobilized volumes (even more than 10 9 m 3 ) and their unpredictability (Varnes 1978;Hutchinson 1988;IAEG 1990;Cruden and Varnes 1996). A careful territorial analysis should include the identification of potential debris flow source areas for correct land use and risk management, especially in mountain regions where these aspects become fundamental for the resilience of the rural areas and to tackle the effects of climate change.
Geological, geomorphological, hydrogeological and landslide maps (Soeters and Van Westen 1996;Fell et al. 2008;Corominas et al. 2014), supported by direct field observations, are fundamental to detect and delimit zones susceptible to landslide triggering (including debris flow). Several models have been proposed in the scientific literature for evaluating landslide susceptibility by combining geo-environmental factors and landslide spatial distribution (Brabb 1987;Soeters and van Westen 1996;Carrara et al. 1984Carrara et al. , 1999Carrara et al. , 2008Aleotti and Chowdhurry 1999;Guzzetti et al. 1999;Dai and Lee 2001;Chacón et al. 2006;Fell et al. 2008;Reichenbach et al. 2018). In general, these models are based on the identification of those factors that contribute to landslide triggering, by distinguishing between predisposing and triggering factors (e.g. Costa and Jarrett 1981;Hutchinson 1992;Cruden and Varnes 1996;Jakob and Hungr 2005;Hungr 2005;Van Westen et al. 2008;Pereira et al. 2012;Corominas et al. 2014;Iverson 2014Iverson , 1997Hungr et al. 2014 and references herein). Data related to triggering factors represent an important set of input parameters for landslide hazard assessment, while the predisposing factors and landslide inventories play a key role in landslide susceptibility analysis (Dai and Lee 2001;Clerici et al. 2002;Corominas et al. 2003Corominas et al. , 2014Van Westen et al. 2006).
Susceptibility analysis can be assessed through both qualitative (inventory-based and knowledge-driven methods) and quantitative (data-driven methods and physically based models) methods (Carrara et al. 1995;Soeters and Van Westen 1996;Guzzetti et al. 1999Guzzetti et al. , 2006aDai and Lee 2001;Van Westen et al. 2006;Clerici et al. 2002). Inventory-based methods provide a multitemporal landslide distribution (spatial and temporal frequencies) based on historical series and represent a key starting point for hazard mapping and risk assessment. The analysis of past debris flow events provides useful information for forecasting future debris flow, based on topographic, geological and geomorphologic characteristics. The knowledge-driven approaches, or heuristic methods, are based on the expert knowledge of landslide mechanisms that allows the degree of instability, combining geomorphological observations and thematic geological maps to be determined (Abella and Van Westen 2008;Nachbaur and Rohmer 2011). This approach can be applied when the landslide inventory is incomplete or when a preliminary and expedited evaluation is needed.
In data-driven landslide susceptibility methods, statistical and quantitative predictions are made using the records of past landslide events through three different approaches: bivariate, multivariate and active learning statistical analyses. Landslide distribution could define the functional relationships among known or inferred instability factors (Malusà and Mosca 2002;Guzzetti et al. 1999;Huabin et al. 2005;Chacón et al. 2006;van Westen et al. 2008), but results are strongly influenced by the quality and completeness of the inventory (Guzzetti et al. 2006a, b;Reichenbach et al. 2018).
The physically based models are based on the mathematical modelling of landslide failure and a set of numerical parameters that describe the geometry, and the internal and external slope condition. Slope analysis can be done using simple limit equilibrium models, such as the infinite slope model, or more complex approaches like kinematics analysis or numerical modelling (Montgomery and Dietrich 1994;Rigon et al. 2006;van Asch et al. 2007;Simoni et al. 2008;Baum and Godt 2010;Anagnostopoulos and Burlando 2012;Alvioli and Baum 2016).
For susceptibility mapping purposes, the choice of the methods should be driven by the availability of the input parameters required for the analyses and their mutual interactions. Recently, Bonetto et al. (2021) proposed an application of the Rock Engineering System (RES; Hudson 1992) to debris flows analysis. The RES approach has also been used for other landslide types (Mazzoccola and Hudson 1996;Kim et al. 2008;Rozos et al. 2008;Tavoularis et al. 2017Tavoularis et al. , 2021Xiao et al. 2020;Pokharel et al. 2021 and references herein) and implemented in neural network approaches (Wang et al. 2014;Meten et al. 2015). Bonetto et al. (2021) combined inventory, expert evaluation and data-driven methods for a quantitative evaluation of the debris flow propensity at a basin-scale RES methodology were then used to evaluate the Debris flow Propensity Index (DfPI) by quantifying and scoring the mutual interaction between predisposing parameters. Starting from these assumptions and promising results, in this paper, the authors propose an upgraded version of the Bonetto et al. (2021) methodology applied to debris flow by considering new parameters for the DfPI determination and by implementing the procedure in GIS environment. The predisposing parameters are encoded into two interaction matrixces to consider outcropping lithologies and Quaternary deposits and the DfPI values are mapped onto a 5 × 5-m grid cell resolution. The procedure has been tested on the same sector of the western Italian Alps (Upper Susa Valley) as in Bonetto et al. (2021) for a direct comparison and validation of the results.
The susceptibility map obtained is compared with those obtained by two methods available in the scientific literature: the Regional Qualitative Heuristic Susceptibility Mapping (RQHSM) method proposed by Soeters and van Westen (1996) and the Likelihood Ratio (LR) method (Lee 2004;Regmi et al. 2010;Sujatha et al. 2013;Kanungo et al. 2011;Akgun 2012). The limitations and strong points of the methodologies are discussed by comparing susceptible areas with an available debris flow source area database. Furthermore, the differences between the three methods are quantitatively assessed by using the "Prediction Rate Curve" method (Chung and Fabbri 2003).

Study area: the Upper Susa Valley
The study area is the Upper Susa Valley, in the Western Italian Alps (Fig. 1A). Many basins of this valley are affected by recurrent debris flow events, which occurred especially during late summer and fall seasons (Tiranti et al. 2008(Tiranti et al. , 2014(Tiranti et al. , 2016. These phenomena have caused many damages to the crucial infrastructures, which link the valley (and more in general Turin) to France, and to the urbanized areas, which mainly developed on debris fans. This valley is drained by the Dora Riparia River, a left-hand tributary of the Po River, and has been carved by glaciers at the end of Pleistocene in units belonging to the Penninic domain of the Western Alps.
The Western Alps result from the collision between European and Adriatic plates following subduction and closure of their interposed Piemonte-Liguria oceanic basin (e.g. Dal Piaz 2010a, b and therein references). The upper Susa Valley exposes a tectonic stack of continental margin (Ambin Massif Auct., Pre-Piedmont and Briançonnais units) and oceanic Piemonte-Liguria units ( Fig. 1; Polino et al. 2002;Piana et al. 2017). The Ambin massif Auct. comprises two pre-alpine complexes: the Clarea and the Ambin complexes, resting at lower and [Continental margin units] 8-Re Magi Unit; 9-Chaberton-Grand Hoche-Grand Argentier unit; 10-Valfredda Unit; 11-Gad Unit; 12-Vallonetto Unit; 13-Permo-mesozoic succession of the Ambin Complex; 14-Ambin Complex; 15-Clarea Complex; 16-Gypsum and tectonic carbonate breccias; a -alluvial deposits; b -major landslide upper structural levels respectively (Polino et al. 2002;Malusà et al. 2002;Malusà and Mosca 2002). The Clarea complex consists of micaschist and paragneiss embedding metabasite and orthogneiss of dioritic composition. The Ambin complex is formed by gneiss and quartz-micaschist derived from volcanoclastic protoliths, several types of micaschist and bodies of aplitic gneiss. This complex is overlaid by a Permo-Mesozoic cover consisting of conglomeratic quartzite with pink quartz clasts and levels of sericitic schist (upper Permian), massive quartzite (lower Triassic) and a succession of dolomitic marbles and calcareous schists with local intercalation of carbonate breccias (Triassic to Cretaceous). The Vallonetto unit is characterized by a Briançonnais-type succession including marble and dolomitic-marble of Triassic age followed by Jurassic to Cretaceous calcschist. The Chaberton-Grand Hoche-Grand Argentier unit, belonging to the Pre-Piedmont zone, is formed by a thick dolomitic succession of Norian age followed by Rhaetian-Hettangian calcareous schists and then unconformable overlaid by prevailing calcschists with phyllites and beds of breccias (Jurassic to Cretaceous). The Piemonte-Liguria Oceanic units are formed by thick sequences of Upper Jurassic (?)-Cretaceous calcschist containing levels of micaschist and phylladic schist and embedding bodies of ophiolites (serpentinite, metabasite). Masses of gypsum and of carbonate breccias (Carniole Auct.) occur along the main tectonic contacts.
During the Quaternary, the Susa valley was carved with typical U-shaped cross section by the action of the glacier and, in its upper part, glacial deposits are recognizable on the valley flanks. Following the retreat of glaciers, several deep-seated gravitational slope deformations (DSGSD) developed on the valley flanks.

Data availability
The data used in this paper can be classified into two main categories: (1) landslide inventory and (2) thematic GIS-based maps. Data were collected using freely accessible and available datasets in national and regional geodatabases. Since the RES approach applied to debris flow requires the identification of environmental predisposing factors, thematic GIS-based maps were drawn starting from the available spatial datasets. Compared to the original version proposed by the authors (Bonetto et al. 2021), in this study, ten predisposing factors were considered for the definition of the DfPI (Table 1): bedrock lithology, quaternary deposits, lineament density, slope angle, curvature, elevation, slope aspect, channel network, landslide activity and land use. For the sake of simplicity, these factors are grouped into geological, geomorphological and hydrogeological parameters and land use.

Landslide inventory and source area archive
Debris flow source areas were available in the RiskNat-Arpa Piemonte inventory that contains events occurred (among other) in the Upper Susa Valley from 1990 to 2015. The original source area map was updated by manually adding recent events, easily identifiable from the analysis of the past aerial photographs between 2017 and 2018. An amount of 846 points of debris flow source area were recorded in 390 km 2 in the study area (Fig. 2).

Predisposing factors
Geological parameters: bedrock lithology, quaternary deposits, lineament density The geology of Upper Susa Valley was derived by the Foglio 132-152-153 Bardonecchia at 1:50,000 scale (Polino et al. 2002) of the official Italian Geological Cartography project (Fig. 3). The shapefiles of bedrock lithologies, quaternary deposits and landslides were converted into raster formats with 5 × 5-m cell resolution. Based on the differences in strength and permeability of rocks and deposits and the presence of landslides, five classes for bedrock lithology and five classes for quaternary deposits were identified (Bonetto et al. 2021) (Table 2).
Geological discontinuities (i.e. fractures, faults and foliations) induce a decrease in the mechanical strength of rock mass and increase the capacity to produce loose debris (Ferrero et al. 2016;Caselle et al. 2020;Umili et al. 2020). A regional map of the rock fracturing degree with a traditional survey is not feasible, and consequently, the general condition of the rock masses can be described by remote extraction of the main lineaments (Tripathi et al. 2000;Jordan et al. 2005;Vaz et al. 2012;Bonetto et al. 2015;Umili et al. 2018). Lineament extraction can be realized through automatic or manual approach from digital elevation model or orthophotos. In the automatic approach, the linear traces are detected on a DTM by algorithms based on principal curvature values (Bonetto et al. 2015(Bonetto et al. , 2017 while in manual approach the user manually detects and tracks the lineaments. In this paper, manual lineament extraction was performed from visual interpretation of the orthophotos to increase lineament extraction accuracy, only focusing on bedrock outcrop areas. The track density map was derived ( Fig. 3B) by using the Line Density tool on QGis. This tool allows the measurement, within a given circular area, of the line density for each raster cell. This measure is obtained as the sum  Geomorphological parameters: slope angle, curvature, elevation, slope aspect, landslide activity The geomorphological layers were obtained from the DEM analysis using specific tools freely available in QGis software. The resulting maps are shown in Fig. 4.
Slope angle is the result of the combined influence of many agents (Huma and Radulescu 1978;Carrara 1983;Maharaja 1993;Rozos et al. 2008) and is one of the predisposing factors capable of triggering debris flows. Slopes ranging between 20 and 45° (Takahashi 1981;Hungr et al. 1984;Rickenmann and Zimmermann 1993) are characteristic values for source area location. Five different slope classes were selected (Fig. 4A) for the classification of the study areas: The terrain curvature is the curvature of a line formed by intersecting a plane with the terrain surface. Operatively, the curvature value is the reciprocal of the radius of curvature of the line. Debris flows generally start where curvature is concave (Wieczorek et al. 1997;Delmonaco et al. 2003) and the flow can be channelled into gullies. Consequently, a distinction between concave (negative values), convex (positive values) and flat surface (values near zero) was made (Fig. 4B). Elevation (Fig. 4C) and slope aspect (exposition) (Fig. 4D) do not directly influence the slope failure, but they are the result of tectonic activity and erosion process related to climatic condition (Rozos et al. 2008). The elevation was distinguished into classes with an elevation step of 500 m. In alpine regions, source areas are usually located at high elevation where deposits are concentrated.
The slope aspect reflects the exposition that is responsible for different local microclimatic conditions and solar exposition during the day. The exposition was classified with a step of 60° starting from the north (0 value).
The landslide activity map (Fig. 4E) shows the landslide distribution (not only debris flow but also every landslide type) and their state of activity: (i) active, (ii) quiescent, (iii) inactive and (iv) not defined. In these areas, the overall rock mass and deposit conditions might increase the loose material availability that could be mobilized during a flow event.
Hydrological parameter: distance from channel network Debris flow events require a channel to flow downstream. The flow erodes the channel leaves and banks, which are important sources of material available during the event. Usually, the areas prone to supplying material are located at distances lower than 200 m, while beyond these values the probability decrease (Rozos et al. 2008). For this reason, five different buffer zones were created to identify the distance along the river where material can be mobilized during the flow: (i) 0-50 m, (ii) 50-100 m, (iii) 100-150 m, (iv) 150-200 m, (v) > 200 m (Fig. 5).
Land use condition: land use The land use describes the vegetational, mechanical and hydrological characteristics that control the slope stability (Glade 2003;Reichenbach et al. 2004). The land use influences the soil behaviour during rainfall and the magnitude of potential mobilizable material. The dataset, in raster format, provided by the regional archive allows five classes to be identified: (i) grassland, (ii) lakes, (iii) high forest, (iv) low forest, (v) rock and deposits (Fig. 6). This classification reflects the different soil conditions related to erodibility or resilience to the impact of rainfall.

The RES methodology
The RES was proposed by Hudson in 1992. The methodology is based on a matrix approach that allows the numerical encoding of mutual interactions between the predisposing factors arranged along the diagonal terms (Pi) of the matrix (Fig. 7). Off-diagonal terms (I ij ) are scored with values from 0 (no interaction) to 4 (critical interaction) using the Expert Semi-Quantitative method (ESQ) (Harrison and Hudson 2006;Vagnon et al. 2015). The contribution of each parameter to the debris flow triggering is described by the weighting coefficient a i : where C is the sum of the values in each row (Cause -C) and it represents the influence of the parameter P i on the system, E is the sum of the values in each column (Effect -E) and it represents the influence of the system on parameter P i .
The debris flow susceptibility index (DfPI) is given by: where a i is the weighting coefficient calculated for each parameter using Eq. 1, and P iK corresponds to a specific value between 0 and 4 attributed to each class of the identified predisposing factors. The 0 value represents the most stable conditions (lower debris flow susceptibility) while the 4 value represents the most favourable conditions for debris flow triggering. The P iK value describes the weight assigned to each predisposing factor based on its propensity to induce instability (Mazzoccola and Hudson 1996; Rozos et al. 2008Rozos et al. , 2011. (1) In this paper, an upgraded version of the Bonetto et al. (2021) approach is proposed to develop a GIS-based debris flow susceptibility map for the whole Upper Susa Valley by using a 5 × 5-m resolution grid. Curvature, land use and landslide activities were considered as predisposing factors in addition to those of the original version: lithology, fracture network, quaternary deposits, slope and channel network. Two interaction matrices were created to separately analyse the mutual interaction between the bedrock lithology (matrix A) or deposits (matrix B) and the other parameters.
The GIS-based DfPI has the same range of values of the original DfPI version (from 0 to 100). Five susceptibility classes were defined by using a modified version of Brabb's susceptibility scale: low (0-20),

RES results
The off-diagonal terms of each matrix were coded considering the mutual influence between the parameters and a i were calculated respectively for matrix A (bedrock lithologies) and matrix B (quaternary deposits) (Tables 3 and 4). By applying Eq. 2 to each grid cell and analysing separately bedrock outcrop areas and quaternary deposits (Table 5), the map of the whole study area was obtained (Fig. 8). The RES method highlights that debris flow susceptibility zones are concentrated at high elevation and where talus deposits are present. In addition, a large concentration of high susceptibility areas near the channel network and along the slope is emphasized. The NW sector has the areas with the highest susceptibility values compared to the rest of the valley. In fact, these sectors are characterized by the presence of deformed dolomitic rocks affected by intense post-metamorphic brittle deformation, as highlighted by the presence of the high amount of talus deposits.
The basins analysed in Bonetto et al. (2021) and their correspondent global DfPI are highlighted in Fig. 8. The results highlight that (i) qualitatively, there is a good agreement between the global DfPI and the grid DfPIs and (ii) the grid DfPI allows the direct identification of the areas capable of triggering debris flow events, providing great advantages especially for the planning of risk management activities.

Comparison with other susceptibility methods and quantitative validation of the RES
Qualitative and preliminary analysis of the proposed RES approach is promising for the identification of the areas of the Upper Susa Valley most susceptible to debris flow triggering. However, for evaluating the reliability of this procedure, highlighting its potentialities and limitations, RES was firstly compared with two well-established susceptibility methods available in the scientific literature: the Regional Qualitative Heuristic Susceptibility Mapping (RQHSM) and the Likelihood Ratio (LR). Then, for a quantitative comparison and for evaluating the effectiveness of the RES, the RQHSM and the LR, the prediction rate approach (Chung and Fabbri 2003) was used. Figure 9 shows the logical steps followed for the quantitative comparison between the susceptibility mapping methods analysed in this work.

Fig. 5 Channel network and buffer with distance from channel
The RQHSM -regional qualitative heuristic susceptibility mapping The RQHSM is a heuristic method that allows defining a susceptibility index (SI) by scoring the predisposing factors using geological and geomorphologic criteria (Soeters and van Westen 1996;Riopel et al. 2006;Blais-Stevens et al. 2010, 2012, 2014 and references herein). For each parameter, different classes were defined, according to the propensity of triggering debris flows (Table 6). Following the equation proposed by Blais-Stevens and Behnia (2016), SI was calculated as: with geology (G), slope angle (S1), slope aspect (S2), drainage system (D) and plan curvature (C). The resulting debris flow susceptibility map is shown in Fig. 10. For a direct comparison with DfPI values, SIs were multiplied by 100 and classified following the same criteria.
The RQHSM highlights the areas near the drainage system with a high probability of failure, considering that the onset of debris flows is usually triggered in steep streams. In this case, only the materials into the stream and along the riverbanks are detected as susceptible at triggering phenomena.

The likelihood ratio (LR)
The likelihood ratio (Lee 2004;Lee and Pradhan 2007;Demir et al. 2015 and references herein) is a statistical method that correlates environmental conditions with landslide areas and extends the landslides spatial occurrence in similar setting areas.
Based on the assumption that future landslides will occur under the same conditions as past landslides, the statistical method allows defining which factors, or combination of factors, play a fundamental role in the landslide initiation. With the LR, it is possible to evaluate the relationship between the dependent parameter (landslide occurrence) and the independent parameter (such as geological, geomorphological and hydrogeological features) and retrieve a ratio between the landslide occurrence probability and the non-occurrence probability calculated for each class factor. For debris flow analysis, these terms correspond to the ratio between: where the landslide occurrence ratio is the ratio between the number of landslides in i-th class and the total amount of landslide in study (4) FR = Landslide occurrence ratio Area domain  area, the area domain is the ratio between the area of i-th class and the total area in our case study. If the FR is greater than 1, it means that factor class has a high correlation with the event and vice versa. Operatively, for evaluating the FR, the landslide inventory described in the "Landslide inventory and source area archive" section, was randomly divided into two sets in the proportion 70% (593)-30% (253): the training and the validation sets (Fig. 11). The training set was used to build the statistical model, while the validation set was used to validate the results.
Using Eq. 4, the FR values were calculated for each layer class ( Table 7). The susceptibility index for each cell was obtained as the sum of all the FRs calculated for each selected parameter: and the resulting debris flow susceptibility map is shown in Fig. 12. SI values were scaled up to vary between 0 and 100 for a direct comparison with DfPIs and SIs from the RQHSM. The LR method highlights three areas with high susceptibility values (> 80) located in the southwestern and central parts of the study area. These areas correspond to the high-altitude zones at the peaks and the crests of the mountain chain and to large areas occupied by detrital materials and talus deposits. In these areas, the structural setting of the rock masses (made of prevailing marbledolostone and calcschist) and their mechanical proprieties are favourable for producing a huge amount of loose material. Along the slope and in presence of the forests, the susceptibility values are low while the presence of talus deposits is noted as the most critical areas. These low values are attributed to high slope areas and zones in proximity to the channel network focusing on the areas with coarse deposits and high fractured rocks.

Model validation
In the previous section, a visual and qualitative comparison was made between the RES and the other models. In this section, quantitative analyses were performed to assess the validity of the three susceptibility methods. However, the assessment of the model robustness is always a challenging task. Many authors have proposed different methodologies for the comparison between forecasted results and observed data. In this study, the susceptibility maps were validated by using the Prediction Rate Curve (Chung and Fabbri 2003). This approach is based on the direct comparison between the source area estimated from susceptibility maps and the real source area from debris flow inventory. The validation set consists of 253 landslide source areas (30% of total source areas inventory) previously excluded from the statistical model analysis. The source area was overlaid on the susceptibility maps (Fig. 13) and the validation phase was performed plotting the LSI value in x-axis and the cumulate percentage of landslides on y-axis (Fig. 14).
The use of the same validation dataset for all three methodologies allows the direct comparison between the curves obtained and the evaluation of their robustness. The validation step verifies that the maximum number of landslides was included in the highest susceptibility classes. Results show that the landslide areas predicted by the RES method in the high-extreme susceptibility range (50-100) are 96% while in the same range, the performance of the RQHSM and LR methods is 82% and 77%. This means that in the RES methodology most landslides were predicted in areas characterized by high-extreme susceptibility while, in contrast, in the other methods more landslides were predicted in areas of low-high susceptibility. This fact reflects a lower prediction capability.

Discussion
Using predisposing parameters selected from national and regional geodatabases, debris flow susceptibility maps were developed using the RES, RQHSM and LR methods. Analysis of the prediction curves (Fig. 14) shows the trend differences between the analysed methodologies: the RES appears to be more accurate in predicting susceptible areas based on debris flow database analysis. To compare and evaluate the models, the percentage of the areas classified as susceptible is reported in Fig. 15. In the RES and RQHSM methods, the differences in the percentage of susceptible areas are negligible (0.2|0.8 Low, 13.2|17.7 Medium, 47.0|43.8 High, 33.3|32.2 Very high, 6.3|5.5 Extreme). These two methodologies classify most of the territory with high and very high susceptibility while, on the other hand, the LR approach identifies a different proportion of areas with respect to the other two methods (25.3 Low,44.2 Medium,16.8 High,10.2 Very high,3.3 Extreme). In this case, the LR method identifies few areas of high susceptibility by concentrating high values in a small portion of the area. The difference is also appreciable comparing the maps in Fig. 13, which shows the overlay of the susceptibility map and the location of the source areas used in the validation process.
A common feature of all three methods is that the spatial distribution of the most susceptible areas corresponds to the zone of maximum elevation with the presence of debris along the slope and steep sections. The RQHSM and RES methods identify the channel network as playing a significant role in the triggering phenomena, while the LR method highlights highly fractured outcrop rocks as a key factor. The low susceptibility zone was detected in the alluvial planes and flat areas, but the same susceptibility class in the LR was also detected along the slope and near channels.
The applicability of different methods to debris flow susceptibility depends on several factors (e.g. amount and quality of the data) but established standards and codes of practice are not available for the choice of the most appropriate method for landslide susceptibility evaluation.
The LR method is based on the evaluation of the relationships between predisposing factors (thematic layers) and the distribution of debris flow source areas collected in the past years (landslide inventory). Thus, a detailed landslide archive is a necessary input to apply the methodology. In remote mountain areas, where it is difficult to collect information, some events may be missed, and this issue can lead to an underestimation of the potential triggering areas with a consequent decrease in forecasting quality. The application of remote sensing survey can redress this problem, but the precise time of occurrence and the definition of magnitude remain undetectable. This information is not required for landslide susceptibility, but it is fundamental for hazard or risk assessment. Moreover, the visual detection of source areas has some limitations due to the misclassification of landslide events and underestimation of the number of events recorded. The growth of vegetation and action of snow and glaciers could partially or totally delete the evidence of past events. In addition, the complex terrain morphologies may not allow the proper identification of the source area if the quality of aerial photography is not sufficiently detailed.
The evaluation of landslide susceptibility requires geological, morphological and statistical input data to process the analysis. Much effort is required to collect and validate the necessary input data, which are not always available from open-access geo-environmental databases. Statistical methods need the available spatial-temporal datasets of landslide events. The RES and RQHSM methodologies are essentially based on predisposing parameters derived from geological, morphological and hydrogeological thematic layers.
In the RES methodology proposed, starting from available databases, several geological, geo-structural, morphological and hydrographic parameters were considered to quantify their mutual interaction and to define a debris flow susceptibility map. In particular, the parameters for the different lithological classes of the bedrock and for their degree of fracturing were quantified in the regional-scale tectonic setting, and for the Quaternary deposits.
The RES method, as applied in this study, proved to be low-cost and time-saving and allowed sectors with different propensities to triggering debris flow to be identified. Fig. 11 The validation (circles) and the model (yellow triangles) sets in which the landslide inventory was randomly divided

Conclusions
Susceptibility analysis in mountain areas is the first step in the risk assessment procedure. It is useful for local authorities to define potential areas that could be affected by debris flow phenomena. Different methodologies for susceptibility analyses are available in the scientific literature, but which to choose depends on the scale of analysis required and the amount and quality of available data. In this work, an innovative GIS-based application of Rock Engineering System (RES) was used for debris flow susceptibility mapping. Starting from the definition of a global DfPI for a single debris flow basin, as suggested by Bonetto et al. 2021, in this paper, a debris flow susceptibility map of the Upper Susa Valley (W Alps, Italy) was developed with a 5 × 5-m grid resolution. Important updates to the Bonetto et al. (2021) approach were proposed: (i) new predisposing factors, such as land use, landslide activity and other geomorphological aspects, were considered; (ii) two interaction matrices were proposed for considering the mutual interaction of the bedrock lithology or deposits with the other parameters; (iii) GIS-based mapping with evaluation of DfPI for each cell of the grid.
The susceptibility map obtained was compared with the RQHSM and LR methods. Numerical quantification of the validity of the susceptibility forecasting was performed by using the Prediction Rate Curve model, comparing the susceptibility model results with an available debris flow source inventory. In general, there was a good agreement between the forecasted and observed source areas. However, the RES-based map appears to be the most efficient and robust in the detection of source areas, since it was able to predict 96% of the source areas that fall into the high-extreme susceptibility range (50-100).
The study demonstrates that the application of the RES method offers an opportunity for initial debris flow susceptibility screening at medium and large scale, using available and open-access data, and meets the needs of authorities for land use management and planning. Further studies will be devoted to a comparison with other more sophisticated and complex methodologies (e.g. multivariate statistical approach or artificial neural network methods) developed in recent years.

Funding
Open access funding provided by Università degli Studi di Torino within the CRUI-CARE Agreement.

Conflict of interest The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.