Derivation of critical rainfall thresholds for shallow landslides as a tool for debris flow early warning systems

Real-time assessment of debris-flow hazard is fundamental for developing warning systems that can mitigate risk. A convenient method to assess the possible occurrence of a debris flow is to compare measured and forecasted rainfalls to critical rainfall threshold (CRT) curves. Empirical derivation of the CRT from the analysis of past events’ rainfall characteristics is not possible when the database of observed debris flows is poor or when the environment changes with time. For debris flows and mud flows triggered by shallow landslides or debris avalanches, the above limitations may be overcome through the methodology presented. In this work the CRT curves are derived from mathematical and numerical simulations, based on the infinite-slope stability model in which slope instability is governed by the increase in groundwater pressure due to rainfall. The effect of rainfall infiltration on landside occurrence is modelled through a reduced form of the Richards equation. The range of rainfall durations for which the method can be correctly employed is investigated and an equation is derived for the lower limit of the range. A large number of calculations are performed combining different values of rainfall characteristics (intensity and duration of event rainfall and intensity of antecedent rainfall). For each combination of rainfall characteristics, the percentage of the basin that is unstable is computed. The obtained database is opportunely elaborated to derive CRT curves. The methodology is implemented and tested in a small basin of the Amalfi Coast (South Italy). The comparison among the obtained CRT curves and the observed rainfall amounts, in a playback period, gives a good agreement. Simulations are performed with different degree of detail in the soil parameters characterization. The comparison shows that the lack of knowledge about the spatial variability of the parameters may greatly affect the results. This problem is partially mitigated by the use of a Monte Carlo approach.


