A quantitative approach for debris flow inception and propagation analysis in the lead up to risk management

This paper presents a quantitative analysis of debris flows in weathered gneiss through a methodology to identify key aspects in the lead up to risk management. The proposed methodology considers both the triggering and propagation stages of landslide using two physically based models “Transient Rainfall Infiltration and Grid-Based Regional Slope-Stability Analysis” (TRIGRS) and “Smoothed Particle Hydrodynamics” (SPH), respectively. The TRIGRS analyses provided information about the location and initial volume of potential landslides. The SPH model, adopting the initial triggering volumes as input data, allowed the back analysis of the propagation stage in terms of both main pathway and depositional area. Both models can be easily implemented over large areas for risk assessment and are able to provide interesting information for the design of risk mitigation structures. Clearly, the rigorous implementation of these models requires the use of geotechnical data obtained from in situ and laboratory tests. When these data are not available, literature data obtained for similar soils for genesis and stress history with those studied can be used. The applicability of the methodology has been tested on two debris flows which occurred in 2001 and 2005 in the province of Reggio Calabria causing extensive damage involving various lifelines. The model results have been validated with the real events, in terms of both triggering/inception areas and debris fans using two dimensionless indices (Itrig and Idep). The indices compared, respectively, the real triggering/inception area with the simulated one and the real depositional area with the numerical one. For analysed phenomena, the values of Itrig were higher than 90% while Idep assumes values higher than 70% which support the applicability of the proposed methodology.


