Evaluation of shallow landslide-triggering scenarios through a physically based approach: An example of application in the southern Messina area (northeastern Sicily, Italy)

Abstract. Rainfall-induced shallow landslides are a widespread phenomenon that frequently causes substantial damage to property, as well as numerous casualties. In recent~years a wide range of physically based models have been developed to analyze the triggering process of these events. Specifically, in this paper we propose an approach for the evaluation of different shallow landslide-triggering scenarios by means of the TRIGRS (transient rainfall infiltration and grid-based slope stability) numerical model. For the validation of the model, a back analysis of the landslide event that occurred in the study area (located SW of Messina, northeastern Sicily, Italy) on 1 October 2009 was performed, by using different methods and techniques for the definition of the input parameters. After evaluating the reliability of the model through comparison with the 2009 landslide inventory, different triggering scenarios were defined using rainfall values derived from the rainfall probability curves, reconstructed on the basis of daily and hourly historical rainfall data. The results emphasize how these phenomena are likely to occur in the area, given that even short-duration (1–3 h) rainfall events with a relatively low return period (e.g., 10–20~years) can trigger numerous slope failures. Furthermore, for the same rainfall amount, the daily simulations underestimate the instability conditions. The high susceptibility of this area to shallow landslides is testified by the high number of landslide/flood events that have occurred in the past and are summarized in this paper by means of archival research. Considering the main features of the proposed approach, the authors suggest that this methodology could be applied to different areas, even for the development of landslide early warning systems.


Introduction
Landslides triggered by rainstorms occur in many parts of the world and cause significant damage and loss to affected people, organizations and institutions as well as to the environment (Glade, 1997;Nadim et al., 2006;Petley, 2012). Within this category of natural disasters, shallow landslides (in particular debris flows)can pose a serious threat to life or property, in particular due to their high velocity, impact forces and long runout, combined with poor temporal predictability (Jacob and Hungr, 2005). These phenomena consist of sudden mass movements of a mixture of water and granular material that rapidly develop downslope, eroding the soil cover and increasing their original volume (Iovine et al., 2003). Due to their high destructiveness, the study of these processes is an important research topic that can support decision makers in developing more detailed land-use maps and landslide hazard mitigation plans.
However, it is not simple to predict the probability of the occurrence and magnitude of shallow landslides, considering the complexity of the phenomenon, mostly related to the variability of controlling factors (e.g., geology, topography, climate and hydraulic conditions). In this respect, a relationship between triggering events (i.e., rainfall) and landslide occurrences needs to be established. To evaluate this causeeffect relationship, an approach widely used in the literature relies on the definition of empirical thresholds (Caine, 1980;Reichenbach et al., 1998;Wieczorek and Glade, 2005;Guzzetti et al., 2007). An empirical threshold defines the rainfall, soil moisture or hydrological conditions that, when reached or exceeded, are likely to trigger landslides (Reichenbach et al., 1998). Rainfall intensity-duration thresh-L. Schilirò et al.: Evaluation of shallow landslide-triggering scenarios olds for the possible occurrence of landslides are defined through the statistical analysis of past rainfall events that have resulted in slope failures, and can be classified on the basis of the geographical extent for which they are determined (i.e., global, national, regional or local thresholds) and the type of rainfall information used to establish the threshold (Brunetti et al., 2010).
Nonetheless, the reliability of empirical thresholds is generally affected by quality and availability of the historical data. In fact, adequate historical data on landslides and simultaneous rainfall are in most cases available only for a relatively short period, which may not be sufficiently significant from a statistical point of view. Furthermore, rainfall intensity and duration alone may not be able to capture most of the uncertainty related to landslide triggers (Peres and Cancelliere, 2014). Another drawback of the empirical rainfall thresholds is a general lack of spatial resolution. This aspect cannot be neglected if we consider that the terrain factors which control the onset of instability during a rainfall event can vary spatially to such an extent that, from a theoretical point of view, the rainfall threshold can be different for each landslide (Lo et al., 2012). Finally, further criticisms are based on the observation that it is not the amount of precipitation but the (largely unknown) amount of water that infiltrates and moves into the ground, and causes slope failure (Guzzetti et al., 2008).
For this reason, in recent years different models have been developed to define physical (process-based) thresholds. Specifically, these models can determine the amount of precipitation needed to trigger slope failures, as well as the location and time of expected landslides, using spatially variable characteristics (e.g., slope gradient, soil depth and shear resistance) with a simplified dynamic hydrological model that predicts the pore-pressure response to rainfall infiltration (Montgomery and Dietrich, 1994;Wilson and Wieczorek, 1995;Terlien, 1998;Frattini et al., 2009). These models, although they are challenging to apply over large areas where a detailed knowledge of input parameters is very difficult to acquire (Berti et al., 2012), are usually calibrated using rainfall events for which rainfall measurements and the location and time of slope failures are known.
In this paper we propose an approach based on TRIGRS (Baum et al., 2002), a physically based model that predicts the timing and distribution of shallow, rainfall-induced landslides combining an infinite slope stability calculation with a transient, one-dimensional analytic solution for porepressure response to rainfall infiltration. This model has been used in order to define different shallow landslide-triggering scenarios in the study area (located SW of Messina in northeastern Sicily, Italy) by varying the rainfall input on the basis of the results deriving from the analysis of the historical rainfall data. Prior to this stage, the model was thoroughly validated through the back analysis of the disaster which occurred in the same area on 1 October 2009. On that day, a heavy rainstorm triggered several hundreds of shallow land-slides, causing 37 fatalities and severe damage to buildings and infrastructure (Schilirò and Esposito, 2013). Given the nature of the event, it can be considered particularly representative of the studied phenomenon and, thus, suitable for testing the reliability of the physically based model. The paper is organized as follows: after a brief description of the study area and the 1 October 2009 event (Sect. 2), a summary of the landslide/flood events that have occurred in the past is reported (Sect. 3). Then, the methods used for the parameterization of the TRIGRS model and the analysis of the historical rainfall data are outlined (Sect. 4). Afterwards, the results of the back analysis of the 2009 event and the evaluation of possible future triggering scenarios are provided (Sect. 5) and discussed, along with the main features of the proposed approach (Sect. 6). Finally, in Sect. 7 the main conclusions are summarized.