Introduction
Rainfalls with peculiar characteristics of intensity and duration may trigger debris flows (DF).These events are particularly dangerous for a number of reasons.DF may travel for long distances, as water flows and the steep slopes induce high velocities.They may easily impact the alluvial fans of mountain torrents that, especially in recent years, have been intensively urbanized.The density effect on the impact force of the DF makes it much higher than the one due to water flows and thus damages to buildings are proportionally greater.Predicting DF is quite difficult because there are no premonitory signs and the lag time between the occurrence of the triggering rainfall and the impact of the DF in the vulnerable downstream area is very short.
The mitigation of DF hazards may be achieved building structural countermeasures, such as check dams and retention basins.However, in some cases, the steep slopes of the interested areas, or the lack of space, makes it difficult to build structural countermeasures.Moreover, the impact of these works on the landscape may be rather high.This problem is particularly noted in areas with high environmental and historical value.For the aforementioned reasons, in many cases, non-structural countermeasures, such as warnings through real-time hazard assessment and civil protection measures are more suitable in reducing the risks.Economic Published by Copernicus Publications on behalf of the European Geosciences Union.
reasons can also influence the choice, as non-structural countermeasures are less expensive than structural ones.Due to the short lag time, warning systems must rely on forecasted and nowcasted rainfall.The warning is given when the forecasted rainfall overcomes a critical threshold.
The most common approach adopted in literature for the assessment of the critical rainfall thresholds (Caine, 1980;Wieczorek and Glade, 2005;Guzzetti et al., 2008;Brunetti et al., 2010) is based on the elaboration of data sets of historical events.The first limit of these approaches is due to the fact that they may only be adopted for the basins where a certain amount of inventoried DF events is available for the derivation of the empirical threshold line.Another drawback of empirically derived thresholds is that they cannot anticipate how DF hazards may change in response to changing environments, for example, land use changes, large forest fires' occurrence and sediment availability's fluctuation.The latter may have been determined, for example, by a past debris flow.
In order to overcome these limitations, it is possible to estimate rainfall critical thresholds through a model that reflects the physics of the phenomenon and provides the link between rainfall and possible DF hazards.
A DF can be triggered by several mechanisms.For example, triggering factors may be rainfalls, earthquakes, rapid snow melting, volcanic eruptions, sudden release of water from a lake, or a combination of these factors.In case of rainfall-triggered DF, two main types of formation mechanisms can be distinguished.The first one is the result of a progressive entrainment of bed sediments into a water flow (e.g.Takahashi, 1991;Berti and Simoni, 2005).The second type happens when a rainfall-triggered landslide evolves into a proper flow causing a DF or a mud flow (e.g.Johnson, 1984;Iverson et al., 1997;Baum et al., 2010).Landslidetriggered DF may involve higher sediment volumes and consequently are more dangerous than those resulting from continuous erosive processes (Malet et al., 2005).The study here presented is addressing only to this second type of triggering mechanism.Following the classification proposed by Hungr et al. (2001) the proposed tool is suitable for DF and mud flows triggered by shallow landslides or debris avalanches.
Theoretical models of rainfall-triggered landslides are usually based on the infinite-slope stability analysis in which slope instability is governed by the increase in groundwater pressure due to rainfall.These models are usually implemented in discrete landscape cells and compute the security factor for each cell.Many of the approaches proposed in literature are based on the hypothesis of steady groundwater flow conditions (Montgomery and Dietrich, 1994).Casadei et al. (2003) coupled an infinite slope stability model with a dynamic hydrologic model inspired by the Topmodel (Beven and Kirkby, 1979).Simoni et al. (2008) proposed a model (GEOtop-FS) that computes soil moisture and matric suction within soil layers by numerically integrating the Richards equation in a 3-D scheme.Additionally, Bathurst et al. (2006) developed another hydrological model (SHETRAN model), including an infinite slope stability module.Lanni et al. (2012) proposed a new hydrological model in which a dynamic topographic index is used to describe the transient lateral flow.In this model the lateral flow initiates only if hydrological connectivity is determined by rainfalls exceeding a threshold value.
A simplified model was proposed by Iverson (2000) to assess short-term pore water response to rainfall in the hypothesis of vertical infiltration.Baum et al. (2010) developed a model of the infiltration process using a two-layers system that consists of an unsaturated zone above a saturated zone, and implemented this model in a geographic information system (GIS) framework.Comparison of model results with observed scars of DF formation areas (Godt et al., 2008) has shown that the model of Iverson (2000) is more effective for regional shallow landslide hazard maps.However, the number of false positives and false negatives in the predicted unstable cells is still high.The time needed to reach potential instability in the unsaturated zone for a given rainfall was investigated by Godt et al. (2012).Examining four hypothetical hillslopes, they showed how the timing of potential slope instability depends on the soil water characteristics of the hillside's material.
The models that compute the safety factor, cell by cell, usually overestimate the potential instabilities because a single instable element is not going to move if it is surrounded by stable elements.In order to overcome this limitation, some authors (Lehmann and Or, 2012) developed techniques for the simulation of the cascade load redistributions: initially localized failures evolve into successive failures propagating across the hillslope.
The long computational times may prevent the operational use of distributed slope stability models for the prediction of the possible occurrence of a DF.In addition, if quantitative precipitation forecasting ensembles are used, the computational runs should be multiplied.Some studies are reported in literature in which the slope stability models are run offline for the derivation of critical rainfall threshold (CRT) for the triggering of landslides.Using a physically based model, Frattini et al. (2009) provided an estimate of the percentage of potentially unstable areas in which failure could be triggered with a certain probability.De Vita et al. (2012) derived the CRT simulating the effect of different rainfalls on a specific slope that was representative for the local topographic and geological conditions.Stability modelling was based on the simplified Bishop method, and the water pressures were computed solving the 2-D form of the Richards equation.
With the aim of giving DF warnings, it is not necessary to know the distribution of instable elements along the basin but only if a DF may affect the vulnerable areas in the valley.The capability to reach the downstream areas depends on many factors linked with the flow discharge, the topography, the solid concentration, the rheological properties of the debris mixture, the occurrence of liquefaction of the sliding mass, the entrainment of sediments from the bed and banks of the channel.Many of these factors are not time dependent.The most important rainfall dependent factors are flow discharge and correlated total debris volume.In the present study, the total volume that is unstable, and therefore available for the flow, is considered to be the governing factor from which it is possible to assess whether a DF will affect the downstream areas or not.
Many studies demonstrated the importance of the amplification of the debris volume due to sediment entrainment from the bed and bank of the channel during the runout process (e.g.Iverson et al., 1997;Takahashi, 2009;Papa et al., 2004;Quan Luna et al., 2012).Because of the great increment of DF volume a huge volume of sediment could be produced even if the initiating debris flow is of a small size.This aspect is crucial for a correct DF hazard mapping.Nevertheless, in the present study, attention is drawn to the link between the rainfall characteristics and the possible triggering of DF.The flow propagation and deposition are not simulated.The possible occurrence of a landslide-triggered DF is investigated, aiming at identifying the rainfall amounts that are critical for those events, thus providing a useful tool for DF early warning systems.
The modelling approach presented is based on the simulation of a large number of cases covering the entire range of rainfall intensity, rainfall duration and antecedent rain and considering different combinations between the three of them.The total debris volume, available for the flow, is computed in each simulation.From this, a database is built up in order to obtain rainfall threshold curves.When operating in real time, if the observed and forecasted rainfall exceeds a given threshold, the corresponding probability of DF occurrence can be estimated.Warning for possible DF occurrence may be given congruently with these results.