Introduction
Rainfall-induced debris flows are very rapid to extremely rapid phenomena characterised by a well-defined triggering area, flow path and depositional area. Debris flows may be triggered by a slide; they occur on established paths, usually gullies and first-or second-order drainage channels, and are characterised by a strong entrainment of material and water from the flow path. As these phenomena cause significant economic and human losses, the assessment of spatial and temporal landslide forecasting on a large scale is the key for assessing landslide risk and designing risk mitigation measures. The zoning of landslide hazard makes it possible to provide landslide hazard maps (Borrelli et al. 2018;Cascini et al. 2015;Ciurleo et al. 2016Ciurleo et al. , 2017Ferlisi and De Chiara 2016;Mandaglio et al. 2015;Mandaglio et al., 2016 a, b) and suggests the areas where protection measures are necessary to preserve the elements at risk.
The approaches used to pursue this aim are generally based on the correlation of historical records of landslide occurrence with the identification of landslide triggering and predisposing factors. Amongst the possible approaches available in scientific literature, quantitative modelling can be pursued by physically based models that, in case of debris flow analysis, should focus on both the triggering and propagation stages. During the triggering stage, the soil generally follows a solid-like behaviour and is studied by continuous mechanics; during the propagation stage, the soil-water mixture flows in the channel and the soil particles undergo significant displacement.
Most of the methods proposed in literature are based on grids, either structured (finite differences) or unstructured (finite elements and volumes). Generally, the methods used to solve small deformation problems (typically of pre-failure and failure stages) are based on the Lagrangian approach characterised by a mesh attached to the material and deforming with it. The Lagrangian approach is able to model history-dependent material behaviour but it suffers from mesh distortion in large-deformation problems (Soga et al. 2016). In order to avoid numerical problems associated with mesh distortion, large deformation problems (typically of the post-failure stage) are often based on methods using the Eulerian formulation of the problem. The Eulerian description is characterised by a computational mesh that is kept fixed in space while mass moves through it, but it is difficult to use with historydependent constitutive models (Soga et al. 2016) The coupled Eulerian-Lagrangian method is an arbitrary Lagrangian-Eulerian method used with the aim to capture the strengths of the Lagrangian and Eulerian methods in modelling large-deformation problems (Qiu et al., 2011;Soga et al. 2016) but its applicability, especially over a large area, is often limited by greater computational time.
Some of the interesting and promising meshless numerical methods available for modelling large deformation problems include material point method (MPM), smooth particle hydrodynamics (SPH), particle finite-element method (PFEM), finiteelement method with Lagrangian integration points (FEMLIP) and element-free Galerkin (EFG) (Soga et al. 2016).
The choice of the most appropriate methods for the analysis of triggering and propagation stages also depends on the objectives to be pursued. Particularly, for studies aimed at landslide susceptibility, hazard and risk assessment over large areas, it is necessary to contemplate the use of 3D models intended here, for simplicity, models capable of simulating motion along a 3D surface (McDougall 2017).
The majority of research papers focus on landslide triggering (e.g., Sorbino et al. 2010;Schilirò et al. 2016;Lizárraga and Buscarnera, 2019) or propagation (e.g., Pastor et al. 2009Pastor et al. , 2014Pirulli and Sorbino 2008), and one numerical model able to quantitatively study, at the same time, landslide triggering and propagation over large area is still not available.
Recently, some papers have combined two different models for analysing both triggering and propagation phases over large areas (e.g., von Ruette et al. 2016;Fan et al. 2017). However, these papers do not provide both the depth and velocity of the propagating mass, in different sections and in different time, necessary for a quantitative assessment of debris flow risk and for the design of risk mitigation structures.
In this context, the paper presents a physically based approach for quantitatively analysing debris flow triggering and propagation with the future aim to apply this approach in the whole study area or in areas with similar geo-environmental conditions for a quantitative assessment of debris flow susceptibility, hazard and risk zoning and the evaluation of debris flow height and velocity to be used as input data in computing impact forces needed in the design of risk mitigation structures.
The proposed approach is based on geological and geomorphological surveys and on geotechnical investigations carried out in a study area frequently affected by debris flows and consists in using the "Transient Rainfall Infiltration and Grid-Based Regional Slope-Stability Analysis" (TRIGRS) for predicting the shallow landslide source areas of debris flows and the "Smoothed Particle Hydrodynamics" (SPH) for reproducing the propagation stage.
These two methods are introduced for modelling debris flow problems over large areas because, despite the drawbacks that are inherent in the original assumptions of the models, both methods have been proven to be robust, efficient and reliable tools, throughout their application to practical cases worldwide, as reported by, e.g. Godt et al. (2008), Salciarini et al. (2006Salciarini et al. ( , 2017, Sorbino et al. (2010) and Ciurleo et al. (2019b) in the successful application of TRIGRS for landslide inception, and by, e.g. Basu et al. (2011 and Pastor et al. (2009Pastor et al. ( , 2014 in the successful application of SPH for landslide propagation.
The reliability of the combination of these two methods and its results have been assessed against the capability of the physically based approach to reproduce the triggering area, propagation path and depositional area of two debris flows, which occurred in 2001 and 2005 in the municipality of Scilla, Calabria, Southern Italy. These can be considered representative of the most typical landslides in the study area and surrounding zones.
Although it can be argued that other methods might have similar features, the main aim of the paper is to show that a combination of a limit equilibrium method (TRIGRS) and a meshless numerical method (GeoFlow SPH) is useful in simulating the mechanics of landslide motion during the failure, propagation and deposition stages in the frame of risk assessment and for designing risk mitigation structures.

Proposed methodology
The methodology proposed herein for the analysis of debris flows starts with the identification of shallow landslide events, in terms of classification, geometry and triggering rainfall data. It consists of three stages (stage I, II and III) and is based on a series process that considers the output of each previous stage as the input to the stage that follows (Fig. 1). Stage I consists of creating and compiling a spatial database, and aims to identify the range of variation of input parameters to be used in the second and third stages of analysis. Stage II, using as input data those derived from stage I, is employed to analyse debris flow triggering areas. Stage III consists of analysing the propagation path and deposition zone starting from the information gathered in the previous stages of analysis. To perform these analyses, in stages II and III, the use of physically based models is necessary in order to obtain significant results.
In detail, stage I starts with the collection of geological and topographic data available in the studied area. Referring to geotechnical and rheological data, stage I is pursued by combining and summing up the data regarding the soils involved in debris flows, available in scientific literature, with in situ investigations and laboratory tests specifically carried out for soils outcropping in the study area. In this stage, all geotechnical data needed for the reconstruction of soil cover thickness, index properties of investigated soils, shear strength, hydraulic properties and rheological properties of soils should be summarised and critically analysed. To do this, the authors suggest performing several in situ tests (e.g. continuous borehole drillings, seismic tomography, standard penetration tests, dynamic penetration tests) and several laboratory tests (e.g. classification tests, shear strength tests, hydraulic tests). Referring to the rheological characterisation of soil-water mixtures, performing laboratory tests (e.g. viscometer and/or rheometer tests) in order to identify rheological parameters is suggested. It is worth highlighting that the selection of a particular rheological model is difficult and several models can provide similar results . Basically, different models dedicated to the behaviour of geomaterials are available in scientific literature (Pastor et al. 2009) and the most commonly used are the frictional-turbulent "Voellmy" resistance and the visco-plastic Bingham-type resistance relationship.
The main aim of this stage is to identify both the range of variation of the main geotechnical properties of the soils affected by debris flows and the rheological law capable of simulating the soil water mixture behaviour. In this paper, the debris flow triggering stage (stage II) is studied by the Transient Rainfall Infiltration and Grid-Based Regional Slope-Stability model (TRIGRS). It is a deterministic distributed model, physically based, widely used for computing the triggering areas of rainfall-induced shallow landslides in different parts of the world. It couples an infiltration model with an infinite slope stability model assumption, generally verified because debris flows have shallow slip surfaces compared to their length. The infiltration model is based on the linearised solution of Richards equation (Iverson 2000;Baum et al. 2002) and the one-dimensional infinite-slope stability model to compute factor of safety on a cell-by-cell basis. During the analysis, the program makes it possible to vary the following input parameters: precipitation intensity, topographic data, soil cover thickness, initial water-table depth, saturated vertical hydraulic conductivity, saturated hydraulic diffusivity, a four-parameter soil-water characteristic curve (θ r , θ s , α, K sat ), cohesion for effective stress, shear strength angle and total unit weight of soil.

Original paper
In this stage, several parametric analyses have to be carried out varying the input parameters in the range identified in stage I.
In output, TRIGRS allows us both to identify unstable (FS ≤1) from stable cells (FS> 1) and to quantify the initial volume of debris flow in the source areas.
The best performance of the model is chosen on the basis of a dimensionless index (I trig ), for the triggering stage, defined by Ciurleo et al. (2019a) and reported as follows: where A TL are the landslide source areas according to the landslide inventory (observed source areas), and A UTL are the areas computed as unstable located within the A TL . The validation of the analysis, consisting of the best analysis of debris flow source areas, depends on the value assumed by I trig . Particularly, according to Ciurleo et al. (2019a), I trig lower than 80% cannot be accepted because a welldefined back analysis should be at least capable of predicting 80% of the debris flow source area. If I trig is less than 80%, the validation process is not satisfied and the geotechnical data used in the analyses must be verified (Fig. 1).
It should be noted that, when I trig is higher than 80%, stage II ends quantitatively providing the initial volume of debris flow source area (V 0 ), which can be considered one of the most important input parameters of stage III.
In stage III, the propagation stage is studied by the numerical Smoothed Particle Hydrodynamics (SPH) and allows the assessment of the final volume and the shape of debris fan, the propagation path, travel distance, the thickness and the velocity of flowing mass.
The SPH method is a fully mesh-free method where the absence of a computational grid or mesh allows us to easily model flows with a complex geometry or flows where large deformations occur.
In order to carry out the numerical analyses using the SPH code, the following input parameters are necessary: initial volume of the detached mass (V 0 ); discretisation of the initial mass; topographic data of the study area; rheological law and parameters; erosion law and parameters.
Moreover, numerical simulations can be performed using both the digital elevation model (DEM) and the digital surface model (DSM), with pixel size equal to 5 m, 2 m and 1 m (respectively 1:5000, 1:2000 and 1:1000 scale), for taking into account the influence of the topographic data and the presence of obstacles (Agliardi and Crosta 2003;McDougall 2017;Guo et al. 2020).
In order to calibrate landslide runout models, methods based on optimisation theory and statistics have been recently performed (Aaron et al. 2019).
For the calibration of models, it is worth defining quantities that can be used to compare simulation results to field data.
In the proposed methodology, we used the debris fan area mapped in the landslide inventory as quantity to validate the numerical analyses.
In this regard, when debris flow travel distance, propagation path, shape and area of debris fan are known, the validation of numerical analyses can be based on the dimensionless index (I dep ) defined for debris fan, as follows: where A TDF is the debris fan mapped in the landslide inventory and A SDF is the numerically computed debris fan located within the A TDF . Once debris flow is simulated, the values of debris flow height and debris flow velocity can be reconstructed in a significant section and can be used for evaluating the impact forces on the designed protection structure following standard or advanced procedures as proposed by Gioffrè et al. (2017).

Studied area and landslide phenomena
The study area is located in Favazzina, a hamlet of the Scilla municipality (Reggio Calabria, Italy). Two recent landslides were classified as debris flows that started as translational landslides involving the shallow material and weathered part of the metamorphic bedrock and occurred upslope to the channels. The mobilised masses began to move in steep channels and the bed became subject to soil erosion that caused an addition of further soil material to the flow.
The most significant landslide events recorded in the area occurred in May 2001 and March 2005, causing extensive damage and involving various lifelines (Fig. 2).
The 2001 and 2005 debris flows initiated from the slide from a steep bank and are characterised by the small volume of the initiating slide; the bulk of the material involved in the debris flow events originated from entrainment from the path.
Particularly, on 12 May 2001, two translational shallow landslides were triggered at the head of the Favagreca river, Fig. 2, respectively, at 567 m a.s.l. (triggering area A 01 ) and 558 m a.s.l. (triggering area A 02 ) above sea level in correspondence with two incisions. At about 300 m a.s.l., the two unstable masses joined, forming one main channel before affecting the SNAM station of the methane pipeline, the main SS 18 road and the railway causing the derailment of the Turin-Reggio Calabria intercity train.
On 31 March 2005, a similar phenomenon was activated in the valley near Favazzina. It was characterised by three translational shallow landslides, Fig. 2, at about 370 m (triggering area A 03 ), 242 m (triggering area A 04 ) and 170 m a.s.l. (triggering area A 05 ), respectively, that evolved into debris flow. This debris flow caused serious damage to the transport infrastructures (the SS18 state road and the railway) and the derailment of the ICN Reggio Calabria-Milan intercity train.
From a geometrical point of view, for the event on 12 May 2001, the only available information is that the triggering areas presented a prismatic geometry with a slip surface located at a depth of about 1.5 m (Bonavina et al. 2005).
With regard to rainfall data, the available information is daily precipitation from the Scilla rain gauge with reference to landslide triggering dates on 12 May 2001 and 31 March 2005.

Application of the proposed methodology
Stage I-topographic, geological, mechanical and rheological data A topographical analysis, which contemplates the use of geometrical parameters such as the area of the basin and the slope of the fan, allowed the authors to provide more information on the type of dominant flow for the studied drainage basins. This analysis, used in scientific literature for debris flow susceptibility mapping (Wilford et al. 2004;Bertrand et al. 2013;Church and Jakob 2020), consists in comparing the Melton Index with the slope of the fan representative of the dominant transport process. The Melton Index provides information about the gravitational energy of the basin as it represents the difference in elevation between the Referring to the available topographical maps, for a realistic simulation and assessment of debris flows, the DEM or DSM can be used. DEM and DSM, with pixel size equal to 5 m, 2 m and 1 m, were available for the Favazzina case study and they have been used to obtain both topographical and geometrical characteristics of the slope and the exposed elements.
The study area is characterised by a Palaeozoic crystalline basement that shows intense and deeply weathered conditions (Fig. 4). Particularly, according to Borrelli et al. (2012) and Gioffrè et al. (2016), highly and moderately weathered rocks (respectively class IV and III) crop out in the middle-lower portions of the slope (below 300 m a.s.l.) while completely weathered rocks (class V) prevail in the upper portion of the slope, above 300 m a.s.l. About 60% of the study area is covered by residual, colluvial and detrital soils (gneiss of class VI) which are the soils potentially affected by debris flows.
Finally, along the main streams through the area, gravellysandy loose alluvial deposits lie, while, between the sea and the base of the slope, beach deposits outcrop made up of sands and gravels.
It should be pointed out that the weathering grade map was drawn up after the 2001 and 2005 debris flows. In its present configuration, it shows that the triggering areas of the 2001 debris flow involved gneiss of class VI, soil still visible on site. On the contrary, in correspondence with the 2005 debris flow triggering areas, gneiss of class VI, V-IV and III outcrop, this is probably due to the effect of the 2005 debris flow that carried away the majority of the outcropping class VI.
From a tectonic point of view, according to Gioffrè et al. (2016), the study area is crossed by several fault segments mainly NE-SW and WNW-ESE. The main fault system, NE-SW trending, is arranged in a north-westward stepwise system that controls the morphology of the Favazzina slope. The second fault system, WNW-ESE trending, is more ancient than the previous one and morphologically less evident. This fault system is responsible for the formation of the hydrographic network of the study area (e.g. Favagreca channel) that flows towards the Favazzina coastal plain.
In the present study, in order to better characterise the soils involved in debris flows from a geotechnical point of view, a large site investigation campaign (Fig. 4) was carried out. It consisted of 5 continuous boreholes, 10 seismic refraction tomographies, the collection of 9 undisturbed cubic samples and 5 undisturbed samples. The site investigation campaign enriches the information already available in the study area, consisting of 3 continuous core drilling boreholes B1 (10) , B2 (10) and B3 (10) (Fig. 4(f)). Figure 4(f) shows that residual, colluvial and detrital soils (class VI) reach, in the upper part of the slope where the 2001 debris flow occurred, a depth ranging from 1.4 to 3 m; at greater depth boreholes prevalently show the presence of completely weathered gneiss (such as class V) up to 19 m of depth (B2 (10) ).
Starting from the above-mentioned data and in order to prevalently characterise residual, colluvial and detrital soils that are more susceptible to debris flow occurrence, the new site investigation campaign consisted of 5 continuous borehole drillings and Q5 samples as well as several Q3 samples were taken.
Regarding seismic tomographies, and in order to identify the transition between gneiss of class VI and less weathered gneiss, borehole logs (B1-B5) were overlapped with seismic prospects (Fig. 5).
All boreholes perfectly were located on the seismic prospects except borehole B4 that is projected on the section. In particular, borehole B1 shows a clear transition between gneiss of class VI and less weathered gneiss at a depth of 4 m, borehole B2 shows the same transition at a depth of 1.2 m and boreholes B3 and B5 show this difference at respectively 3.5 m and 3 m. In the abovementioned cases, p wave velocity assigned to class VI assumed values up to 500 m/s. Borehole B4, located at the head of the 2001 debris flow, shows a clear transition between gneiss of class VI and less weathered rock at a depth of 0.8 m (Fig. 5). This borehole is carried out upstream to the section L2T2 (Fig. 4a). Section L2T2 shows value of thickness equal to 1.8 m and this difference shows that the thicknesses of residual, colluvial and detrital soils vary on the slope direction.
Furthermore, borehole B4 is located at the same elevation of boreholes B1 (10) , B2 (10) and B3 (10) and especially borehole B2 (10) that is the deepest one shows that up to 19 m only class V prevalently has been found in the borehole. This aspect allows us to identify a p wave velocity equal to 2000 m/s as the transition between prevalently class VI-V to class V-IV.
P wave calibration carried out in available surveys allows us to identify the depth of class VI prevalently and the transition between the previously mentioned class and less weathered rocks (Fig. 6).   Fig. 6(b)), the depth of class VI ranges from 2 m to a few decimetres (L8T2, L9T1, L10T1), except for L8T1 that shows a depth of class VI up to 6 m.
Regarding section L10T1, this does not show class VI except for an insignificant decimetric thickness; this is compatible with its location on the map; indeed the map of weathering grade shows that, for now, class VI gneiss is not present in that position (Fig. 4).
For sections L8T1 and L2T2, the presence of a depth of 6 m and 8 m of class VI respectively is due to low slope gradients (for section L8T1) and to a morphological depression filled by soils (for section L2T2). This depression is located at a distance between 8 and 18 m on the section L2T2 and was filled by soils as a first mitigation measure after the 2001 debris flow. In order to reproduce debris flow occurrence, these two values have been neglected from the analysis.

Original paper
To sum up, gneiss of class VI has a depth of about 2 m or less than 2 m (as shown in most of the sections) in the study area.
Standard and advanced laboratory tests are underway on the specimens trimmed from undisturbed samples for the geotechnical characterisation of the soils making up the most superficial layer susceptible to debris flows. In particular, direct shear tests, triaxial tests in saturated and unsaturation conditions, the Richards pressure plate tests, traditional and controlled Fig. 5 Comparison between core appearance, interpreted weathering class and seismic velocity profiles suction oedometer tests, diffractometric analysis (XRF) and analysis using the scanning electronics microscope (SEM) are being carried out. On the Q3 samples, several laboratory tests as classification and mineralogical tests have been performed.
Mineralogical classification has been carried out by diffractometric analyses, the first two XRF tests on specimens taken from samples collected on boreholes B4 and B5 at depth ranging from 0.5 and 1.0 m. The obtained results show two perfectly overlapping curves with a high presence of quartz followed by feldspars (albite and microcline), phyllosilicates (biotite) and kaolinite that is a typical composition of weathered gneiss (Biondino et al. 2020).
The classification tests carried out on Q3 samples allow us to classify the soil (Fig. 7), according to the Unified Soil Classification System (USCS), as prevalently silty sand (SM); the fine fraction is prevalently not plastic and when it shows a modest plasticity, the values of LL range from 34.4 to 60.7% and PI from 7 to 15.3% and are located below the A line. Two samples, respectively B2_S2 and B5_S1 (Fig. 5), can be classified, according to USCS, as clay sand (SC), the fine fraction is plastic, the values of LL range from 34.0 to 41.9% and PI from 11 to 16% and are located above the A line. Only one sample located in the borehole B2 at a depth ranging from 0.6 to 1 m, can be classified, according to USCS, as organic clay with high medium plasticity, and it shows a value of LL equal to 60% and PI equal to 26%.
The grain size distributions (in red in Fig. 7) show a wide range of variation from prevalently "coarse" to prevalently "fine" material.
Moreover, in order to evaluate the percentage in number of particles rather than in weight, the procedure proposed by Moraci et al. (2012) has been applied. In particular, each grain size distribution has been discretised by N pairs of diameters D i (D 1 , …, D i , …, D N ) and for each pair the corresponding frequency in the mass ΔP mi has been obtained. The discretised soil mass grain size distribution is divided into two parts and the soil numerical percentage grain size distribution, characterised by D i and ΔP ni values, is obtained from the soil weight percentage grain size distribution through the relation proposed by Musso and Federico (1983), Eq. 3, and considering that the specific gravity is the same for all grains: where ΔP ni is the numerical frequency of the particle fraction with average diameter D i . The grain size distributions in number (in black in Fig. 7) rather than by weight show that all curves have a value of cumulated numerical frequency P (%) which falls within silt and clay and finally the numerical frequency of clay oscillates from about 30 to 90 % (Fig. 7).
This indicates that the fine matrix constitutes the solid skeleton and coarse particles are inclusions inside the solid skeleton. For this reason, the rheological behaviour of the solid water mixture has been characterised by using viscometer laboratory tests on the fine fraction of the soil.
These soils have a value of saturated unit weight (γ sat ) varying from 19 to 22 kN/m 3 : void index (e) variable from a minimum of 0.65 to a maximum of 1.15 and soil porosity (n) varying from 0.4 to 0.55 (Antronico et al. 2010).
Since the laboratory campaign aimed at characterising class VI soils under total and partial saturation conditions has not yet been completed, a compilation of data available in literature on the same or similar formations has been used for the mechanical and hydraulic characterisation.
Referring to shear strength properties, the data available from literature show a range of variation of effective cohesion from 0 to 5 kPa and a shear strength angle from 30°to 44° (Antronico et al. 2010;Schilirò et al. 2015).
Regarding hydraulic properties, the authors considered the values of saturated conductivity (K sat ), ranging from 1.3E-06 m/s to 6.6E-05 m/s, and saturated volumetric water content θ s , ranging from 0.38 to 0.4, obtained by Cascini et al. (2006) and Schilirò et al. (2015) for gneiss similar for genesis and stress history with those studied.
Referring to the rheological behaviour of the soil water mixture, the Bingham rheological model has shown considerable potential also in the study of the debris flow which occurred in the Favazzina area .
In order to verify the rheological behaviour of the most superficial soils involved in the 2001 and 2005 debris flows, the data already available were collected and further laboratory tests were carried out with a rotary viscometer on mixtures consisting of the fine fraction of these soils and water in different solid concentration by volume.
The tests performed by Moraci et al. (2017) on these soils investigated a range of solid concentration by volume ranging from 20 to 45%; further tests have been carried out by varying the solid-liquid volumetric concentration from 30 to 45%, for two different soil samples (Rs1 and Rs2).
Viscometer tests have been interpreted using the Bingham rheological law because it is the most efficient in the simulation of debris flows in the study area Ciurleo et al. 2020).
The equation at the base of the model is reported as follows: where τ 0 (Pa) is the yield stress and μ is the dynamic viscosity (Pa s). It is generally accepted that yield stress and dynamic viscosity exponentially increase by increasing the solid concentration in volume Cv according to laws listed below: Both values are linked to two empirical coefficients α i and β i obtained by a regression analysis on experimental data for each one of the tested specimens.
The data interpretation shows that Eqs. 5 and 6 well interpolate the variation of τ 0 (Pa) and μ (Pa s) varying the solid concentration by volume of the mixture (Fig. 8). Particularly, the regression analysis has provided the following values for Rs2: α 1 =0.0002 Pa and β 1 =0.27; α 2 =0.001 Pa s and β 2 =0.21.
Therefore, the previous equations are: while with reference to sample Rs1, the obtained data can be well located on the law of τ 0 -C v and μ-C v already identified by Moraci et al. (2017) with the following values α 1 =0.251 Pa and β 1 =0.132; α 2 =0.0112 Pa s and β 2 =0.163, respectively (Fig. 8).
Therefore, the equations become: Landslides 18 & (2021) 2081 In both cases, using Eqs. 7-10, it is also possible to obtain τ 0 and μ for other values of solid concentration by volume Cv.

Stage II-triggering analyses
The input data of TRIGRS are as follows: the digital elevation model (DEM) with pixel size of 5 m; soil cover thickness that, according to in situ investigations and the available on site evidence immediately after debris flows, has been considered equal to 1.5 m; geotechnical properties (hydraulic and mechanical) of class VI weathered gneiss, identified in stage I and summarised in Table 1.
The saturated hydraulic diffusivity (D 0 ) was calculated according to Grelle et al. (2014) and Schilirò et al. (2015) using the formula below: where K s is the saturated hydraulic conductivity, H the average soil thickness (assumed equal to 1.5 m where the 2001 and 2005 debris flows occurred) and Sy is the specific yield that, according to Johnson (1967), and Loheide et al. (2005), has been assumed equal to 0.34 for the analysed soils. Regarding pore water pressure regime, due to the lack of data and, according to Ciurleo et al. (2019b), the water table was assumed at the ground surface in the upper part of the slope (zone "A" in Fig. 9) and at the contact between class VI and less weathered gneiss in the remaining part of the study area (zone "B" in Fig. 9). This assumption derived from the results obtained by Ciurleo et al. (2019b), demonstrating that if the water table is not considered located at 0 m from ground surface in zone "A" (a buffer zone constituted by three contour lines, each equal to 5 m, from the secondary road), one of the triggering source areas (triggering area A 01 , Fig. 9) which occurred in 2001 cannot be analysed by the model. Several parametric analyses in saturated conditions have been performed by TRIGRS using different values of shear strength data. The main aim is to identify the combination of geotechnical parameters to best fit the 2001 and 2005 debris flow occurrence.
The overall analysed cases, Table 2, have been implemented considering the average values of hydraulic parameters (K sat = 1.79E−05 m/s, D 0 = 7.92E−05 m 2 /s, and θ sat = 0.4) and a constant value of soil cover thickness equal to 1.5 m both gathered in stage I.
With reference to shear strength properties, starting from the range of values gathered in stage I (Table 1), five parametric analyses (Table 2) for each landslide event (12 May 2001 and 31 March 2005) have been carried out considering (i) the average value of shear strength angle (φ′=38°) and varying the cohesion value between 1 and 2.5 kPa (average value), and (ii) the average value of cohesion (2.5 kPa) and varying the shear strength angle from 30°to 38°(average value).
The results obtained for the 2001 and 2005 landslide events (Table 3) show that for both cases, (i) the best fit between the source areas which occurred on 12 May 2001 and on 31 March 2005 and the numerical analyses was achieved for case 4 respectively obtaining a I trig (2001) = 90.5% and I trig(2005) = 94.9%, (ii) in these cases, for the 2001 debris flow, the computed triggering volumes are 900 m 3 and 1125 m 3 (triggering volumes V 01 and V 02 , respectively, in Fig. 10), for the 2005 debris flow, the computed triggering volumes are 325 m 3 , 187 m 3 and 761 m 3 (triggering volumes V 03 , V 04 , and V 05 in Fig. 10).

Stage III-propagation analyses
In the third stage, the SPH model used as input data the triggering volumes identified in stage II, topographical data, rheological parameters and erosion rate (Table 4), considering the erosion law of Hungr (1995) by the "growth rate" (Es) that represents the bed-normal depth eroded per unit flow depth and unit displacement, defined as: where V fin is the final volume of the mobilised material, V 0 is the initial volume of the debris flow and d is the travelled distance by debris flow. Table 4 summarises the values of τ 0 and μ b obtained by using Eqs. 7-10 considering value of Cv ranging from 50 to 55%.
Stage III consisted of an iterative analysis using different rheological parameters (function of Cv), erosion coefficients (Es) and topographical data (DEM and DSM). In particular, the authors started from a DEM with pixel size equal to 5 m for coherence with the topographical data used by TRIGRS, and performed several analyses using DSM and DEM (with pixel size equal to 1 m, 2 m and 5 m) and varying Cv and Es.
The 2005 debris flow is a landslide phenomenon characterised by three triggering masses moving on the same track and forming, at the end of the propagation path, a triangular depositional fan.
For this debris flow, the obtained numerical results show that, independently on the scale of analysis, and for fixed values of Cv and Es, the analyses carried out with DSM always show both the best reproduction of the triangular shape of the 2005 debris fan (Fig. 11) and higher values of I dep than those performed with DEM. Furthermore, with reference to the scale of analysis, it is evident that DSM at 1:5000 scale overestimates the propagation path ( Fig.  11(a)) that appears significantly larger than the real propagation zone. As expected, this overestimation tends to be reduced passing from 1:5000 scale to 1:2000 scale, where the propagation path appears more coherent with the real propagation zone and I dep assumes the highest value equal to 78.42% if DSM is considered (Fig. 11(b)). This value is higher than those obtained considering a DSM with pixel resolution of 1 m (Fig. 11(c)). Moreover, the deepest deposit formed at the end of the propagation path in correspondence of the depositional area ( Fig. 11(a-c)) shows values significantly different passing from a topographical scale 1:5000 to 1:1000. Figure 12 summarises the results obtained in Fig. 11 and clearly shows that I dep obtained by using DSM data (continuous line) is higher than that obtained by DEM data (dotted line).
Focusing only on DSM, and in order to verify that a pixel resolution of 2 m is better than 1 m, several numerical analyses have been carried out varying the value of Cv and fixing Es (Fig.   13). These numerical analyses confirmed that, for varying Cv, the topographical data at 1:1000 scale implies that some smoothed particles come out of the propagation path and this effect is more evident by reducing the value of solid concentration by volume Cv. Indeed, a value of Cv equal to 55% does not produce overestimation of the propagation path independently on the DSM scale. On the contrary, a Cv equal to 55% does not allow us to reproduce the debris fan presenting an I dep changing from 44.23 to 45.34% respectively for DSM with pixel resolution of 1 m and DSM with pixel resolution of 2 m.
Focusing only on DSM with pixel resolution of 2 m, several analyses have been carried out varying the values of Cv and Es (Fig. 14). The highest value of I dep is obtained by considering a value of Cv = 50% and Es = 0.0018 m −1 .
The 2001 debris flow is a landslide phenomenon characterised by two triggering masses that started from two separate morphological hollows and joined at an altitude of 300 m a.s.l. continuing the propagation along the same track and, at the end of the propagation path, formed a rectangular depositional fan.
Starting from the consideration gathered by the 2005 debris flow analyses, the only used topographical data is the DSM at 1:2000 scale (Fig. 15). The figure shows that the highest value of I dep equal to 72.0% is obtained by considering a value of Cv=54% and Es= 0.0015 m −1 .

Discussion of the results
With reference to stage I, the creation of a database, the overlapping between seismic tomographies and borehole logs allowed the reconstruction of a clear transition between gneiss of class VI and less weathered rocks at a depth ranging from a few decimetres to less than 2 m where p wave velocity assumed values up to 500 m/s. Furthermore, the deepest borehole B2 (10) showed that up to 19 m only class V prevalently has been found. At this depth, seismic refraction prospects carried out in the study area measured a p wave velocity equal to 2000 m/s according to scientific literature (Hunt 2006).
Regarding grain size distribution, significant information has been provided when passing from grain size distribution by weight to number distribution. Indeed, in the first case, grain size distribution showed a wide range of variation from prevalently "coarse" to prevalently "fine" material, while for grain size distribution by number it was found that all curves have a value of cumulated numerical frequency P (%) falling within silt and clay fraction. For the propagation analyses, the high presence of fine fraction allowed us to use viscometer tests on fine fraction and the obtained results were representative of rheological behaviour prevalently of solid water mixture.
The proposed methodology can be used both if an accurate in situ and laboratory geotechnical characterisation is available and if a mechanical and hydraulic characterisation from scientific literature already exists. Clearly, the latter is possible only after checking the consistency between the soils investigated in scientific literature and those studied.
Referring to rheological behaviour of the soil water mixture, the viscometer tests carried out on the specimen from Rs1 provided values of τ 0 and μ different to those obtained for the specimen from Rs2. This is probably due to the clay fraction (CF). Indeed, the specimen from Rs1 is made up of a silt fraction = 31.4% and a clay fraction = 54.2% while the specimen from Rs2 shows values of  Regarding the Bingham rheological model, it was observed that the increase in CF produces higher values of yield stress when compared with the sample with smaller CF, but this trend is not as clear when dynamic viscosity values are considered (Fig. 16).
Referring to stage II of the proposed methodology, the validation index I trig showed that the best simulation of the 2001 and 2005 debris flow triggering areas (I trig higher than 90%) has been obtained for a value of effective cohesion equal to 2.5 kPa and a value of shear strength angle equal to 30°.
For stage III, the numerical model used appears sensitive to the roughness of topographical data, indeed if all other factors are equal, in agreement with McDougall (2017), a rougher surface derived from a detailed topographical data (e.g. DSM with a pixel  Original paper resolution of 1 m) leads to a shorter modelled runout distance and, in some cases, a very rough sliding surface can also destabilise depth-averaged flow models. To overcome this problem, in SPH, the input topographical data has been filtered in order to smooth out the roughness and, fixing other parameters, the best obtained results have been gathered considering a DSM with pixel resolution of 2 m that is the most suitable for the analysed debris flows. For the 2001 and 2005 debris flows, since the SPH model considered the value of Cv constant for the whole debris flow, that is actually subjected to entrainment of material and water during the downhill path, it was decided to use different values of Cv in order to identify the equivalent value that best simulates the 2001 and 2005 debris flows in terms of I dep . I dep assumed values higher than 90% for the 2005 debris flow and higher than 70% for the 2001 debris flow, respectively, considering Cv=52% and Es=0.0018 m −1 and Cv=54% and Es=0.0015 m −1 .
For the 2005 debris flow, if topographical data and growth rate (Es) are kept constant, the I dep value tends to decrease as the solid concentration by volume Cv increases; and fixing Cv, the I dep value undergoes a slight increase as the growth rate increases (Fig. 14); while for the 2001 debris flow, considering DSM with pixel resolution of 2 m, the I dep prevalently assumes higher values increasing growth rate (Es) and, fixing growth rate, I dep shows a variable trend as Cv varies (Fig. 15). The variation of debris flow height and velocity in time is sensitive to different aspects linked to the distance (from the considered section), the volume of detached mass and the geometry of the propagation path.
Indeed, on the one hand, for the 2005 debris flow, characterised by three different detached masses with different volumes (Fig. 10), debris flow height and velocity reached their maximum value at 6 s (time of first impact) and rapidly decreased, due to the position of the O-O′ section located at a distance of about 130 m from the lowest detached mass with the highest volume and influenced by a channelised flow (Fig. 17). On the other hand, for the 2001 debris flow (Fig. 18), the maximum velocity was reached after 50 s, the   detached masses (with about the same volume) are located at a distance over 1000 m, and the maximum height was not reached at 50 s (time of first impact) but at about 130 s. This is due to the fact that in this case the O-O′ section is located after a flat zone; it is not directly impacted by a channelled flow and part of the energy of the debris flow was dissipated in the flat upstream area.
By comparing Figs. 17 and 18, it is clear that in both cases the maximum velocity corresponds to the first impact velocity while the maximum debris flow height assumes different values depending on the topography of the propagation path and the location of the O-O′ section as reported above.
In both cases, the obtained values in terms of debris flow height and velocity are coherent with those numerically computed by Pastor et al. (2009) for other debris flows phenomena in other basins.

Conclusions
In order to quantitatively assess debris flow susceptibility and hazard in the frame of risk management, a new methodology based on a combination of limit equilibrium method over large areas, for the triggering phase, and a numerical model for the propagation analyses, is proposed.
The first stage of the methodology was the creation of a database by merging the data already available in the studied area with those obtained by in situ and laboratory tests. Stage II, using as input data those derived from stage I, was employed by means of TRIGRS used to identify triggering volumes that represent some of the most important input data for stage III (propagation analyses). Numerical simulations of the propagation path and deposition zone have been performed by SPH using as input data the initial triggering volumes gathered in stage II, different topographies and rheological data.
The obtained results highlighted the role played by in situ and laboratory tests that allowed us to better characterise the outcropping soils, identifying the transition between class VI and less weathered gneiss, and reasonably assert that, for the propagation analyses, the results of viscometer tests carried out on fine fraction can be considered representative of the rheological behaviour of the soil water mixture. Furthermore, a combination of viscometer tests with grain size distribution of the two investigated samples allowed us to show that increasing the clay fraction the rheological behaviour of the soil is affected by two empirical coefficients (a i and β i ) of the Bingham rheological law.
Therefore, a sound prediction of the propagation path and depositional area must include at least the knowledge of the grain size distribution and the rheological behaviour of the soil water mixture.
The parametric analyses also allowed us to identify both the values of cohesion and shear strength angle for which failure occurs (c′=2.5 kPa, ϕ′=30°) and the range of variation of the solid concentration by volume Cv = 50 -54% and the growth rate Es = 0.0015 -0.0018 m −1 to be used for the susceptibility and hazard assessment over the entire study area.
Moreover, the proposed methodology enabled us to obtain, in any section of the propagation path, debris flow height and velocity in time that represent the principle information for designing risk mitigation structures through protective structures (passive) in relation to dynamic and static stresses. The design of protective structures requires the calculation of the impact force, function of several variables including debris flow height and velocity. The latter are the characteristics of debris flow phenomena necessary for the evaluation of the impact force that, according to Gioffrè et al. (2017), can be calculated by adding the static component, which depends only on the flowing mass height, to the dynamic component, that depends on both the flowing mass height and its velocity.
Due to the great potential of the proposed methodology, both in terms of zoning landslide hazard and for the design of risk mitigation works, the authors suggest its use in contexts different from the one in which it was tested. Clearly, in order to successfully apply the formulated methodology, a set of robust and reliable input data is required for its application to other case histories.