General features of the study area and the 1 October 2009 event
The study area ( Fig. 1) is located south of Messina (northeastern Sicily, Italy) at the NE termination of the Peloritani Mountain Belt, that represents the southern border of the Calabrian-Peloritan arc. This chain is composed of different metamorphic units (Kabilo Calabride Complex) of prealpine age and later involved in Hercynian and Alpine orogenic processes, and involved in tectonically overlapping the sedimentary Maghrebidian units (Lentini et al., 2000). Since the Late Miocene, the opening of the Tyrrhenian Basin led to the formation of an extensional fault system that involved and re-oriented some of the former structures. These faults, generally oriented NE-SW, have influenced the development of this region during the Pleistocene-Holocene (Antonioli et al., 2003;Di Stefano et al., 2012), resulting in a landscape characterized by steep slopes eroded by torrent-like straight watercourses, with alluvial conoids and debris flow fans along the valleys. A thin (0.5-2 m) layer of colluvial deposits or coarse-grained regolith overlays the majority of the slopes, where small outcrops of marine terraces, documenting the different uplift rates during the Late Pleistocene, can be found. Three orders of terraces can be distinguished at approximately 185, 135 and 95 m a.s.l. (Catalano and De Guidi, 2003), whereas there is no evidence of fluvial terraces. Catchments generally have small dimensions and markedly elongated shapes, with a short time of concentration and direct discharge into the sea. In particular the study area, which has an extension of about 8 km 2 and a maximum height of about 700 m a.s.l., comprises three main catchments (Giampilieri, Divieto and Racinazzo torrents), highly affected by the heavy rainstorm which occurred on 1 October 2009. This event was characterized by an extremely high rainfall amount in a few hours, for instance the Santo Stefano di Briga monitoring station, which is one of the rain gauges closest to the study area ( Fig. 2a), recorded very high rainfall peaks (e.g., 18.5 mm  Istituto Geofisico), whose locations are shown in the sketch on the upper left (the red square represents the study area). The shaded areas indicate the landslide timing during the event, on the basis of witness reports: 1 -first landslide events in Giampilieri (17:00-17:15 UTC); 2 -beginning of the most critical stage, with the triggering of hundreds of shallow landslides in Giampilieri and the surrounding areas (18:00 UTC); 3 -landslides in Molino (18:30-18:45 UTC); 4 -end of the landslide events (20:00 UTC); (b) accumulated rainfall between 18:00 and 21:00 UTC, based on radar (satellite) data (from Melfi et al., 2012).
of cumulated rain between 17:00 and 17:10 UTC) for a total of 225 mm of rain falling in just 7 h (i.e., between 14:00 and 21:00 UTC). The analysis of satellite data (Fig. 2b) highlights the marked localization of the weather system, also confirmed by the low rainfall values recorded in two rain gauge stations approximately 20 km from the study area (Antillo and Messina Istituto Geofisico monitoring stations, see Fig. 2a). Furthermore, it is worth noting that in the days preceding the event, the area was affected by two intense rainfall events: one on 16 September and one in the night between 23 September and 24 September. The cumulative rain in this period was approximately 300 mm and thus the total rainfall from 15 September to 1 October amounted to about 500 mm, which is 80 mm higher than the average October-December precipitation (421 mm), calculated from 79 years of historical pluviometric data recorded at Santo Stefano di Briga monitoring station. These data are directly available on the Sicily Region website (http://www.osservatorioacque.it). As a consequence of such an extreme event, numerous shallow landslides were triggered; according to the classification by Hungr et al. (2014), such landslides can be classified as debris flows (Fig. 3a) and debris slides (Fig. 3b) that frequently evolved to debris avalanches ( Fig. 3c) (Trigila et al., 2015). According to Schilirò et al. (2015), between 15:00 and 17:00 UTC the critical conditions rapidly developed due to a large increase in precipitation. Afterwards, first landslide events occurred: for instance, a witness asserts that the main debris flow in Giampilieri occurred between 17:00 and 17:15 UTC. Finally, after a further rainfall peak (at approximately 18:00 UTC), many other shallow landslides were triggered in Giampilieri and the surrounding areas. In particular, a devastating debris flow suddenly rushed down the Racinazzo valley, crushing buildings, infrastructure and killing 14 people in Scaletta Marina. On the other hand, similar events were recorded slightly later in Molino, between approximately 18:30 and 18:45 UTC, due to the storm motion towards the inner areas. Landslide occurrences were reported until 20:00 UTC, after which no remarkable event occurred in the area. However, the experiences reported by witnesses, along with the damage to buildings, both indicate very fastmoving debris flows.

Previous flooding and related events
The regional climatic setting of northeastern Sicily is typical Mediterranean, whereby rainfall is concentrated during the wet season (October-April), which is when extreme rainfall events generally occur. However, the local orographic control, coupled with marine effects, can highly influence the occurrence and magnitude of such events. Furthermore, the particular drainage network of the area, characterized by loworder streams with high gradients and short lengths, increases the energy of runoff waters during rainfall events, favoring the erosion processes and the transport of the loosened material, even of large dimensions. For these reasons, this part of Sicily has been affected by recurring flood/landslide events in the past. According to the results of archival research, based on the review of technical reports from local authorities, newspapers, local church archives etc., about 46 landslide/flood events have occurred since the 17th century, most of which occurred during the autumn-winter period (Fig. 4). In this respect, it is worth mentioning the "quadruple rainstorm" that affected the whole Messina area on 13 November 1855, during which more than 100 people lost their lives due to the triggering of countless landslides and widespread flooding (Cuppari, 1856). Extreme rainfall events have also affected this region in recent years; in particular, an event similar to the 1 October 2009 one occurred just 2 years before, on 25 October 2007. On this day, approximately 134 mm of rain fell in the area, and this event also featured extremely high intensity peaks (i.e., 29.1 mm in 10 min). The damage to buildings and infrastructure was remarkable; however, unlike the 1 October event, there were no casualties (Aronica et al., 2008). Heavy rainstorms also hit the area after the 1 October 2009 event, i.e., in the night between 9 and 10 March 2010 and on 1 March 2011. During these events, several landslides were triggered in Scaletta Marina and Giampilieri, convincing the authorities to declare a state of natural disaster for the villages of the southern Messina area. However, even though the increase of recorded events during the last 20 years also depends on the increasing number of sources of information, it is important to emphasize that, in the same period, the landslide risk exposure of the area increased substantially due to the enlargement of the urban area as a consequence of poor land-use planning (Del Ventisette et al., 2012). Finally, recent studies (Bonaccorso et al., 2005;Arnone et al., 2013) indicate an increasing trend of extreme, short-duration rainfall events over the last few decades in Sicily, especially in coastal areas.

Methodology
As mentioned above, the approach proposed in this paper consists of two stages: 1. back analysis of the 1 October 2009 event 2. evaluation of different triggering scenarios.
With regards to the first step (Fig. 5), we compared the safety factor (SF) map obtained by using TRIGRS, a wellknown regional, physically based stability model (Sect. 4.1), with the inventory map of the landslides triggered during the event. In this respect, it is important to note that for each of more than 700 mapped landslides, identified through the analysis of high-resolution aerial orthophotos, integrated by field surveys in the days after the event, the landslide deposit has been distinguished from the source area (Trigila et al.,  2015). Therefore, in order to achieve more accurate assessment, only the source areas have been used for the comparison with the numerical simulations. The input parameters of the numerical model (i.e., soil thickness, spatial rainfall rate, hydraulic conditions and geotechnical properties of the colluvial deposits) have been evaluated by means of lab tests, em-pirical methods and numerical simulations (Sect. 4.2). In the next step, we used specific rainfall inputs for the reconstruction of different triggering scenarios, whose return periods (RP) have been defined through a statistical analysis of the historical rainfall data available for the study area (Sect. 4.3).

Theoretical basis of the TRIGRS model
TRIGRS (transient rainfall infiltration and grid-based slope stability model) is a Fortran program designed for modeling the timing and distribution of shallow, rainfall-induced landslides. It combines a transient, one-dimensional analytic solution for pore-pressure response to rainfall infiltration with an infinite slope stability calculation. In the original version (Baum et al., 2002), the infiltration model was based on Iverson's (2000) linearized solution of Richards' equation, with implementation of complex storm histories, an impermeable basal boundary at finite depth and a simple runoff routing scheme (Savage et al., 2003;Salciarini et al., 2006). Introducing a time-varying rainfall input on the ground surface I Z (t), the pressure head response (Z,t) can be computed using the following input parameters (locally variable throughout the model): slope, soil layer depth d lb , depth of the initial steady-state water table d wt , long-term (steadystate) surface flux I ZLT and saturated hydraulic conductivity K s . However, this solution is appropriate for initial conditions where the hillslope is tension-saturated (Fig. 6a). In

L. Schilirò et al.: Evaluation of shallow landslide-triggering scenarios
the second version (Baum et al., 2008), the TRIGRS model was expanded to address infiltration into a partially unsaturated surface layer above the water table by using an analytical solution of the Richards' equation for vertical infiltration (Fig. 6b). TRIGRS uses four hydrodynamic parameters to linearize the Richards' equation through the unsaturated zone, according to the hydraulic model proposed by Gardner (1958): the saturated (θ s ) and residual (θ r ) water content, the saturated hydraulic conductivity (K s ) and a specific parameter linked to the pore size distribution (α G ). If the amount of infiltrating water reaching the water table exceeds the maximum amount that can be drained by gravity, TRIGRS simulates the water-table rise, comparing the exceeding water quantity to the available pore space directly above the water table or capillary fringe, and then, for each time step, applies the water weight at the initial top of the saturated zone to compute the new pressure head (Baum et al., 2010). For the calculation of the safety factor in the unsaturated configuration, the pressure head is multiplied by the effective stress parameter, as suggested by Vanapalli and Fredlund (2000): where θ is the soil water content. This approximation has application to a generalized effective stress law and represents a simplified form of the suction-stress characteristic curve (Lu and Godt, 2008;Lu et al., 2010). Considering the presence of a thin layer of colluvial deposits over the metamorphic bedrock, in this work we assume a finite depth for the impermeable basal boundary and initially unsaturated conditions for the regolith.

Parameterization of the numerical model
In order to define the input parameters of the TRIGRS model, different methods and techniques were used. To estimate the spatial variation of soil thickness, the model proposed by Saulnier et al. (1997) was applied to the study area (Fig. 7a). This model correlates soil depth to the local slope angle according to the following equation: where h i is the soil thickness computed at pixel i, h max and h min are the maximum and minimum soil thickness values measured in the area, α i is the slope value at pixel i, while α max and α min are the maximum and minimum slope values encountered in the study area. The maximum and minimum values of slope and soil thickness, which were measured within the source areas of the shallow landslides triggered during the 1 October 2009 event, are equal to 58-17 • and 1.5-0.5 m respectively, and they can be considered reliable since the 2009 landslides mostly involved the entire soil profile. Although the model proposed by Saulnier et al. (1997)  relies heavily on geomorphological simplifications, it is frequently used to estimate a spatially distributed soil depth field in basin scale modeling (e.g., Salciarini et al., 2006). To reproduce the spatial rainfall distribution of the 1 October 2009 rainstorm, the conditional merging technique (Ehret, 2002;Pegram, 2003) has been chosen as the interpolation method. In this approach, the information from the satellite radar is used to condition the spatial rainfall field obtained by the interpolation of rain gauge measurements. Although there are numerous deterministic methods for estimating spatial rainfall distribution (e.g., Thiessen polygon, inverse distance weighting, polynomial interpolation), geostatistical methods are commonly preferred because they allow not only for spatial correlation between neighboring observations to estimate values at ungauged locations to be accounted for, but also more densely sampled secondary attributes (i.e., weather radar data) to be included with sparsely sampled measurements of the primary attribute (i.e., rainfall) to improve rainfall estimation (Mair and Fares, 2011). In particular, meteorological satellite radars give a large-scale vision of precipitation fields compared to scattered point estimates from rainfall gauges. In this study, the precipitation rate maps, deriving from the processing of EUMET-SAT (European Organisation for the Exploitation of Meteorological Satellites) satellite data, were used. These maps, that were made available by the National Center of Aeronautical Meteorology and Climatology (CNMCA) of the Italian Air Force, are generated from blending PMW (passive microwave) measurements and IR (infrared) brightness temperatures, coupled with the NEFODINA (DYNAmic NEFOanalysis) software, that allows the automatic detection and classification of convective cloud systems, reducing the underestimation of precipitation (Mugnai et al., 2013). Tenminute rainfall records of six stations (Antillo, Colle San Rizzo, Fiumedinisi, Ganzirri, Messina Istituto Geofisico and Santo Stefano di Briga) have been used as input data, after conveniently converting them into 15 min data for the comparison with the corresponding radar rainfall maps. Thus, sequential rainfall maps (Fig. 7b) have been obtained for the time period of between 13:00 and 21:00 UTC.
The hydraulic properties of the colluvial deposits, the steady-state water-table depth and the initial soil moisture conditions have been estimated using the HYDRUS 1-D model (Šimůnek et al., 1998), a USDA (United States Department of Agriculture) Salinity Laboratory software package which can simulate the water flow into unsaturated porous media resulting from a rainfall event. The software describes infiltration in the vadose zone using a modified version of Richards' equation. In this paper, numerical simulations have been performed for the period 1-30 September 2009 in order to quantify the effect of the 1-month antecedent rainfall on soil moisture conditions. The van Genuchten-Mualem model (van Genuchten, 1980) was chosen as a hy-draulic model to simulate water flow, whereas the hydrodynamic parameters θ s , θ r , α G and K s are predicted from soil grain size distribution using the ROSETTA Lite module (Schaap et al., 2001). This module uses a database of measured water retention and other properties for a wide variety of media. For a given grain size distribution and other soil properties, the model estimates a retention curve (i.e., the relationship between soil water suction and the amount of water remaining in the soil θ) with good statistical comparability to known retention curves of other media with similar physical properties (Nimmo, 2005). Daily rainfall data have been used as input for the model, whereas evapotranspiration is accounted for by inserting the maximum and minimum temperature values recorded during the investigated period into the Hargreaves equation (Jensen et al., 1997). As a lower boundary, a zero-flux condition was assumed due to the presence of impermeable bedrock below the soil cover. A 80 cm soil profile inclined to 38 • (i.e., the average soil thickness and slope observed within the landslide source areas) was chosen as the most representative geometric configuration of the slope prior to the 1 October event.
Laboratory tests have been performed to measure physical and mechanical properties of the colluvial deposits (Table 1). The grain size distribution analyses show a soil composed mainly of coarse-grained particles (gravel and sand) with minor components of silt and clay. With regard to the mechanical parameters, drained triaxial tests have been conducted on three large reconstituted specimens (H = 200 mm, D = 100 mm). To reconstitute each specimen, the soil was compacted inside a mould in different layers of decreasing depth, in order to account for under-compaction. The tested material was sieved, leaving a maximum grain size of 10 mm, imposing 35 % of porosity (i.e., the average porosity obtained from different soil samples) and 8 % of initial water content. The latter value can be considered representative of the inves-  (Aronica et al., 2012;Peres and Cancelliere, 2014;Penna et al., 2014) reported values which vary between 30 • and 40 • for the friction angle and between 0 and 5 kN m −2 for the cohesion; thus the resulting internal friction angle (i.e., 36.3 • ), obtained by assuming a null cohesion, substantially agrees with these values. However, it is important to stress that this difference can depend on both the natural spatial variability of soil shear strength parameters and the type of deposit, characterized by an extremely variable texture resulting from erosion and weathering areas.

Analysis of historical rainfall data
In order to depict different shallow landslide-triggering scenarios in the study area, firstly it is necessary to evaluate the recurrence of specific rainfall events, which can be used as input for the physically based model. Therefore, a statistical analysis of historical rainfall data has been performed.
The hydrological-statistical model is based on the analysis of the maximum values assumed by the chosen hydrological variable (i.e., cumulative rainfall at different time intervals). Depending on the type of considered rainfall event (i.e., prolonged or intense short-duration rainfall), the hydrological variables are identified and then the recurrence of the event can be expressed in terms of return period (RP). In this study the probability model relied on the generalized extreme value (GEV) distribution introduced by Jenkinson (1955). This distribution is a generalized version of the more well-known Gumbel distribution, which is largely used in the study of extreme events. The variables of "cumulative rainfall" (PC n ) i.e., in 1, 2, 5, 10, 30, 60, 90, 120 and 180 days are computed from daily rainfall data by means of the expression with n = 1, 2, 5, 10, 30, 60, 90, 120 and 180, and where j is the progressive number of days that form the analyzed time interval and P i is the rainfall value recorded the ith day.
The maximum values of each variable are extracted, year by year, from the data sets generated in this way, and the parameters of the GEV function are determined from the above values, by applying the probability weighted moments (PWM) method introduced by Greenwood et al. (1979) and subsequently modified by Hosking et al. (1985). Finally, the inversion of the probability function yields the values of cumulated rainfall x for each of the variables (1, 2, 5, 10. . . 180 days) and for different RPs. Then, these values are interpolated with a view to build the rainfall probability curves.
To yield reliable results, this type of analysis requires sufficiently long and continuous time series of rainfall data (at least 20 years of recorded data according to Houghton et al., 2001 andSerrano, 2010). For this reason, daily rainfall data were obtained from Santo Stefano di Briga and Messina Istituto Geofisico rainfall stations, that have been operational since 1925 and 1952, respectively. However, if we consider that the extreme rainfall events which periodically affect the study area are usually of short duration, as in the case of the 1 October 2009 event, it would be extremely interesting to analyze the historical data of maximum hourly rainfall intensity. Unfortunately, these data are not available for the stations mentioned above. For this reason, hourly rainfall data (i.e., cumulated in 1, 3, 6, 12 and 24 h) were used from Alì Terme station that is located approximately 4 km SE of Fiumedinisi station (see Fig. 12b); thus, this station has been considered sufficiently close to be used to assess the recurrence of the rainfall events recorded in Fiumedinisi station.

Back analysis of the 1 October 2009 event
As previously mentioned, before applying the TRIGRS model to back-analyze the 1 October 2009 event, it was necessary to evaluate the soil moisture conditions prior to the event through the HYDRUS 1-D model. According to the simulation results obtained by using the grain size distribution of sample no. 1 (see Table 1), the absence of a steady-state water table within the soil cover can be assumed, whereas in Fig. 8a the resulting volumetric water content trend with depth at four different times (1, 24, 25 and 30 September) is reported. The initial soil moisture (θ = 0.047) is assumed to be near to the residual water content value considering the hot, dry conditions during the preceding summer months. The effect of the September rainfall (Fig. 8b) results in a progressive increase in soil water content that is equal to 0.194 on 25 September (the day after the second rainfall event; see Table 2). It is worth noting that in the lower part of the soil profile the water content progressively increases over the days, due to the advance of the wetting process, resulting in an average value of 0.19 on 30 September, which corresponds to a gravimetric water content (w) of approximately 10.8 % and a degree of saturation (S r ) equal to 54.4 % (on the basis of the physical properties reported in Table 1).
Once the initial soil moisture conditions are estimated, all the input parameters required by TRIGRS can be defined (Table 3). A detailed (2 × 2 m) pre-event digital elevation model (DEM) was used, substantially resampled at the 4 × 4 m resolution due to limitations on computing time. Soil thickness (H ) and rainfall rate (I Z ) vary from cell to cell on the basis of maps obtained through the methods described in Sect. 4.2. According to the available data, an average friction angle of 35 • was used, whereas the cohesion was progressively increased to 3 kN m −2 , so that only very few cells (i.e., 260, which represent about 0.04 % of all grid cells in the study area) result in instability before the beginning of the event. In any case, the chosen value lies within the range of values reported by other authors for the same material (see Sect. 4.2). The depth-averaged soil unit weight (γ n ) is equal to 19.22 kN m −3 given the porosity (35 %), the degree of saturation (54.4 %) and the unit weight of soil solids (26.73 kN m −3 ), while θ s (saturated water content), θ r (residual water content) and K s (saturated hydraulic conductivity) are directly predicted using the HYDRUS-1D model. Given the absence of an initial water table, its depth (d wt ) corresponds in this way to the bedrock-soil interface. To evaluate the α G parameter, that is typical of the Gardner hydraulic model, the conversion formula that was introduced by Ghezzehei et al. (2007) was used, which defines a correspondence between Gardner and van Genuchten-Mualem models through the capillary length approach (Warrick, 1995). On the basis of the results of the same simulations, the I ZLT parameter, that represents the long-term background rainfall rate, was assumed to be equal to the cumulative actual surface flux value (8.49 × 10 −8 m s −1 ). Finally, the saturated hydraulic diffusivity (D 0 ) has been calculated according to where K s is the saturated hydraulic conductivity, H the average soil thickness (80 cm) and S y the specific yield (Grelle et al., 2014). If we consider that the investigated soil can be classified as loamy sand, the specific yield is assumed to be equal to 0.34, on the basis of typical values given by Johnson (1967) (also reported in Loheide II et al., 2005) for each soil textural class. With regard to the comparison between the numerical simulations and the landslide inventory map, Table 4 reports, as well as the number and relative percentage (P U ) of predicted unstable pixels (i.e., SF ≤ 1), the percentage of correctly pre- dicted landslide (P L ) and stable (P S ) pixels between 13:00 and 21:00 UTC. Figure 9 shows the temporal evolution of slope instability at the catchment scale: only 0.2 % of pixels are indicated as unstable at 13:00 UTC (beginning of the rainfall event). After a progressive increase in the following 4 h, the instability rapidly rose between 17:00 and 18:00 UTC (P U : +21.2 %) in correspondence of a rainfall peak, and the critical stage continued until 21:00 UTC, given that P U passed from 39.6 % to 100 % in just 3 h. This temporal evolution of the phenomenon substantially agrees with the witnesses, although during the real event no particular increase of slope instability has been recorded after 20:00 UTC. To fully evaluate the accuracy of the model, a ROC (receiver operating characteristic) curve analysis was performed by comparing the final SF map (21:00 UTC) with the landslide inventory. The ROC curve measures the goodness of the model prediction plotting, for different threshold values, the true positive rate, i.e., the proportion of correctly predicted positive values ("landslide presence") and the false positive rate, i.e., the proportion of negative values ("landslide absence") erroneously reported as positive. The area under the curve (AUC), which varies from 0.5 (diagonal line) to 1, quantifies the predictive capability of the model. According to the results, the SF map correctly classifies 68.8 % of source areas (True Positive) and 75.1 % of stable areas (True Negative) with SF = 1, whereas the AUC is equal to 0.795 (Fig. 10a). However, due to the great variability of the soil texture, three different grain size distributions are available (Table 1), but the back analysis has been performed using the grain size characteristics of sample no. 1 only. This choice can be explained analyzing the simulation results obtained accounting for the grain size distributions of the other two samples. With regard to the parameters derived from HYDRUS 1-D (Table 5), it can be noted that the saturated hydraulic conductivity and consequently the hydraulic diffusivity can differ by up to an order of magnitude, due to the presence of a larger quantity of fine material (silt and clay) in sample no. 2 and no. 3. In addition, the other parameters are substantially similar for all three of the samples. By using these parameters in TRIGRS simulations, a different reconstruction of the 2009 event is obtained (Fig. 10b). Although the maximum number of unstable pixels is similar in the case of sample no. 1 and sample no. 2, in the sec-    ond one, the instability peak occurs far too late compared to the real event (approximately at 05:00-06:00 UTC on the following day). Furthermore, for sample no. 3 the low value of hydraulic conductivity causes not only an even greater delay of the instability, but also a lower number of unstable pixels at maximum instability, which develops in a larger time interval (approximately between 09:00 and 13:00 UTC of 2 October). For these reasons, we consider the parameters obtained according to the grain size characteristics of sample no. 1 as the most representative of the investigated soil. Figure 11a and b show a graphic comparison between cumulative frequency (symbols) and GEV probability function (continuous line), obtained by using the daily rainfall records from Santo Stefano di Briga and Messina Istituto Geofisico stations. As it can be observed, the good fitting between data and probability function confirms the reliability of the applied method. With regard to the probability curves ( Fig. 11c  and d), the comparison reveals that the highest rainfall values are attributed to the Santo Stefano di Briga curves for the same RP. This finding emphasizes that, in the past, this station (the most representative of the sector most severely hit by the 2009 event) has recorded more intense and severe rainfall events than the other one. On the basis of the same curves, the RPs of the rainfall accumulated up to 1 October 2009 have been estimated (Table 6). An estimation has been made also for the rainfall accumulated up to 30 September (i.e., the day prior to the event), but the obtained values infer that the rainfall amount, at both stations, is far from exceptional (estimated RP of 1 year); thus, rainfall prior to the event practically lies within the standard range, in contrast with rainfall accumulated up to 1 October. In this case, while rainfall recorded at Messina Istituto Geofisico continues to be unexceptional (estimated RP of 4-5 years), rainfall accumulated in a single day (1 October) at Santo Stefano di Briga has a RP of 47 years. This means that the event under study was not only strongly localized in space, but also particularly severe in that specific sector. This finding is also substantiated by what has been previously pointed out, i.e., the highest RPs have been obtained for the station with the highest rainfall probability curves. With regard to the analysis of the historical data of maximum hourly rainfall intensity from Alì Terme station, the results show that the fit between cumulative frequency and Saturated hydraulic conductivity 6.6 × 10 −5 1.25 × 10 −5 7.91 × 10 −6 α G (m −1 ) Gardner parameter 11.8 11.1 12.2 I ZLT (m s −1 ) Background rainfall rate 8.49 × 10 −8 5.35 × 10 −8 5.37 × 10 −8 D 0 (m 2 s −1 )

Rainfall probability curves and return period of the 1 October 2009 event
Hydraulic diffusivity 1.55 × 10 −4 3.84 × 10 −5 2.43 × 10 −5   (Fig. 12a) is not as good as in the preceding analyses. However, it is worth stressing that data for intense precipitation are generally scantier than daily values and that the resulting statistical analyses are usually less reliable. The resulting rainfall probability curves (Fig. 12b) define a RP of 78 years for the 1 h rainfall recorded on 1 October 2009 in Fiumedinisi station, a value greater than that estimated for the 1-day rainfall recorded at Santo Stefano di Briga station (47 years). Therefore, the 1 h rainfall event can be classified as an extreme event. Nevertheless, it is particularly interesting to analyze the sub-event of maximum duration equal to 3 h, a time after which major damage was observed in the area. In this case the estimated RP is equal to 26 years (Table 7); this value infers that even if the RP of the 1 h rainstorm shows that it was an extreme event, the 3 h sub-event is characterized by a much lower RP, which suggests its classification to be not severe. However, it is worth noting that the probabilistic analysis is affected by several uncertainties, related to the type of probabilistic model and the definition of the parameters of the model itself. Generally speaking, the uncertainty tends to increase when the sample size (i.e., the number of measurement years) is decreased and the relevant RP is increased.

Evaluation of different triggering scenarios
Once the physically based model was validated through the back analysis of the 1 October 2009 event, different TRIGRS simulations were performed by varying the rainfall input, on the basis of the rainfall probability curves described above. With regard to the other input parameters of the model, those used for the analysis of the 1 October event have been kept. It is important to note that only hourly simulations have been performed. In fact, if we reproduce the 2009 event and compare the simulation results obtained by using the 1-day rainfall value with those gained with the 15 min rainfall maps, it results in a much lower number of predicted unstable pixels in the first case (5405) compared to the second (126 207). Therefore, daily data are not suitable for this type of analysis, due to the excessive underestimation of the instability phenomena. In the hourly simulations, four rainfall values, which correspond to four different RPs (i.e., 2, 4, 10 and  20 years) have been used, in such a way as to evaluate the effect of more frequent rainfall events with respect to the 2009 event. Table 8 shows the number of unstable pixels at the end of the simulated rainfall event and in the next 2 h. In this way, it can be noted that the increase of instability continues up to 1 h after the end of the rainfall, specifically for infiltration rates not less than 15 mm h −1 . This aspect can be explained in that water needs a certain time to reach greater depths, so that even the deeper and flatter grid cells may become unstable. The results also emphasize the importance, in terms of produced instability, of the 1 and 3 h rainfall amount even for relatively low RPs (i.e., 10-20 years). For instance, 3 h rainfall with a RP of 10 years would cause about half (48.9 %) of the slope instability calculated for the entire 1 October event. In the case of 6 h rainfall, a significant level of instability develops only for events with RP ≥ 20 years, while with events of longer duration (12 h) the number of predicted unstable pixels is extremely low, even for relatively high RPs. Therefore, the duration of a rainfall event can produce completely different triggering scenarios for an equal rainfall amount. In this respect, we used two different rainfall inputs (i.e., 45 and 85 mm) whose RP varies depending on the duration of the event (i.e., 1 or 3 h), in order to analyze the variation of SF and the pressure head over time at the same grid cell of interest (Fig. 13). The elevation and slope value of the cell (i.e., 320 m a.s.l. and 38 • ) correspond to the average values observed within the landslide source areas.
The results shows that the SF decrease (and consequently the increase of the pressure head measured at the same depth) is greater and more rapid in the case of 1 h rainfall events, due  to the greater rainfall rate. Furthermore, also in this analysis the increase of instability occurs in the hour following the event. In this sense, the maximum variations are observed in the 85 mm/1 h rainfall, in which the SF and the pressure head pass from 0.96 to 0.9 and from 0.08 to 0.14 m, respectively. Furthermore, in all the four cases a very slow recovery of the preceding stability conditions can be noted over the 2 h after the instability peak.

Discussion
The results obtained from the numerical simulations confirm not only the need of hourly analyses for extreme but short- 3). In fact, considering the high values that characterize the rainfall probability curves of Alì Terme rain gauge station (see Sect. 5.2), we can assert that short-duration rainfall events frequently have a high intensity in this specific area. Therefore the combination of recurring and heavy rainfall events, probably due to specific geomorphological and climatic features that influence the development of localized severe storms, justifies the approximately 40 landslide/flood events that have occurred since the last century in this area. From a methodological point of view, it is worth noting how the approach proposed in this paper (Fig. 14) combines various techniques and methods, optimizing different types of data depending on their availability. For instance, the parameterization of the physically based model can be performed both in the absence and presence of preceding reference events. In the former case, only the geotechnical parameters and the soil thickness are needed, while in the latter, the process used for the back analysis of the 1 October 2009 event can be applied to any other event. Here, in particular, it is important to stress that the comparison with a preceding landslide event allows the reliability of the model to be increased, as long as a comprehensive and detailed event-based landslide inventory exists. With regard to the evaluation of the initial soil conditions, in this study the HYDRUS 1-D model has been used in relation to the 1-month antecedent rainfall. However, some recent studies have investigated the linkage between soil moisture and landslide occurrence by using soil moisture data derived by in situ (Baum and Godt, 2009;Hawke and McConchie, 2011) and satellite sensors (Ray and Jacobs, 2007;Ray et al., 2010), and this type of measurements, if available, can be used to define the input parameters of the model. The last step of the approach concerns the definition of the rainfall input to be used for the evaluation of a specific triggering scenario. By means of the statistical analysis of hourly rainfall data, different rainfall values with different RPs may be used to depict different scenarios, and change the initial soil conditions and the duration of the rainfall input. However, as emphasized in the discussion of the results, establishing which is the critical rainfall duration that triggers the shallow landslides cannot be straightforward, because completely different slope stability conditions can be obtained with different combinations of single rainfall inputs within the same rainfall event, even with the same initial soil conditions. For this reason, and considering the chance to use also soil moisture data, a possible application of the approach proposed here could be to develop an early warning system based on rainfall thresholds, identified by the physically based model validated according to the process described above. In this way, a model in which the initial conditions are constantly updated could depict more consistent and reliable triggering scenarios by using any rainfall event forecasted for the next hours.

