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

Abstract

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 (I_trig and I_dep). 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 I_trig were higher than 90% while I_dep 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. 2016, 2017; Ferlisi 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 history-dependent 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), finite-element 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. 2009, 2014; Pirulli 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. (2006, 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), Cascini et al. (2014) and Pastor et al. (2009, 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.

Fig. 1

Flow chart of the proposed methodology

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 (Moraci et al. 2017). 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.

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:

$$ {I}_{trig}=\frac{A_{UTL}}{A_{TL}}\times 100 $$

(1)

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 well-defined 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₀), 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₀); 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:

$$ {I}_{dep}=\frac{A_{SDF}}{A_{TDF}}\times 100 $$

(2)

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).

Fig. 2

2001 and 2005 debris flows in Favazzina

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₀₁) and 558 m a.s.l. (triggering area A₀₂) 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₀₃), 242 m (triggering area A₀₄) and 170 m a.s.l. (triggering area A₀₅), 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 maximum elevation of the basin and the elevation of the fan apex divided by the square root of basin area.

Figure 3 shows the empirical relationship proposed by Bertrand et al. (2013) which, based on a database of 620 basins subject to flow phenomena, separates fluvial processes (hyper-concentrated flows, solid transport, etc.) from the landslide phenomena of debris flow type by means of the green curve.

Fig. 3

Classification of the type of flow (modified from Bertrand et al. 2013)

With reference to the geometries of the basin affected by the 2001 and 2005 debris flows, the obtained points shown in Fig. 3 confirm that in both cases they are basins affected by debris flows.

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.

Fig. 4

Weathering grade map (Gioffrè et al. 2016 mod.) and location of in situ surveys. Geological legend: 1, gneiss of class II; 2, gneiss of class III; 3, gneiss of class IV; 4, gneiss of class V; 5, colluvial and detrital soils (class VI); 6, coastal and alluvial deposits (Holocene); 7, terraced marine deposits (Middle Pleistocene); 8, marine coarse sandstone deposits (Upper Pliocene-lower Pleistocene); 9, landslide debris (Holocene); 10, normal fault; 11, fault with undetermined kinematics. In situ campaign 2018, C1–C9 cubic samples; seismic refractions (L_iT_i): (a) detailed view of the area of borehole B4 and seismic refractions L1T1, L1T2, L2T1 and L2T2; (b) detailed view of the area of borehole B5 and seismic tomography L7T1; (c) detailed view of the area of borehole B4 and seismic tomography L4T1; (d) detailed view of the area of borehole B3 and seismic tomographies L8T1 and L8T2; (e) detailed view of the area of borehole B2 and seismic tomographies L9T1 and L10T1; (f) borehole logs of previous in situ surveys 2010

Finally, along the main streams through the area, gravelly-sandy 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₍₁₀₎, B2₍₁₀₎ and B3₍₁₀₎ (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₍₁₀₎).

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).

Fig. 5