The shallow landslide triggering model
The process by which shallow landslides are triggered by rainfalls and develop into DF is very complex and some assumptions are done to simulate the process with a simple model.The soil is modelled as a single layer discarding the possible effects of the layered soil stratigraphy.The contribution of the vegetation to soil stability is not taken into account.The stability conditions are evaluated for any single computational cell without taking into account the possible interactions between neighbouring cells.For example the single computational element that is simulated to be unstable, may not move if it is surrounded by stable elements.Not taking into account this mechanism may overestimate the instable volume.The model does not investigate the mechanisms governing the development of a local instability into a debris flow (e.g. the occurrence of liquefaction).The assumptions are made that (1) all the unstable volume may potentially contribute to the generated DF and that (2) if the total volume available for the flow overcomes a given threshold, the resulting DF will affect the downstream areas.
The possible triggering of a debris flow is simulated, in a generic element of the basin, by an infinite slope stability analysis (Taylor, 1948;Iverson, 2000).At any depth from the surface (Z), and at any time (t), the factor of safety (F S ) is computed by the ratio between the resisting Coulomb friction and the driving stresses induced by gravity: where α is the slope degree, Z is the vertical coordinate, positive downward, c is the soil cohesion, ϕ is the angle of internal friction, γ s is the depth averaged soil unit weight, γ w is the unit weight of groundwater and ψ(Z, t) is the underground water pressure head that depends on the vertical coordinate and time (t).When a critical value of F S is reached (e.g.F S = 1) the soil over the Z depth is considered unstable.
Many observed DF events have been triggered by longterm low intensity rainfall followed by a short-term heavy rainfall (Crozier, 1989;Wieczorek and Glade, 2005).As a consequence, the triggering groundwater pressure is calculated by superimposing the effect of an antecedent rainfall and an event rainfall.The groundwater pressure response to antecedent rainfall is used as the initial condition for the time-dependent computation of the groundwater pressure response to the event rainfall.
If the antecedent rainfall has a sufficiently low intensity and long duration, steady-state conditions are reached and the direction of the groundwater flux may be assumed to be parallel to the slope.Under this condition, the groundwater pressure, at the initial condition (t = 0), can be calculated by where d is the water table depth, measured in the Z direction, in steady-state conditions.Following Montgomery and Dietrich (1994), the mass conservation equation of groundwater flow gives the following: where Z T is the depth of the impermeable bed, (I z ) steady the infiltration rate at ground surface in the normal slope direction, in steady conditions, K x the hydraulic conductivity in the parallel slope direction, A the drained catchment, b the width of the slope element along the direction tangent to the local topographic contour.Following the approach of Iverson (2000), the short-term response to rainfall may be assessed assuming vertical infiltration.An analytical solution was proposed of a linearized Richards equation, valid in the assumption of almost saturated initial conditions.The boundary conditions are transient groundwater vertical flux equal to zero at great depths below the water table and water entry at ground surface governed by Darcy's law.In these conditions, the water pressure heads are given by (Iverson, 2000) where T is the duration of the event rainfall, ψ(Z, 0) the groundwater pressure head at the beginning of the event rainfall, I Z the infiltration rate at ground surface in the normal slope direction, K Z the hydraulic conductivity in the normal slope direction and R(T * ) defined as follows: in which: where D 0 is the maximum characteristic diffusivity, governing the transmission of pressure heads when the soil is close to saturation.The diffusivity D 0 is estimated with where C 0 is the change in volumetric water content per unit change in pressure head when the soil is close to saturation.The soil water retention curve is described by means of the analytical equation proposed by van Genuchten (1980) and C 0 is estimated through the derivative of the van Genuchten equation.