Conclusions
In this study, we introduce an approach for the analysis of shallow landslide-triggering scenarios that uses the TRIGRS code, a physically based model which describes the stability conditions of natural slopes in response to specific rainfall events. As a first step, the model has been validated through the back analysis of a reference landslide event, i.e., the disaster that occurred in the southern Messina area on 1 October 2009. Comparing the results of the numerical simulation with the 2009 landslide inventory, it turns out that the model is able to reproduce the reference event quite well, both in terms of temporal evolution and spatial distribution of slope instability, identifying the areas most affected by shallow landslides. It is worth stressing that the model has been accurately parameterized through different methods and techniques, with specific focus on the evaluation of the spatial pattern of the triggering storm and initial soil moisture conditions.
Once the physically based model has been validated, different triggering scenarios have been reconstructed by varying the rainfall input, on the basis of the rainfall probability curves obtained through a statistical analysis of historical rainfall data. The results indicate that the 1-day simulations strongly underestimate the landslide events triggered by extreme but short-duration rainfall events (such as the 1 October one). With regard to the hourly analyses, it results that even recurring (with a RP of 10-20 years) 1 and 3 h rainfall events can lead to significant slope instabilities. This feature confirms the destabilizing effect of recurring rainfall events in the study area, justifying the high number of landslide/flood events that have occurred in the past. With regards to the proposed approach, the use of different techniques allows its application to different case studies, on the basis of the data availability. Furthermore, if we consider the possibility of depicting constantly updated triggering scenarios, this approach could be used to develop specific landslide early warning systems in order to support decision makers in both risk prevention and emergency response.