Shallow landslide initiation on terraced slopes: inferences from a physically based approach

ABSTRACT In the last years, great efforts have been made to improve the assessment of the temporal and spatial occurrence of rainfall-induced shallow landslides. Therefore, in this paper we used a physically based stability model (TRIGRS) in order to reproduce the landslide event occurred in the Monterosso catchment (Cinque Terre, Eastern Liguria, Italy) on 25 October 2011. The input parameters of the numerical model have been evaluated taking into account the land-use setting and paying specific attention to the evaluation of the spatial variation of soil thickness on terraced areas. The resulting safety factor maps have been compared with the inventory map of the landslides triggered during the event. The simulation results, which have been obtained also considering four different spatial resolutions of the digital terrain model, emphasize the influence of land use in shallow landslide occurrence and indicate the importance of a realistic spatial variation of soil thickness to enhance the reliability of the model. Finally, different triggering scenarios have been defined using hourly rainfall values statistically derived from historical data. The results indicate the proneness of the area to shallow landsliding, given that rainfall events with a relatively low return period (e.g. 25 years) can trigger numerous slope failures.


Introduction
Shallow landslides triggered by rainfall are a widespread hazard that frequently results in considerable damage to infrastructure and human losses in many mountainous regions of the world, especially in areas characterized by the widespread presence of soil cover and subject to heavy precipitation (Aleotti and Chowdhury 1999;Dai et al. 2002;Lin et al. 2006). In order to predict the occurrence of such events, in the last years many attempts have been made to establish a relationship between rainfall and landslides, in particular by means of physically based numerical models. Despite numerous slope stability models have been developed (e.g. Montgomery and Dietrich 1994;Iverson 2000;Borga et al. 2002;Liao et al. 2011;Formetta et al. 2014;Ho and Lee 2016;Thiery et al. 2017), the fundamental controls leading to slope failure driven by rainfall are still not well quantified (Borja et al. 2012), and thus the improvement of current models is still an important research topic (Chang et al. 2008). The difficulty of building up reliable mathematical models lies in the numerous variables involved in the triggering process, such as spatial and temporal rainfall variability, mechanical and hydraulic soil properties, slope morphology, vegetation coverage, initial soil suction and moisture (Greco et al. 2010). However, since it is challenging to comprehensively analyse all these parameters, it makes sense to focus on those factors that mainly contribute to the onset of the instability. In this respect, land use/land cover is widely considered as a relevant factor for the triggering of shallow landslides (Glade 2003), since it strongly influences different parameters, such as hydrological and mechanical properties of the soil, thickness of the debris cover. Different studies have been performed to evaluate the role of land use in shallow-landsliding, mainly through the use of statistical methods (e.g. Baeza and Corominas 2001;Beguer ıa 2006;Piacentini et al. 2012;Galve et al. 2015;Trigila et al. 2015;Persichillo et al. 2016). These methods are based on conceptual models that describe the functional relationships between instability factors and the past and present distribution of slope failures (Carrara 1983). Bivariate and multivariate approaches are the most common statistical techniques, though probabilistic (e.g. Lee 2005;Goetz et al. 2011Goetz et al. , 2015Youssef et al. 2016;Castro Camilo et al. 2017) and Machine-Learning approaches (e.g. Lee et al. 2004;Ermini et al. 2005;Marjanovi c et al. 2011;Pham et al. 2017) are also used in this type of analyses. Such approaches are very sensitive to the type and quality of the factors chosen for the susceptibility analysis, and the lack of suitable expert opinion can produce unreliable results (Soeters and Van Westen 1996). Additionally, being data-driven, a statistical model built up for one region cannot readily be extrapolated to the neighbouring area. On the contrary, the main drawback of physically based modeling is the difficulty to gather the input parameters over large and complex areas (Carrara et al. 2008;Chang and Chiang 2009). However, a possible solution is to calibrate the parameter values through back-analysis of preceding major landslide events, after which event-based landslide inventories are generally available (Casadei et al. 2003;Li et al. 2011). In addition, the possibility to consider different triggering inputs makes physically based models particularly suitable for depicting different landslide scenarios, which represents the first step toward a proper hazard analysis.
For these reasons, a well-established model, as TRIGRS (Baum et al. 2008), was used. This model predicts the shallow landslide occurrence by means of the transient, one-dimensional (1D) analytic solution for pore-pressure response to rainfall infiltration with an infinite slope stability calculation (Savage et al. 2003). In particular, the model was applied for the reconstruction via back-analysis of the landslide event occurred on 25 October 2011 in a coastal basin located in Eastern Liguria (Northern Italy), for which a comprehensive event-based landslide inventory has been specifically prepared. For the parameterization of the numerical model, we used a multi-methodological approach recently applied to a similar case study in Southern Italy (Schilir o et al. 2015a(Schilir o et al. , 2015b; however, we properly modified this approach in order to take into account the specific landuse of the study area, characterized by the widespread presence of agricultural terraces. Specifically, we performed several numerical simulations considering two different methods for the evaluation of the spatial variation of the soil thickness and four different spatial resolutions of the digital terrain model. In addition, since the investigated area has already been affected by flood/landslide events in the past, after evaluating the predictive capability of TRIGRS through comparison with the 2011 landslide inventory we also performed different numerical simulations for depicting future triggering scenarios, on the basis of rainfall values resulting from a statistical analysis of the historical rainfall data.

General features of the study area
The study area is located in Eastern Liguria, along the coast of Cinque Terre (Italy) (Figure 1(a)). The Cinque Terre are renowned worldwide as a typical example of man-made landscape, characterized by widespread century-old agricultural terraces retained by dry stone walls. For this reason, this area was declared 'World Heritage Site' by UNESCO in 1997 and successively it was included in the 'Cinque Terre National Park'.
From a geological point of view (Figure 2(a)), the bedrock is mainly composed of a sandstoneclaystone flysch (Riomaggiore Banded Sandstones lithofacies, Macigno Fm., Tuscan Nappe) with rare claystones with limestones and silty sandstones turbiditic rocks (Canetolo Shales and Limestones Fm., Sub-ligurian Domain). Towards W, pelitic rocks (Mt.Veri Shales, External Ligurides), jaspers (Mt. Alpe Jaspers Fm., Internal Ligurides), serpentinites and gabbros (Internal Ligurides) outcrop, respectively (Abbate 1969;Regione Liguria, 2005). On slopes, the sedimentary bedrock is covered by a thin (0.5-2.5 m) eluvial-colluvial soil cover, usually characterized by mixtures of gravel and sand with minor components of silt and clay (Cevasco et al. 2013b(Cevasco et al. , 2014; instead, the ophiolitic slopes often lack of soil mantle or it is very thin . The shallow soils and debris covers were largely reworked in the past centuries by humans for the building of terraces. The geomorphological features of the study area, which are controlled by the vicinity of the watershed to the coast, are typical of most coastal basin of the easternmost Liguria. In particular, these catchments have small area (few square kilometres or less than 1 km 2 ), high gradient of slopes and short streams with ephemeral hydrological regime. More than 40% of the slope angles are ranging between 30 and 40 . The characteristics above, along with the urbanization of the Monterosso village, which mainly developed on the floor of the Rio delle Rocche and Canale Pastanelli valleys, make Monterosso particularly exposed to flood risk. In addition, the conditions imposed by the coverage of the final tracts of both the Rio delle Rocche and Canale Pastanelli, for about 400 m upstream of the marina, increase the risk of flooding (Figure 1(b)). Peculiar land-use features, typical in Cinque Terre and Liguria, characterize the study area. Despite the slopes were almost completely terraced during the past centuries (about 71% of the study area) for vineyards and oliveyards, only 19% of the total terraced areas (about 14% of the study area) is still cultivated (Figure 2(b)). A very similar ratio between cultivated terraces vs. total terraced area characterizes the Vernazza catchment (Cevasco et al. 2013a), close to the study area ( Figure 3). This point indicates how the social changes occurred in the second half of 1950s affected the land use at Cinque Terre. On abandoned terraced areas, vegetation spreads quite rapidly. At first, a vegetation consisting mainly of scrubs grows up sparsely. After a time span of about 25-30 years, a dense vegetation constituted by forest tree species (mainly pine, mixed woods, chestnut woods, holm oak woods) or by Mediterranean scrub covers abandoned terraces .
The climate of the study area is Mediterranean, particularly mild due to geographical and morphological peculiarities. In general, summers are hot and dry, winters are mild and autumn is the rainiest season. Based on the analysis of 1954-2012 data from the Levanto rain gauge (4 km NW from the study area), the mean annual precipitation is 1033 mm, whereas the rainiest month is October, with a mean precipitation of 154 mm . The mean annual precipitation tends to increase with altitude moving from the coast to the watershed (Pedemonte 2005). In particular, in the proximity of the culminations of the watershed bounding the study area to the N the mean annual precipitation reaches about 1200 mm.
Due to self-regenerating storm cells or to persistent cyclonic Tyrrhenian circulation (Crosta 1998;Cevasco et al. 2009), extreme precipitation can occur, between late summer and mid-autumn, along the coast of eastern Liguria (van Delden 2001;Silvestro et al. 2012Silvestro et al. , 2015Rebora et al. 2013;Buzzi et al. 2014;Cevasco et al. 2015;Cassola et al. 2016). Usually these rainstorms have short duration (<24 h), but rainfall intensities can reach or exceed values of 100 mm h ¡1 . 2.2. The 25th October 2011 event and previous flooding events On 25 October 2011 extreme precipitation hit a wide area between east Liguria and north-west Tuscany. In particular, a violent storm system with a 'self-healing' structure originated over the Cinque Terre area in very short time, caused extreme rainfall between 7:00 and 17:00 (Universal Time Coordinated, UTC) along the coast and the inland Vara valley (east Liguria), and then it moved to the Lunigiana (north-western Tuscany). Based on data from Regione Liguria (2012), at the Monterosso rain gauge a cumulative rainfall of 382 mm was recorded, with maximum rainfall intensities reaching 90 mm/h, 195 mm/3 h and 350 mm/6 h ( Figure 4). Higher values, regarding both the cumulative rainfall and the rainfall intensities, were recorded at the Brugnato rain gauge (Vara valley, about 10 km NE from the study area), with cumulative rainfall of 538 mm and maximum rainfall intensities reaching 153 mm/h, 328 mm/3 h and 472 mm/12 h, respectively. Hundreds of shallow landslides and floods affected the Tyrrhenian coastal basins between Bonassola and Manarola and the Magra river basin (Cevasco et al. 2012(Cevasco et al. , 2013aD'Amato Avanzi et al. 2013, causing 13 casualties and very severe damage to villages and infrastructures. In the study area, the Monterosso village was affected by a catastrophic mud/debris flood. Enormous amount of earth and debris, mobilized on slopes and charged by torrents downhill, overwhelmed the central streets of Monterosso as consequence of the filling of the final culvert of the Canale Pastanelli. The maximum thickness of deposited materials reached about 3 m in the centre of the village. One casualty occurred and shops and restaurants were devastated. After the 25th October 2011 rainstorm, a detailed rainfall-induced landslide inventory was specifically compiled for the study area through analysis of high-resolution aerial photographs and field surveys. The air photo analysis was carried out on digital geo-referenced orthophotos, provided by Liguria Regional Administration, taken on 11 November 2011 by the Air Service of Remote Sensing and Monitoring of Civil Protection of Friuli Venezia Giulia Regional Administration (ground resolution from 3 to 50 cm, according to the altitude). The high resolution of aerial photographs allowed to achieve great detail in the mapping of landslides.
A total of 260 shallow landslides were mapped in the study area, more than 214 of which occurred in the Rio delle Rocche and the Canale Pastanelli catchments. The shallow landslides triggered on 25 October 2011 initiated as debris slides (Cruden and Varnes 1996) and in most of cases evolved into debris avalanches or, in few cases, into debris flows (Hungr et al. 2014). Shallow earth and debris covers (from 0.5 m to about 2.5 m thick) and, sometimes, portions of the weathered and fractured bedrock were involved. According to the landslide classification by Hungr et al. (2014), we considered as debris avalanches 'very rapid shallow flows of partially or fully saturated debris on a steep slope, without confinement in an established channel' and as debris flows 'very rapid to extremely rapid flow of saturated non-plastic debris in a steep channel'. The total area affected by landslides was about 6.5 ha, corresponding to 1.2% of the study area, whereas the average density of landslides was 47.6 landslides/km 2 . Landslide areal extent ranged between few tens of square metres up to few thousands square metres. Earth and debris mobilized on slopes or along channels by shallow landslides increased the sediment transport along the drainage network and also the erosive energy of streams. The inability to dispose flowing of water/sediment mixtures of the final culvert of the Canale Pastanelli, played a fundamental role in originating the mud/debris flood that affected the Monterosso village.
However, in the last century, another disastrous flood occurred on 2 October 1966 at Monterosso. Data about this event were collected on the 'Italian archive of historical information on landslides and floods' (Guzzetti et al. 1994;Guzzetti and Tonelli 2004), Italian newspapers and local chronicles. The flood occurred in the evening of the 2 October 1966. The central streets of Monterosso were flooded by water and earth, spreading on basements, shops, bars and sweeping away towards the sea the boats standing on the beach. Newspapers reported images showing the central street of the village ( Figure 5(a)), as well as the area between the railway and the marina, flooded by water and debris. From the comparison of these images with those depicting the same areas after the 25 October 2011 flood one can deduce that the effects of the two floods were of comparable size ( Figure 5(b)). Another sign testifying the similarity between the effects of these two floods was imprinted on the wall of the S. Giovanni Battista Church, located just upstream the marina, indicating the maximum level (about 2.5 m) reached by water and debris in both floods ( Figure 5(c)).
The flood event occurred in 1966 caused severe damage to people, buildings and infrastructures. In particular, 50 displaced persons and 300 homeless were estimated (AVI database, http://sici.irpi. cnr.it/). The same database reports that rainfall lasted 14 h continuously. Since no rainfall data for that event were available for the Monterosso rain gauge, we referred to the nearest Levanto rain gauge. The rainfall started on 30 September 1966 and lasted until 4 October 1966. A total rainfall of 263 mm was recorded in these five days but 142 mm were recorded just from 2 October 1966 to 3 October 1966. Maximum intensity rainfall of 85 mm/6 h, 127 mm/12 h and 152 mm/24 h were recorded during that event. By comparing rainfall data of 1966 and 2011 floods, it results that the 2011 event was the most severe.

Methodology
The work described in this paper can be divided in three stages. First, we performed a backanalysis of the 25 October 2011 event by comparing the resulting triggering scenarios with the inventory of the landslides occurred during the event. Afterwards, we performed a statistical analysis of the historical rainfall data available for the study area in order to evaluate the recurrence of specific rainfall events. Finally, once the physically based model was validated through the back analysis of the reference landslide event, we reconstructed different triggering scenarios by using rainfall inputs, whose return periods (RPs) have been defined through the preceding statistical analysis.

Theoretical basis of TRIGRS model
TRIGRS (Transient Rainfall Infiltration and Grid-based Slope Stability model) is a physically based model devoted to the spatial and temporal prediction of rainfall-induced shallow landslides over large areas. The model performs an infinite-slope stability calculation accounting for the pressure head response C(Z,t) to a time-varying rainfall input on the ground surface I Z (t). Specifically, the infiltration process is simulated considering the 1D analytic solution of Richards' equation as well as described by Iverson (2000). In the most recent version of the program (Baum et al. 2008), TRIGRS was expanded to address infiltration also into unsaturated soils, then assuming a two-layer system consisting of a saturated zone (with the possible presence of a capillary fringe) and an unsaturated zone that extends to the ground surface ( Figure 6). For the linearization of Richards' equation within the unsaturated zone, TRIGRS considers the Gardner (1958) hydraulic model, which is based on the definition of four main hydrodynamic parameters: the saturated (u s ) and residual (u r ) water content, the saturated hydraulic conductivity (K s ) and a specific parameter linked to the pore size distribution (a G ). According to the hydraulic model, the unsaturated zone absorbs part of the water that infiltrates through the ground surface, while the remaining water accumulates above the initial water table. If the amount of infiltrated water 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. On the basis of the new water table level, the model then computes the new pressure head value (Baum et al. 2010). Finally, the Safety Factor at depth Z is calculated on the basis of a series of input parameters locally variable throughout the model: where d is the slope angle, c' is the soil-cohesion for effective stress, '' is the soil friction angle for effective stress, g w is the unit weight of groundwater and g is the total unit weight of soil.

Evaluation of TRIGRS input parameters and assumptions in the light of a preceding case study
In view of back-analysing the 25 October 2011 event, we evaluated the input parameters of TRIGRS through an approach based on the application of different methods and techniques. This approach has been already tested (Schilir o et al. 2015a; in a similar case study occurred in Southern Italy on 1 October 2009. On that day, a heavy rainstorm triggered several hundreds of shallow landslides in the southern Messina area (north-eastern Sicily), causing 37 fatalities and severe damage to buildings and infrastructure, mostly near the village of Giampilieri (Schilir o and Esposito 2013). The choice of using the same approach can be then explained considering the similarities between the two case studies. In fact, the geological and geomorphological contexts are very similar given the presence, in both areas, of a thin eluvial-colluvial soil cover above an impermeable bedrock and of small coastal basins having steep slopes with torrent-like straight watercourses. Furthermore, the triggering events are comparable to each other, being both related to a violent storm system with a 'self-healing' structure which caused extreme, short duration rainfall events with extremely high intensity peaks. The typology of triggered landslides is substantially the same ( Figure 7) and probably reflects similar triggering conditions, also considering the similarities between the involved materials, both characterized by high content of coarse grained particles with minor components of silt and clay (Table 1). Precisely in order to investigate the failure mechanisms of the soil cover in a wide range of initial conditions (i.e. porosity and water content), in a recent study ) specific laboratory flume tests have been conducted on the material sampled near Giampilieri, and the development of the failure process was also numerically simulated with SLIP, a simple physically based model that predicts the triggering of shallow landslides on the basis of the soil   characteristics and the rainfall amount (Montrasio 2000;Montrasio and Valentino 2008). The results of the experimental tests substantially agree with the outcome of the numerical simulations, since the mean variation between real and predicted failure time is in the order of a few minutes, even for the longest tests . Furthermore, it is interesting to note that not only the shallowest part, but also the inner part of the soil profile appeared saturated after the failure ( Figure 8). This feature fits with the hydrological model implemented in SLIP, which assumes a water flow through the soil macropores that generates a progressively saturation of non-adjacent volumes of soil (Zhang and Zhang 2009): if the raining process persists, the saturated portions of soil extend and become continuous ( Figure 9) and the sliding process starts when a relatively wide continuous stratum of saturated soil has formed. Although no measurements of pore water pressure have been recorded during these tests, in previous works (i.e. Montrasio andValentino 2007, 2008), experimental measurements of suction and pore pressure have been recorded during flume tests under different rainy conditions, confirming that shallow soil layers reach a zero value of suction well before the deeper ones and the first slip occurs only after the matric suction has dropped to zero everywhere. This condition corresponds to a state of complete saturation of the soil layer that is in contact with the impermeable bed (Montrasio and Valentino 2007). The reliability of the SLIP model was also substantiated by the results obtained for the back-analysis of the 2009 event at the catchment scale . In this respect, it is worth noting that these results are similar to those obtained with TRIGRS: this point suggests the similarities between the two models, despite a different way to simulate the infiltration process. However, it is important to stress that SLIP bypasses the transient analysis of the infiltration process assuming the formation, as a consequence of a rainfall event, of the final condition of a perched water table akin to that resulting from the TRIGRS simulations. Therefore, on the basis of the remarks pointed out above, we have defined the TRIGRS input parameters (Table 2) by using different types of data available for the Monterosso area. A detailed (1 £ 1 m) pre-event digital elevation model (DEM), deriving from airborne LiDAR data, has been used as basic information for the numerical simulations. The original DEM has also been resampled to 2, 4 and 8 meters by using the bi-cubic polynomial interpolation technique, with the aim of evaluating the simulation results at different resolutions. Accurate geological and land-use maps (Figure 2), resulting from air photo analysis, field surveys and site investigations, have been used for the evaluation of the spatial variation of some specific parameters. For instance, the estimation of soil thickness ( Figure 10) has been performed by correlating soil depth to the local slope angle (Saulnier et al. 1997): where h i and a i are, respectively, the soil thickness and the slope measured at pixel i, whereas h max , h min , a max and a min are the maximum and minimum thickness-slope values encountered in each terrain unit, identified from the overlap of the geological and land-use maps ( Figure 11). The soil thickness ranges have been defined according to field measurements performed in the neighbouring Vernazza basin (Cevasco et al. 2013b(Cevasco et al. , 2014, since this basin has quite similar geological and geomorphological features to the study area. However, the spatial variation of the soil thickness can be Figure 9. From top to bottom: evolution of saturation process through the soil macropores (from Montrasio and Valentino 2008) highly influenced by the widespread presence of agricultural terraces, which represent one of the most relevant peculiarities of the Cinque Terre area. For this reason, we performed an alternative soil thickness estimation by using a morphometric approach. Since the available DTM is enough detailed to reproduce the morphological variations due to agricultural terraces, the basic idea was to  analyse the stepped slope morphology by means of quantitative morphometric indexes to account for the 'cyclic' soil thickness variation that occurs moving from the downslope-facing wall of a terrace (maximum thickness) to the innermost part of the flat slope (minimum thickness) featuring the proper terrace. For this purpose, the Topographic Position Index (TPI) algorithm (Weiss 2001) has been tested to analyse the terrain morphology at a very low wavelength, with the final aim of zoning the terraced areas in terms of soil thickness. The TPI algorithm performs an automated landform classification by comparing the elevation of each DTM cell with the surrounding ones. The resulting values are strongly dependent on the size of the considered neighbourhood and the shape of the search window. Starting from the TPI values, it is possible to partition the territory according to the Slope Position Index (SPI), which provides for the classification of the TPI continuous values according to specific threshold values that are based on the standard deviations of TPI. Assigning a SPI value to each cell allows to define its relative position with respect to the topographic structure, i.e.for examplewhether a DTM cell element belongs to a flat area, or represents a portion of a ridge, a valley or a slope (low, middle and upper). For this case study, after several attempts to find the 'best fit' search window in terms of both radius and shape, it was possible to select a circular window with a 6 m radius. The results have been validated by comparison with the actual topography. Furthermore, after having classified the TPI values in terms of SPI, it was possible to partition the terraced slopes in zones (namely 'valleys', 'lower slope', 'mid slope', 'upper slope', 'ridge' and 'flat slope') ( Figure 12) having different average soil thicknesses. The values of soil thickness related to each SPI class have been estimated based on the knowledge of the technique used to build the agricultural terraces. The resulting map is represented in Figure 13: however, it is important to note that we used the TPI method for the terraced areas only, while we continued to employ the Saulnier's method for those not terraced, assuming that these areas are not influenced by human actions and the soil thickness then depends on the natural slope conditions. In addition, since TRIGRS is specifically aimed at the analysis of shallow landslide initiation on slopes, the beaches and the urbanized areas, as well as the gullies and the cliffs, have been excluded from the thickness calculation both with the TPI and the Saulnier's method. An evaluation of the consistency of the two resulting soil thickness maps has been also made by using a 1 £ 1 m post-event DEM, which is likewise available for the study area. In detail, we compared the pre-and post-event DEMs in order to calculate the elevation difference into the landslide  areas ( Figure 14). In this respect, it is important to note that for each mapped landslide the deposit has been distinguished from the source area through the analysis of high-resolution aerial orthophotos, integrated by field surveys in the days after the event. At this point, we can assume that the elevation difference calculated within only the source areas coincides with the actual soil thickness, considering that in numerous cases the bedrock was exposed after the landslide occurrence.
Therefore, we compared the average elevation difference calculated for each of 260 source areas with the corresponding average soil thickness value estimated with the TPI and Saulnier's method at different spatial resolutions. The results, that have been expressed in terms of Mean Absolute and Root Mean Square Error (Table 3), indicate that the mean error is always higher (+ 20-30 cm) in the case of maps obtained by using the Saulnier's method. The difference between the two methods is very low just with the 8-m resolution, probably due to the excessive coarseness that does not allow to correctly reproduce the terraces' spatial pattern. However, it is important to stress that estimating the spatial variation of the soil thickness at the catchment scale can be difficult and often unreliable, since this parameter is strongly influenced by different interplaying factors, such as underlying lithology and vegetation cover (Catani et al. 2010). For this reason, an MAE ranging between 0.49 and 0.67 m (for the soil thickness map obtained with the TPI method) can, however, be considered a reasonably good result, especially in relation to the complexity of the territorial context typical of  the Cinque Terre area. In this sense, we performed a further analysis in order to describe the variation of MAE and RMSE as a function of the agricultural terrace type (i.e. still cultivated, recentlyabandoned or long-time abandoned terrace) and the presence-absence of underlying thick debris covers. The first distinction has been made for taking into account the potential influence of the vegetation cover on the estimation of the soil thickness, considering that the vegetation is much sparse in the case of still cultivated and recently abandoned terraces. According to the results (Table 4), the highest errors have been obtained for the areas characterized by the presence of thick debris covers, where the difference between estimated and real soil thickness is generally higher than 70 cm. On the contrary, such difference can be also lower than 40 cm for the terraced areas where the underlying debris cover is absent. However, it is important to note that the error calculated for the areas where the debris cover is present may be affected by the fact that in different landslide detachment zones the sliding surface did not develop at the contact with the underlying bedrock, but within the soil profile. Therefore, in such areas the real soil thickness can be underestimated, and the comparison with the value estimated with the TPI and Saulnier's method may be not completely truthful. Nonetheless, also in this analysis the mean error is always higher for the maps obtained with the Saulnier's method, whereas there does not seem to be any particular influence of the vegetation cover, considering that there is not any remarkable difference of the MAE and RMSE by varying the terrace type, both in the presence and absence of the underlying debris cover. For the estimation of the initial soil moisture conditions during the 25 October 2011 event, we employed the same approach used in the Giampilieri study case (Schilir o et al. 2015a): therefore, we considered the 1-month antecedent rainfall fallen in the study area by using the HYDRUS 1D model ( Simu nek et al. 1998), which can simulate the water flow into unsaturated porous media resulting from a rainfall event. 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 1.5 m soil profile inclined to 30 (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 2011 event. Therefore, a differentiation between terraced and not terraced areas was not made since just very few source areas are located in not terraced areas; however, it is important to stress that these areas represent just a limited portion of the entire study area compared with that covered by agricultural terraces (i.e. 24% vs. 71%, see Figure 3(a)). The initial soil moisture was assumed to be near to the residual water content value (u r ) considering the hot, dry conditions during the preceding summer months. In fact, HYDRUS 1D has been also used not only for predicting Table 4. Mean absolute (MAE) and root mean square (RMSE) errors obtained from the comparison, at different spatial resolutions, between the average elevation differences calculated for each mapped source areas and the corresponding average soil thickness values estimated with the TPI and Saulnier's method. In this case, the analysis has been performed for six different zones, identified on the basis of the agricultural terrace type (i.e. still cultivated, recently-abandoned or long-time abandoned terrace) and the presence-absence of underlying thick debris covers. u r , but also u s (saturated water content) and K s (saturated hydraulic conductivity) on the basis of the soil grain size distribution by means of the ROSETTA Lite module (Schaap et al. 2001). As regards the grain-size data, we used the information from Cevasco et al. (2014) considering the above-mentioned similarities between the Monterosso and Vernazza basins. According to the results of HYDRUS simulation, it was hypothesized the absence of an initial water table at the beginning of the 25th October 2011 event: thus, its depth (d wt ) corresponds to the bedrock-soil interface. The I ZLT parameter, that represents the long-term background rainfall rate, was assumed to be equal to the cumulative actual surface flux value, obtained by using the van Genuchten-Mualem model (van Genuchten 1980) for the simulation of the water flow. To evaluate a G parameter, that is typical of Gardner hydraulic model, we applied the conversion formula introduced by Ghezzehei et al. (2007) which defines a correspondence between Gardner and van Genuchten-Mualem models through the capillary length approach (Warrick 1995). However, in order to verify that the triggering of the landslides occurred during the 2011 event is actually related to an uprising of the perched water table (as assumed by TRIGRS), we extended the HYDRUS simulation also including the 25th October rainfall event. In this way, it was possible to describe the variation of water content (Figure 15(a)) and pore water pressure (Figure 15(b)) in response to the triggering rainfall. Starting from the initial soil moisture condition induced by the 1-month antecedent rainfall (black line in Figure 15(a)), the effect of such event results in a progressive increase of the water content within the unsaturated part of the soil, until reaching the complete saturation just in the lower portion between 13:00 and 15:00 UTC. At this point, a perched water table forms and progressively rises, up to a depth of approximately 70 cm at the end of the rainfall event. The same trend has been confirmed also for another extreme rainfall event occurred in the same area between 30 September and 4 October 1966 (see Section 2.2). However, it is important to note that in this case the initial soil moisture condition was wetter than that of the 2011 event, due to higher rainfall amounts prior to the event itself (Figure 16(a)): this point then explains how a similar final water table depth has been obtained from the simulation of the two events (Figure 16(b)), even though the 1966 rainfall event was less severe than the 2011 one (i.e. 263 vs 382 mm).

MAE (m) RMSE (m) MAE (m) RMSE (m) Zone
As regards the physical and mechanical properties of the colluvial deposits we again took into account the information from Cevasco et al. (2014). For the soil unit weight and friction angle we attributed an average value of 16.33 kN/m 3 and 31 , respectively, considering the small variability of the values among the samples. On the contrary, we decided to calibrate the cohesion starting with an initial cautionary zero value, considering the high content of coarse grained particles within the eluvial-colluvial cover. Assuming that plant roots provide additional cohesion to the soil as a function of the vegetation type (e.g. Mickovski and van Beek 2009;Burylo et al. 2011), the calibration has been iterated until a minimum number of cells resulted as unstable before the beginning of the event. On the basis of the land-use map, the cohesion then varies between 4.5 kPa (in the case of recently abandoned terraces, where the vegetation is sparse) to 6 kPa, in correspondence of forest and scrubland areas (Table 2). However, it is important to note that the chosen values lie within the range of values reported by other authors for the same material (Mondini et al., 2009;Mazzuoli and Berardi 2016). At the same time, in the case of terraces that have been abandoned for a long time and forest/scrubland areas, we decided to attribute a value of saturated hydraulic conductivity up to one order of magnitude higher than that obtained through HYDRUS 1D simulations (i.e. 4.63 £ 10 ¡6 ms ¡1 ), assuming that a much dense vegetation reduces the runoff, then favouring the seepage process. This assumption relies on the outcome of the comparative analysis between erosional processes and land use conducted by Cevasco et al. (2013a) in the Vernazza basin. In fact, the results of such analysis show that the land-use class most affected by runoff corresponds to the cultivated terraces: the frequency of such process progressively decreases passing from recently to long-time abandoned terraces, reaching its minimum in forest and scrubland areas. As a consequence, we attributed the following three different values of saturated hydraulic diffusivity: where K s is the saturated hydraulic conductivity, H is the average soil thickness (1.5 m) and S y is the specific yield. If we consider that the investigated soil can be classified as sandy loam, the specific yield has been assumed equal to 0.29, on the basis of typical values given by Johnson (1967) for each soil textural class. Finally, with regard to the rainfall input during the landslide event, different hourly rainfall maps were made available by the Regional Environmental Agency (A.R.P.A.L. -Agenzia Regionale di Protezione dell'Ambiente Ligure). These maps, which cover the entire event (i.e. 7-17 UTC), have been obtained by merging radar data and measurements from monitoring stations.

Statistical analysis of historical rainfall data
In order to evaluate the RP of landslide-triggering rainfall events, we performed a statistical analysis of maximum hourly rainfall intensity data. For this aim, we used data recorded by the Levanto monitoring station, that is the rain gauge closest to the study area for which this type of data is available for a long time span (i.e. between 1954 and 2012). 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 hourly intervals). Specifically, we used the generalized extreme value (GEV) distribution (Jenkinson 1955), which is widely used in extreme event frequency analysis rather than the Gumbel distribution, as the literature increasingly suggests that the distribution of extreme events may be more heavily tailed (Fowler and Kilsby 2003). The cumulative distribution function of the GEV distribution is as follows: where m, s and ξ are referred to as the location, scale and shape parameters, respectively. These parameters have been determined by applying the probability weighted moment (PWM) method (Hosking et al. 1985), on the basis of the maximum values of each 'cumulative rainfall' variable (i.e. 1-, 3-, 6-, 12-and 24-h rainfall) extracted, year by year, from the dataset. Finally, the inversion of the probability function yields the values of cumulated rainfall x for each of the variables and for different RPs. Then, these values have been fitted by a power law distribution in order to build the rainfall probability curves.

Back-analysis of the 25th October 2011 event
Once the input parameters of the numerical model have been defined, the 25th October 2011 event has been reproduced and the simulation results, expressed in terms of safety factor (FS), have been compared with the actual source areas of the landslides triggered during the event (Figure 17). In detail, we performed a ROC (Receiver Operating Characteristic) curve analysis in order to quantify the predictive performance of the numerical model. In fact, this type of analysis measures, for different thresholds (i.e. different values of FS), the proportion of positive values (i.e. landslide presence) that are correctly identified as such (TPR; True Positive Rate or sensitivity) together with the proportion of negative values (i.e. landslide absence) that are erroneously reported as positive (FPR; False Positive Rate or fall-out). The area under the curve (AUC) is the value that summarizes the expected performance: the larger the area, the better the accuracy of the model. According to the results (Table 5), the higher values of AUC have been obtained with the mean resolutions (2 and 4 m). It can be also noted that with the exception of the 1-m ones, the simulations which consider the spatial variation of the soil thickness based on the topographic position index (TPI) method are slightly better than those using the soil thickness resulting from the Saulnier's method ('SAU'). In particular, the best AUC result corresponds to the TPI ç 4-m simulation: for this reason, we analysed the temporal evolution of the instability with these conditions (Table 6). In this respect, calculating the relative percentage (P U ) of predicted unstable pixels, it results that only 0.6% of pixels are indicated as unstable at 7:00 UTC (beginning of the rainfall event). After a small increase in the following 4 h, the instability rapidly rose after 11:00 UTC, progressively increasing until the end of the event. Therefore, almost all the landslides would have occurred in just 5 h, since the P U passed from  6.8% to 100% between 11:00 and 16:00 UTC, and the landslide event can be considered as concluded at 16:00 UTC, given the absence of further increases in slope instability after this time. It is worth noting how the final instability peak occurs in correspondence of the maximum rainfall peak, which took place approximately between 13:00 and 14:00 UTC (Figure 4).

Evaluation of future triggering scenarios
After the back analysis of the 25th October 2011 event, we used the same model to predict the slope stability conditions of the study area by assuming different rainfall inputs, deriving from the statistical analysis of historical rainfall data described above. In this respect, the resulting rainfall probability curves ( Figure 18) emphasize how the 2011 event can be classified as an exceptional event, considering the high RP calculated for each rainfall duration (Table 7). On the basis of these results, which confirm those obtained in a preceding analysis ( With regard to the other input parameters of the model, those used for the TPI-4 m simulation of the 2011 event have been kept: in this way, by assuming the same initial conditions it was possible to perform a straight comparison with the reference landslide event. Table 8 shows the importance, in  terms of produced instability, of the 6 and 12 h rainfall amount, also for a RP equal to 25 years. In this respect, the 12-h rainfall would cause about half (50.3%) of the slope instability calculated for the entire 25th October 2011 event. In the case of 6 h rainfall, a similar level of instability develops only for events with RP equal to 75 years. On the contrary, for rainfall events having shorter duration (1 and 3 h), the number of unstable pixels is quite low even considering high RPs. However, in each rainfall condition it can be noted a certain delay between the end of the rainfall event and the instability peak: in fact, the maximum number of unstable pixels can be observed 4-8 h later, depending on the rainfall length (i.e. the less the duration, the greater the delay). This point can be explained considering the time needed for water to reach the greater soil depths, and this delay is greater for short rainfall events, considering the reduced duration of the seepage process during the event itself.

Discussions
The back-analysis performed by using TRIGRS emphasizes that the 2-m and 4-m DEMs are the most suitable for the reconstruction of the 2011 landslide event, as they yield the highest AUC values. In fact, if we compare the TPR, the FPR and the accuracy (i.e. the ratio between the total number of correct predictions and the total number of observations) obtained, for each performed simulation, by imposing a threshold FS value equal to 1 (Table 9), it can be noted that: (1) The SAU ç 1-m simulation is able to predict the 2011 landslide occurrence significantly better than the corresponding TPI simulation, given the highest TPR. However, the resulting accuracy (73.8%) is lower than those obtained from the 2-m and 4-m simulations due to the higher FPR. All these false alarms can be explained considering that the high resolution Note: The parenthesized number indicates the delay (in hours) between the end of the rainfall event and the instability peak. P 2011 represents the percentage of the maximum number of predicted unstable pixels compared to those obtained in the back-analysis of the 25th October 2011 event. implies an extreme heterogeneity of the local slope which, in turn, determines an undue spatial variation of the soil thickness, which then influences the model output ( Figure 19); (2) The TPI ç 2-m and 4-m simulations show higher accuracy values compared to the corresponding SAU simulations because the number of false alarms is lower, even though the TPR values are similar. As regards the TPI simulations, it can be also noted that the 2-m simulation correctly identifies the landslide source areas slightly better than the 4-m one, and vice versa. However, as stated above, the best result in terms of AUC has been obtained with the TPI ç 4-m simulation: this means that this spatial resolution, coupled with the TPI method, allows to efficiently reproduce the shallow slope failures by keeping the spatial information concerning the pattern of the agricultural terraces ( Figure 19). In this respect, for this spatial resolution we also calculated the TPR, FPR and accuracy by distinguishing between terraced and not terraced areas (Table 10). As it can be noticed, the use of TPI method allows to increase the accuracy in the terraced areas of more than 5 percentage points due to the reduction of FP. With regard to the not-terraced areas, although the TPR is very low, the accuracy is extremely high since just very few source areas are located in such areas; (3) The TPI and SAU ç 8-m simulations tend to underestimate the landslide occurrence, according to the low TPR. This result can be explained considering that the coarseness of the spatial resolution makes it difficult to properly reproduce either the shape of the landslide source areas or the real pattern of the agricultural terraces, unlike the highest resolutions ( Figure 19).
At this point, considering that the best results have been obtained with the 2-m and 4-m resolutions, we performed a further analysis to quantify the number of simulated landslides to changing the percentage of correctly predicted source area. In other words, assuming that a landslide can be considered 'triggered' if at least 50%, 60%,…,100% of pixels within the source area polygon results as 'unstable', we calculated the percentage of landslides properly computed with the two spatial resolutions and considering the two different methods for the estimation of the soil thickness ( Figure 20). The results show that the 2-m simulations identify a slightly higher number of landslides compared to the 4-m simulations, but only for low thresholds of correctly predicted source area (50%-60%). In fact, a rapid reduction of the number of correctly simulated landslides can be noted with increasing the percentage of correctly predicted source area. On the contrary, the reduction is much less pronounced in the case of the 4-m simulations, which show a relatively high percentage of correctly simulated landslides even for high thresholds. In addition, the comparison between TPI and SAU simulations indicates that for the 2-m simulations the best results have been obtained with the Saulnier's method, unlike the 4-m simulations. This point confirms how the TPI method is the most  suitable technique for the estimation of the spatial variation of the soil thickness, at least considering the topographic peculiarities of the Cinque Terre area. As far as the evaluation of future triggering scenarios is concerned, the results show that the highest amounts of unstable pixels have been obtained for the longest rainfall events (i.e. 6-12 h). In fact, considering the relatively low hydraulic conductivity of the soil together with the unsaturated conditions, the partial absorption of water infiltrating at the surface results in damping and smoothing of the rainfall input, which cause the triggering of fairly few cells in response to short rainfall events. In this respect, if we compare these results with those obtained for the case study occurred in Giampilieri in 2009 (see Section 3.2), it can be noted thatin that instancethe wetter initial soil moisture conditions resulted in a high number of unstable cells even in response to short rainfall events (Schilir o et al. 2015a). From this point of view those initial soil moisture conditions have certainly enhanced the triggering of shallow landslides, also considering that the 2009 rainfall event was far from exceptional, at least if compared with that occurred in the Monterosso catchment in 2011. This finding then emphasizes that, within the definition of realistic triggering scenarios for rainfallinduced shallow landslides, it is essential a proper evaluation not only of the expected rainfall event but also of the soil hydraulic properties and the initial moisture conditions.

Conclusions
In this study, we used a physically based model for the prediction of rainfall-induced shallow landslides in a coastal catchment located in the Cinque Terre area (Eastern Liguria, Northern Italy). As a first step, we performed a back-analysis of a landslide event occurred in the study area on 25 October 2011 by using four different spatial resolutions of the digital terrain model and two different methods for the evaluation of the spatial variation of the soil thickness. The results indicate that the numerical simulations reconstruct quite well the space-time evolution of the real event and, more specifically, that the best results have been obtained with the mean resolutions (2 and 4 meters) and considering the TPI method as the technique for the spatial definition of the soil thickness. In particular, the best result in terms of AUC has been obtained with the TPI ç 4-m simulation, which confirms that the spatial information concerning the pattern of the agricultural terraces is maintained even at this resolution. At a broader level, it can be also asserted that the TPI method is slightly better than the SAU method since it allows to reduce the false alarms (FP) without reducing the number of correctly predicted landslide pixels (TP). With regard to the evaluation of the RP of the 2011 event, the results confirm the exceptionality of the event itself. However, the analysis of different triggering scenarios shows that long (6-12 h) rainfall events even with relatively low RP can cause a substantial instability level: this latter point would therefore explain the landslide/flood events that have occurred in the past in the same area. The results of the numerical simulations also emphasize the importance of the specific local conditions (such as, in this case, the widespread presence of agricultural terraces) in influencing the occurrence of shallow landslides. For the parameterization of the numerical model we have then considered that different land-use types affect not only the spatial distribution of the soil thickness, but also the geotechnical parameters and the hydraulic conditions of the soil: in our opinion, these aspects cannot be disregarded in this type of analysis and, if taken into account, allow to obtain reliable results. Although future improvements will definitely concern a better evaluation of the input parameters, from a methodological point of view we believe that the proposed approach enhances the potential of physically based models as tools for hazard assessment and risk management.