Model implementation
In order to assess if a specific basin may give place to the formation of a debris flow, the instability simulation previously described is performed in any computational cell and for any combination of rainfall characteristics.For each simulation, the ratio between the number of unstable cells and the total number of basin cells is computed.Throughout the paper, this value is referred to as failure percentage (FP).As the depth of instable soil is computed for any cell, the total instable volume is also provided as output.
The input variables that feed the model are divided into two main families, static and dynamic.Static variables are the morphological features (A/b, Z T , α) and the soil parameters (c, ϕ, γ s , K x , K Z ,D 0 ).These are considered as stationary at the process scale.The dynamic variables are the rainfall related variables ((I z ) steady , I Z , T ).These are computed, at regular intervals, between user-defined minimum and maximum values.
The effect of the uncertainties of the soil properties are evaluated using a Monte Carlo methodology.Most of the properties (c, ϕ, γ s , K x , K Z , D 0 , Z T ) are considered to be stochastic: instead of defining a value, a range and a PDF is provided.The PDFs express the measurement errors and the uncertainties related to the spatial variability of the soils parameters.The parameters space is sampled using the PDF distribution.Different realizations of the parameters fields are carried out.No cross correlation is considered between different soil properties.Because of the assumptions made, the computation of water pressure at any time and for any cell (Eqs.2, 3 and 4) is not affected by the other cells.As a consequence in the parameter spatial fields, no spatial correlation needs to be considered along the basin, and a random PDF process can be used.For every cell the model randomly selects a properties sample, fitting the user provided PDF.The size of the sample is 10 for every parameter.The geomorphologic properties (A/b, α) are considered deterministic because they depend on the topography, and a low error is assumed in this data.
As an example, the total amount of cases to evaluate for the test bed (Sambuco Basin) is computed: the testing mesh size in the dynamic variables space ((I z ) steady , I Z , T ) is 10 × 50 × 50.The basin DTM (digital terrain model) contains 400 000 cells.As the size of the sample is 10 for each of the seven stochastic parameters, the number of all the combinations of the sampled parameters is 10 7 .To evaluate the stability of the soil 10 different depths are selected along the soil column.As a consequence, the total amount of cases to evaluate for this basin including cells, column depths and soil properties is about 10 18 .The resulting huge computational time prevents practical implementation of the model.To overcome this difficulty, two simplifications are adopted.The soil parameters are considered statistically independent, so not all the combinations of the sampled parameters values are computed.This reduces the number of cases related to the soil properties to 10.Moreover, a sample of the topography is selected randomly.For example if 10 % of the total amount of cell is considered, the number of cases related the topographical grid becomes 40 000.With these simplifications the total amount of cases to be evaluated, for the cited example, become 10 11 .With this amount the model run lasts a few hours on a personal computer with average potentialities.Another important computational reduction comes from the fact that once a cell becomes unstable it is not necessary to compute the rest of the rainfall event durations and intensities.A longer duration will trigger instability as well as a greater intensity.
In order to check the performance of the model, a sensitivity analysis was carried out.For a given rainfall, simulations were performed in which the ratio among the number of computational cells included in the sample and the total amount of basin cells ranges from 0.5 % to 100 %.This ratio is called n c . Figure 1 shows the failure percentages as function of different n c values.As the computational cells are selected randomly, when n c is small, two different runs of the same simulation conditions give different value of FP.However, when n c becomes greater the results converge to the solution obtained with the 100 % of the cells (FP = 2. in example).The simulation results presented in this paper are obtained with n c equal to 20 %.
The simulation results are interpolated to construct rainfall intensity-duration curves for each antecedent rain value (I z ) steady .This gives fixed values of FP or of the total volume of possible DF.

Pressure head response to variation in rainfall duration
After several runs of the model were carried out a particular behaviour was observed, linking duration and total amount of event rainfall.In order to clearly capture this feature a brief theoretical analysis was performed, for a better understanding of the model outcomes.The stability conditions, consequent to different rainfall events having the same accumulated rainfall (and therefore rainfall intensity decreasing with increasing rainfall duration), were explored.After substitution into Eq.( 4) of I z = H /T , where H is the accumulated rainfall, the time derivative of ψ(Z, T ) can be easily obtained for the case of constant H . From the analysis of the sign of the derivative results it appears that the function ψ/Z increases with rainfall duration when the following condition applies: Since this is implicit in the variable T * it was solved numerically.It results that the function ψ/Z increases with rainfall duration when: Substituting Eqs. ( 9) into ( 6), the ratio ψ/Z results to be increasing with rainfall duration when T is lower than a critical value T crit given by Table 1.Soil parameters, slope degree and initial condition used to generate Fig. 2.
Figure 2 shows two examples.The soil parameters and the rainfall characteristics are the same for the two panels, except the D 0 that in panel B is 10 times greater than in panel A. The input data used for these simulations are reported in Table 1.
An important consequence of this behaviour is that, in case of small diffusivities, the water pressures continuously increase at constant accumulated rainfall with increasing rainfall durations until durations are lower than a rather high T crit .For example, in panel A of Fig. 2 T crit = 7.2 h.The expression of T crit given by Eq. ( 10) is similar to the timescale Z 2 T /D 0 that has been indicated by Iverson (2000) as the minimum time necessary for strong normal slope pore pressure transmission from the ground surface to depth Z T .This means that the analytical solution of the Richards equation (Eq.4) proposed by Iverson (2000) can only correctly estimate the effects of rainfalls with duration greater than the critical value.On the other hand, the range of rainfall duration that can be simulated is also limited by another timescale (A/D 0 ) that expresses the minimum time necessary for strong lateral pore pressure transmission.The ranges of possible rainfall durations are reported in Table 2 for different values of D 0 .As a result, for the cited example, when D 0 is lower than 10 −4 m 2 s −1 , only the effects of daily rainfall can be properly assessed by the model.

Study case
The Sambuco Basin is a steep coastal watershed of Amalfi Peninsula, southern Italy (Fig. 3).It covers an area of about 6.4 km 2 .The elevations spread from sea level to 1000 m, with a mean elevation of 422 m a.s.l.The mean slope is 32 • .The main river channel is 4.8 km long with a N-S orientation.The site consists of a set of small and steep catchments covered by a series of pyroclastic deposits dating back to the Somma-Vesuvius volcanic eruptions.
The topographic features (A/b, α) of the basin have been derived through a GIS-based approach over a 4 m ×4 m DTM.