Comparison between core appearance, interpreted weathering class and seismic velocity profiles

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 above-mentioned 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₍₁₀₎, B2₍₁₀₎ and B3₍₁₀₎ and especially borehole B2₍₁₀₎ 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). Figure 6(a) shows that at an elevation between 600 m and 500 m a.s.l (typical of 2001 debris flow source areas) class VI (p wave velocity up to 500 m/s) is about 2 m (L1T1) or less than 2 m (L1T2 and L2T1) deep, except for L2T2 that shows a depth of class VI up to 8 m. On the other hand, at an elevation varying from about 400 m a.s.l to 170 m a.s.l. (typical of 2001 debris flow source areas, 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.

Fig. 6

Seismic refraction tomographies (SRT) carried out in the Study area. (a) SRT performed at elevation from 600 m to 500 m a.s.l. (typical of 2001 debris flow source areas), (b) SRT performed at elevation from 400 m to 170 m a.s.l. (typical of 2005 debris flow source areas)

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.

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 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.

Fig. 7

Comparison between soil weight percentage grain size distribution and numerical frequency of the particle forming soil

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₁, …, 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:

$$ \varDelta {\mathrm{P}}_{\mathrm{ni}}\left({D}_i\right)=\frac{\frac{\Delta {\mathrm{P}}_{\mathrm{mi}}}{D_{\mathrm{i}}^3}}{\sum_{i=1}^N\frac{{\Delta \mathrm{P}}_{\mathrm{mi}}}{D_i^3}} $$

(3)

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³: 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 (Moraci et al. 2017).

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 (Moraci et al. 2017; Ciurleo et al. 2020).

The equation at the base of the model is reported as follows:

$$ \tau ={\tau}_0+\mu \frac{du}{dy} $$

(4)

where τ₀ (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:

$$ {\tau}_0={a}_1\ {e}^{\beta_1{C}_v} $$

(5)

$$ \mu ={a}_2\ {e}^{\beta_2{C}_v} $$

(6)

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 τ₀ (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: α₁=0.0002 Pa and β₁=0.27; α₂=0.001 Pa s and β₂=0.21.

Fig. 8

Trend of yield stress τ₀ and dynamic viscosity μ versus solid concentration in volume Cv

Therefore, the previous equations are:

$$ {\tau}_0=0.0002.{e}^{0.27{C}_v} $$

(7)

$$ \mu =0.001.{e}^{0.21{C}_v} $$

(8)

while with reference to sample Rs1, the obtained data can be well located on the law of τ₀-C_v and μ-C_v already identified by Moraci et al. (2017) with the following values α₁=0.251 Pa and β₁=0.132; α₂=0.0112 Pa s and β₂=0.163, respectively (Fig. 8).

Therefore, the equations become:

$$ {\tau}_0=0.251.{e}^{0.132{C}_v} $$

(9)

$$ \mu =0.0112.{e}^{0.163{C}_v} $$

(10)

In both cases, using Eqs. 7–10, it is also possible to obtain τ₀ 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.

Table 1 TRIGRS input data

The saturated hydraulic diffusivity (D₀) was calculated according to Grelle et al. (2014) and Schilirò et al. (2015) using the formula below:

$$ {D}_0=\frac{K_SH}{S_y} $$

(11)

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₀₁, Fig. 9) which occurred in 2001 cannot be analysed by the model.

Fig. 9

Initial water table locations. Legend: A = 0 m from the ground surface; B = 1.5 m from the ground surface (Ciurleo et al. 2019b mod.)

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₀ = 7.92E−05 m²/s, and θ_sat = 0.4) and a constant value of soil cover thickness equal to 1.5 m both gathered in stage I.

Table 2 Shear strength properties used in parametric analyses

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³ and 1125 m³ (triggering volumes V₀₁ and V₀₂, respectively, in Fig. 10), for the 2005 debris flow, the computed triggering volumes are 325 m³, 187 m³ and 761 m³ (triggering volumes V₀₃, V₀₄, and V₀₅ in Fig. 10).

Table 3 Results of analyses for the 12^th May 2001 and 31^st March 2005 debris flows

Fig. 10

The best back-analysed landslide inception results obtained by TRIGRS

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:

$$ {E}_s=\frac{\ln \Big({V}_{fin}/{V}_{0\Big)}}{d} $$

(12)

Table 4 SPH input data for Cv ranging from 50 to 55%

where V_fin is the final volume of the mobilised material, V₀ is the initial volume of the debris flow and d is the travelled distance by debris flow.

Table 4 summarises the values of τ₀ 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.

Fig. 11

SPH simulations for 2005 debris flow obtained by varying topographical data for the same Cv = 52% and Es = 0.0018 m¹

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).

Fig. 12

Relation between the scale of analysis vs I_dep for two different topographical data (DSM and DEM)

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.

Fig. 13

SPH simulations for 2005 debris flow obtained by varying topographical data (DSM 1:1000 and 1:2000 scale), the solid concentration by volume Cv and fixing Es= 0.0018 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⁻¹.

Fig. 14

SPH simulations for 2005 debris flow obtained by varying the solid concentration by volume Cv and erosion rate Es and fixing a topographical data (DSM 1:2000). Red box shows the best simulation in terms of I_dep

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⁻¹.

Fig. 15

SPH simulations for 2001 debris flow obtained by varying the solid concentration by volume Cv and grow rate Es and fixing a topographical data (DSM 1:2000). Red box shows the best simulation in terms of I_dep

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₍₁₀₎ 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 τ₀ 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 silt and clay fractions equal to 30.8% and 12.8%, respectively. 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).

Fig. 16

Values of Cv vs τ₀ e di Cv vs μ for Rs1 and Rs2 specimens

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 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⁻¹ and Cv=54% and Es=0.0015 m⁻¹.

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).

For the best simulations of the 2005 and 2001 debris flows, corresponding respectively to Cv=50% and Es=0.0018 m⁻¹ and Cv=54% and Es=0.0015 m⁻¹, and considering a representative section O-O′ located upstream of the SS18 road (the first exposed element affected by debris flows), the trends of debris flow height and debris flow velocity versus time are illustrated in Figs. 17 and 18.

Fig. 17

The best SPH simulation of 2005 debris flow and related trends of debris flow height and debris flow velocity versus time at the centre of O-O′ section

Fig. 18

The best SPH simulation of 2001 debris flow and related trends of debris flow height and debris flow velocity versus time at the centre of O-O′ section

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⁻¹ 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.