Estimation of the soil properties of the model
The soils covering the carbonate ridge of the basin are made of ashes, pumice and scoriae, coming from the eruptive activity of the Somma-Vesuvius volcano (eruptions immediately before, and after the 79 AD, until the last one in 1944).These deposits cover many of the slopes of the Campanian  Appennine, and particularly of the studied area.Weathering and pedogenetic processes of such deposits have produced soils with peculiar characteristics, referred to as "andic properties" (Basile et al., 2003).One major problem in drawing up a detailed map of the deposits, and therefore a good prevision of possible instable soil volume, is in identifying the most appropriate method for defining the spatial distribution of the different soil typologies.For the Sambuco Catchment a detailed soil map was produced in the framework of previous studies (Papa et al., 2011a).A synthetic description of the procedure is given below.Firstly 12 morphological units were identified and mapped through DEM manipulations.At the same time, the soil profiles were described in more the 20 sites inside the Sambuco Basin.The observed correlation between surface geomorphology and spatial variability of soils and deposit was then used to associate a typical horizons sequence to each of the 12 morphological units.The soil properties (γ s , ϕ, c, K x , K z ) were estimated from literature data (Basile et al., 2003;Bilotta et al., 2005;Ciollaro and Romano, 1995;Iamarino and Terribile, 2008) for any layer of the profile.The soil properties for any morphological unit were then extrapolated from the soil properties of the layers that constitute the corresponding profile (see Table 3).
A set of more than 500 punctual observations of the soil thickness was used to further divide the 12 morphological units into areas with homogeneous soil depth.In this way 21 geomorphological homogeneous units (Fig. 5) were identified.For each unit the soil depth (Z T ) was then assigned.It is worth noting that, in the model, the soil thickness limits the volume of sediments available for the formation of a shallow landslide and therefore limits the maximum possible DF volume.
The diffusivity D 0 was estimated with Eq. ( 7) where C 0 was estimated through the derivative of the van Genuchten equation.Literature values (Basile et al., 2003;Ciollaro and Romano, 1995) of the parameters of the van Genuchten equation for the pyroclastic soils of Campanian Appennine were used.From this it can be shown that for a saturation degree of 90 %, 1/C 0 equals approximately 2. As a consequence the value of the characteristics diffusivity was estimated as The input static variables of the Sambuco Basin are reported in Table 3.
T crit (Eq.10) is calculated for each basin cell as it depends on the different values of soil depth and permeability in the 21 soil districts.T crit ranges from few minutes to few days.As shown in Fig. 7 T crit is less than 5 h for about 80 % of the basin area and is less than 2 h for only 40 % of the basin.As discussed in the above paragraph, the analytic solution of the Richards equation (Eq.4) implemented in the model, can only correctly estimate the effects of rainfall amounts with duration longer than the critical value.As a consequence we can reasonably rely on the simulation of the effects of rainfall longer than 5 h while the results obtained for shorter durations should be regarded with caution.In these cases a complete solution of the Richards equation should be implemented.

Past DF events
The area has been affected by extreme weather events with catastrophic consequences (Ciervo et al., 2012;De Luca et al., 2010;Esposito et al., 2003;Papa et al., 2011a).The pyroclastic soils covering the carbonate rock of the Campanian Appennine are often affected by DF events (Cascini et al., 2008;Martino and Papa, 2008).On 25 October 1954, an extraordinary rainfall event hit the area of the Amalfi Coast and Salerno (≈ 80 km 2 ); severe flooding and landslides caused 318 fatalities and large-scale damages (Frosini, 1955).The Sambuco Basin was on the west margin of the affected area.Slope failures mainly occurred on the east side of the basin.The reconstruction of the meteorological event made by the Servizio Idrografico e Mareografico Italiano (SIMI, 1921(SIMI, -1996) ) has drawn the map of isohyets.The average daily cumulative rainfall in the Sambuco Basin was estimated to be 321 mm and the maximum hourly intensity 86.1 mm h −1 .Papa et al. (2011a) mapped the tracks of the detachments that occurred in 1954, first geo-referencing the aerial photographs of that year and then checking and correcting the resulting map with field observations.It resulted in 2.8 % of the total basin area being affected by detachments (Fig. 4).
The volume involved by landslides was derived with GIS tools, overlaying the map of the detached area with the map of soil thickness.The estimate closes 300 000 m 3 .
This unstable volume formed a DF that flooded the downstream village of Minori with an estimated peak discharge of about of 58 m 3 s −1 (Papa et al., 2011a).
Another event was observed in 2005.Only 0.3 % of the basin area was mobilized and the generated DF did not reach Minori.
The events occurred in Sambuco Catchment, as well as many of the landslides of the flow type of Campanian Appennine, are initially shallow non-confined slides (or flow) along steep slopes and evolve into in-channel flows as they travel downstream.
Figure 6 gives the granulometric curves of the flow mixture.It may be noted that the content of silt and finer particles is about 20 %.Rheological tests (Martino and Papa, 2008;Martino, 2003) were carried out on samples taken in areas affected by the DF phenomena in May 1998 in Sarno.The mixtures involved in these events have the same main characteristics and genesis as the Amalfi Coast soils.As it is shown in Fig. 6, the granulometry of the samples taken in Sarno is analogous to the ones taken in the Sambuco Catchment.As a consequence, the rheological characterization of the Sarno DF can be used to estimate the properties of the DF mixtures of the Amalfi Coast basins, such as the Sambuco Basin.The rheological tests showed a shear thickening behaviour, described by a Herschel Bulkley equation (Coussot, 1997) with the exponent of strain rate equal to 1.7.This means that the peculiar composition of the pyroclastic soils of Campania Region is responsible for a rheological behaviour that is influenced by the presence of both the fine and the coarse particles.The fine fraction gives rise to a threshold-type cohesive effort and a high viscosity.The presence of coarser sediments gives rise to the collisional behaviour.3; nn = absence of soil cover (e.g.urbanized areas, rocks, streets).(Papa et al., 2011a).The values of the static variables for each of the district (numbered from 1 to 21) are reported in Table 3; nn = absence of soil cover (e.g.urbanized areas, rocks, streets).
With reference to plasticity properties, the samples taken in Sarno catchment are classified as non-plastic or slightly plastic soils (Bilotta et al., 2005).
With regard to the classification proposed by Hungr et al. (2001), the phenomenon under study can be classified as DF, considering the low fine content and plasticity index of the mixture, even if the rheological tests show a combination of the resistance mechanisms of a DF together with the ones of a mud flow.

Critical rainfall thresholds derived by the model
The intensity-duration curves (ID curves) for any fixed value of antecedent rain were derived from the data sets that were obtained with the numerical hydrological, slope stability model described above.Any ID curve corresponds to a fixed value of the failure percentage (FP).Figure 8 shows simulation results (sim01) assuming null antecedent rainfalls.With increasing antecedent rainfall, the ID curves move downward (Papa et al., 2011b).The possibility of assigning a PDF to the soil variables, as described in the above paragraph, is not used in sim01 where the soil parameters are assigned in a deterministic way.
The simulation results are compared to the ID curve relative to the event that occurred in October 1954 (Fig. 8).The months before this event were dry and therefore the comparison was carried out with the results obtained for antecedent rain equal to zero.The input static variables are fixed as explained in the above paragraph, except for the soil cohesion that was incremented by 25 % after calibration.This is the only calibration parameter of the model.Results show that a higher value of the soil cohesion is needed in order to reproduce the observed failure percentage (FP).Possible explanations of this higher value may be the stabilizing effects of the plant roots and of the possible stable elements surrounding the unstable ones.As stated above, both these contributions to soil stability are neglected by the model.
For an event rainfall with duration of about 6 h, the ID curve of the 1954 event approaches the ID curve corresponding to a FP of 3 %.This result is in clear agreement with the observed FP that, as assessed above, was about 2.8.The rainfall duration that is critical for DF formation is 6 h; this result is consistent with the timing of the DF.
It is worth noting that, as discussed in the previous paragraph and due to the assumption of vertical flow (Iverson, 2000), only rainfalls with durations greater than 5 h are correctly simulated by the model in the studied basin.
The 2005 event was also studied.This event followed a moderate rainy period; the total rainfall amount of the month before was 212 mm.The Montgomery and Dietrich (1994) model is used to compute the initial water table depth at the beginning of the event rainfall, depending on the antecedent rainfall.As the event of 1954 was in October, after the dry period, the antecedent rainfall was null in that case.For the event of 2005 the stability condition consequent to a stationary rainfall of 212 mm month −1 was computed with the Montgomery and Dietrich model and a FP of about 1.2 % was obtained.This value is higher than the observed one (about 0.3 %).This result confirms that, as reported in literature (e.g.Godt et al., 2008), the hypothesis of steady groundwater flow conditions is likely to overestimate the failure percentage.
If the antecedent rainfall is neglected, the input of the simulation (and therefore the result) is the same as Fig. 8.We can observe that the event rainfall, with duration of about 8 h, approaches the simulated ID line corresponding to FP of 0.3 %.This is equal to the observed FP of the event.For this event, the modelling results suggest that the effects of the antecedent rainfall are barely noticeable.
The effect of the uncertainties of the soil parameters is investigated.Figure 9 shows the results of three different simulations.Sim01 is the one discussed above in which the soil parameters were assigned with a spatial distribution and in a deterministic way.For Sim02 the soil parameters of the entire catchment are taken from one district: district number 13 of Table 1.This is the district with the higher soil depth in the "Z.O.B. and hollow" morphotype.Sim02 aims at reproducing the frequent situation when only little field data are available and the soil properties of a basin are roughly estimated from few point measures taken in a small area.In Sim02 the soils parameters are assigned deterministically as in Sim01.Sim03 is performed assigning an average value to the soil parameters equal to the fixed value of Sim02, but with a confidence interval variable depending on the parameter uncertainties.The confidence intervals were set, for γ s and ϕ equal to the 10 % of the average while for the Z T , c, K x and K z equal to the 20 % of the average.Figure 9 gives the CRT corresponding to FP of 1 % and 3 %, obtained with sim01, sim02 and sim03.For FP = 3 % the CRT obtained with sim02 is greatly higher than the one obtained with sim01.This is probably due to the fact that we assigned a set of homogeneous soil characteristics to the entire basin that is less likely to produce instabilities than the spatially distributed ones of Sim01.Therefore in order to have the same FP a rainfall with higher intensity or longer duration is needed.The results of sim03 show that if the soil parameters are assigned in a stochastic way the solution is more similar to the one of sim01.This means that the lack of knowledge about the spatial distribution of the parameters underestimates the hazard.The stochastic treatment of the soil parameters can partially mitigate this problem.In the reported example, for the lower FP (FP = 1 %) the difference among sim01 and sim02 is less accentuated and the CRT of sim03 is lower than the one of sim01, thus giving an overestimation of the hazard.
The sensitivity of the model to changes in the diffusivity (D 0 ) was investigated too.The same input data of the simulation discussed above (sim01, Fig. 8) were used in a new simulation (sim04).The only difference between the two simulations is the value assigned to D 0 : D 0 = 2K z for sim01 and D 0 = 0.1K z for sim04.T crit change with D 0 ; in case of sim01 T crit ranges from a few minutes to a few days and 80 % of the basin has T crit less than 5 h.In the case of sim04 T crit ranges from one hour to one month and only 15 % of the basin has T crit less than 5 h while 80 % of the basin has T crit longer than 1 day.Figure 11 compares two simulations.In this case, the accumulated rainfall is reported in the graph instead of the rainfall intensity.Each curve of accumulated rainfall versus rainfall duration corresponds to a fixed value of the generated failure percentage.As expected, when D 0 is smaller, the basin is more stable.In the case of sim01 the considered rainfall durations (from 0 to 24 h) are smaller than T crit for the major part of the basin and the critical accumulated rainfall curves increase with rainfall duration.On the contrary, in the case of sim04, T < T crit for almost the entire basin, and the critical accumulated rainfall curves decrease with rainfall duration.This result is consistent with the observation made in the above Sect.2.2.Attention should be paid on the values of T crit as the simulations of rainfall durations with T < T crit give inconsistent results.

Discussion
Simulated ID curves, corresponding to given FP, may be regarded as CRT.If rainfall overcomes these curves the corresponding FP is expected to occur.Once the CRT curve is fixed the number of false alarms in a playback period can be evaluated.Figure 8   The simulation results have been compared with a rainfall threshold line, derived through the elaboration of empirical data relative to the pyroclastic deposits of Campania Region (Calcaterra et al., 2000).From the comparison with the simulation results, the threshold line proposed by Calcaterra et al. (2000) corresponds to a failure percentage of about 0.2 % (Fig. 8).
Intensity duration lines obtained by the simulations may be directly used as a rule for providing DF warnings, once a threshold is fixed for the failure percentage.
Such a threshold may be fixed by taking into account that when the unstable areas are not wide enough, the mobilized soil is not able to reach the downstream vulnerable areas.When a large number of observed slope instabilities is available, the threshold may be fixed by searching for which of the consequent DF reached the downstream vulnerable areas.In the studied example, we have two observations: the 1954 event (FP = 3 %) that caused great damage to the Minori village and the 2005 event (FP = 0.3 %) that did not reach the downstream village of Minori.Consequently the critical FP should be fixed between 0.3 % and 3 %.
When the historical database is not large enough, or in the case of total absence of historical DF events, the threshold FP can be fixed by simulating the downstream effect of DF having different total volumes.These simulations can be performed through the mathematical and numerical modelling of DF propagation (e.g.Fraccarollo and Papa, 2000;Iverson and Denlinger, 2001;O'Brien et al., 1993;Malet et al., 2005).Due to the entrainment process that can occur in the propagation phase, an initial small volume can develop into DF that eventually involve large volumes.This results into a significant increase of DF hazard.Models that simulate the DF propagation and DF amplification are available in literature (e.g.Egashira et al., 2003;Medina et al., 2008;McDougall and Hungr, 2005;Quan Luna et al., 2012).For any input volume the DF simulation models' output is a map of DF intensity (flow depth and velocity).This map is overlapped with the vulnerability map in order to derive the DF risk map.The threshold of initial DF volume is then chosen depending on the risk level that is politically acceptable.One important advantage of the proposed methodology is that a set of CRT can be identified; each CRT corresponds to a fixed level of DF risk.
Once the volume threshold is fixed, a graph similar to the one shown in Fig. 10 may be used as a rule for DF warnings.In this kind of graph, the simulation results are elaborated in order to show, for any antecedent rain, the ID curve, giving place to a fixed value of the total amount of available debris volume.The ID curve of 1954 event lays between the ID curves corresponding to a total unstable volume of 300 000 and 400 000 thus being in clear agreement with the estimated total volume of the event.

Conclusions
A methodology and a numerical tool have been developed for the computation of rainfall critical thresholds for debris flow early warning systems.The method is suitable to assess the possible occurrence of debris flows and mud flows triggered by shallow landslides or debris avalanches.A simple model has been implemented based on the analytical solution of the Richards equation valid in the hypothesis of vertical infiltration (Iverson, 2000).The model performs stability simulations for any possible combination of rainfall duration and intensity.These simulations are repeated varying the antecedent rainfall whose effects are computed in the hypothesis of steady-state conditions (Montgomery and Dietrich, 1994).
The range of applicability of the analytical solution of the Richards equation proposed by Iverson (2000) has been investigated.It resulted that if the rainfall duration is lower than a critical value the pressure heads unrealistically increase with rainfall duration at constant accumulated rainfall.The threshold value of the rainfall duration for the correct application of the model has been analytically derived.This threshold should be regarded in any application of the Iverson (2000) model.
The developed tool computes intensity duration curves corresponding to a fixed total volume of unstable basin's elements.The assumption is made that the initial shallow landslides have sufficient energy to evolve into a proper flow and inundate the downstream villages only if the total unstable volume is sufficiently high.On this hypothesis the simulated intensity duration curves can be regarded as debris flow critical rainfall thresholds.Given the simplifications usually introduced in slope stability models and the uncertainties in the estimation of the soil parameters, the precise location of the unstable elements along the basin is not easy to assess.For this reason the lumped approach here presented, in which only the sum of the volumes of the unstable elements is estimated, may provide better results.
The CRT derived with the developed tool have been compared to rainfall measurements of past DF events in a studied catchment in the south of Italy (Amalfi Coast).
The comparisons showed a clear agreement for the simulations in which the antecedent rainfall is absent or neglected.It results that the implemented model overestimates the effects of the antecedent rainfall.A more complex model, based on the solution of the complete Richards equation could overcome this problem together with the limitations in the duration of the rainfalls that can be simulated.On the other hand, the computational time would greatly increase and it would be necessary to estimate a greater number of soil parameters.
The application of the model to other study cases with a greater number of observed past events is necessary to better understand the capabilities of the proposed approach.
A preliminary analysis is presented about the effects of uncertainties in the soil parameters.It appears that assigning a probability distribution function to the most uncertain parameters can partially mitigate the negative effects of unknown distribution of soil properties.Further analyses are needed in this direction.
In practical application, the developed tool can provide the link between rainfall characteristics and debris flow initial volume.Afterwards, a simulation of the propagation, amplification and deposition of possible debris flows, having different initial volumes, can associate the debris flow initial volume with a map of areas at risk.One important advantage of the proposed methodology is that a set of critical rainfall thresholds can be identified; each one corresponding to a fixed level of expected debris-flow risk.The derived critical rainfall thresholds can be used to implement early warning systems also in basins where empirically derived threshold cannot be derived because the database of past events' observation is poor or the stability conditions have changed for natural or anthropic reasons.

Figure 1 .Fig. 1 .
Figure 1.Failure percentages (FP), computed by different runs of the model, depending on the ratio (nc) among the number of computational cells included in the sample and the total amount of basin cells.

Figure 5 .
Figure 5. Map of the geomorphological homogeneous districts (Papa et al. 2011a).The values of the static variables for each of the district (numbered from 1 to 21) are reported in Table3;

Fig. 5 .
Fig. 5. Map of the geomorphological homogeneous districts(Papa et al., 2011a).The values of the static variables for each of the district (numbered from 1 to 21) are reported in Table3; nn = absence of soil cover (e.g.urbanized areas, rocks, streets).

35Figure 7 .Fig. 7 .Figure 8 .Fig. 8 .
Figure 7. Threshold values of the rainfall duration for the correct application of the model Fig. 7. Threshold values of the rainfall duration for the correct application of the model.

Figure 10 .Fig. 10 .
Figure 10.Comparison between registered rainfalls and simulated ID curves at fixed values of total instable volumes (thousands of m 3 )

Table 3 .
Values of the static variables for the 21 homogeneous districts of Sambuco Basin.