Rainfall-Induced Shallow Landslide Detachment, Transit and Runout Susceptibility Mapping by Integrating Machine Learning Techniques and GIS-Based Approaches

Abstract

Rainfall-induced shallow landslides represent a serious threat in hilly and mountain areas around the world. The mountainous landscape of the Cinque Terre (eastern Liguria, Italy) is increasingly popular for both Italian and foreign tourists, most of which visit this outstanding terraced coastal landscape to enjoy a beach holiday and to practice hiking. However, this area is characterized by a high level of landslide hazard due to intense rainfalls that periodically affect its rugged and steep territory. One of the most severe events occurred on 25 October 2011, causing several fatalities and damage for millions of euros. To adequately address the issues related to shallow landslide risk, it is essential to develop landslide susceptibility models as reliable as possible. Regrettably, most of the current land-use and urban planning approaches only consider the susceptibility to landslide detachment, neglecting transit and runout processes. In this study, the adoption of a combined approach allowed to estimate shallow landslide susceptibility to both detachment and potential runout. At first, landslide triggering susceptibility was assessed using Machine Learning techniques and applying the Ensemble approach. Nine predisposing factors were chosen, while a database of about 300 rainfall-induced shallow landslides was used as input. Then, a Geographical Information System (GIS)-based procedure was applied to estimate the potential landslide runout using the “reach angle” method. Information from such analyses was combined to obtain a susceptibility map describing detachment, transit, and runout. The obtained susceptibility map will be helpful for land planning, as well as for decision makers and stakeholders, to predict areas where rainfall-induced shallow landslides are likely to occur in the future and to identify areas where hazard mitigation measures are needed.

1. Introduction

Rainfall-triggered landslides are notoriously dangerous phenomena capable of causing serious damage and a notable death toll as well as important economic losses all over the world. Shallow landslides usually occur on steep soil-mantled slopes [1,2,3,4,5,6] as a consequence of heavy or intense rainfalls and often may evolve into potentially catastrophic flow-like movements. Slope failures are very common in hilly and mountainous areas due to the combination of several controlling factors such as slope steepness [7,8,9,10], intense or prolonged rainfalls [2,11,12,13,14], land-use changes [15,16,17,18] and wildfires [19,20]. The consequences of shallow landslides and flow-like phenomena are typically more dangerous in inhabited centers located at the foot of slopes, where the hydrographic network is well developed.

The spatial and temporal hazard posed by flow-like movements is determined by both source-related features (e.g., location and volume) and subsequent runout dynamics (e.g., travelled paths and distances) [21]. Therefore, in addition to the identification of potential landslide source areas, it is crucial to investigate landslide runout (i.e., travel distance) and incorporate it into landslide hazard or susceptibility assessments. For this reason, landslide susceptibility zoning should include information about landslide source areas as well as downslope areas that might be affected by mobilized materials. In recent years, a meaningful improvement in landslide detachment susceptibility evaluation has been gained through robust scientific advances, especially by using statistical approaches [22,23,24,25,26]. On the other hand, runout of flow-like landslides has been studied by using simple flow direction algorithms [27,28,29], empirical-statistical approaches [29,30,31,32], analytical methods [33] and dynamic methods [34,35,36,37]. The integration of landslide detachment and runout dynamics has been accepted by many authors as an appropriate and complete methodology for Landslide Susceptibility (LS) assessment [9,38,39,40,41,42,43,44,45]. However, despite its importance, runout evaluation is not as widespread in literature as the analysis of the potential for landslide detachment [44]. Nowadays, only a limited number of studies include both Landslide Detachment Susceptibility (LDS) assessment and Landslide Runout Susceptibility (LRS) assessment. Most of these studies focused on either landslide susceptibility mapping or mobility assessment by considering them separately. The coupled prediction of the potential for landslide detachment and runout remains a challenge for researchers [46].

In this study, a combined methodology to assess susceptibility of shallow landslide failure and runout has been developed. Such a procedure is based on the integration between LDS and LRS assessment using, respectively, Machine Learning (ML) statistical-probabilistic models and Geographical Information System (GIS)-based tools. This methodology has been applied to the Cinque Terre National Park, Liguria, north-western Italy; hereafter CTNP), where risk posed by flow-like movements is very serious. ML techniques and GIS-based tools were adopted based on their reliable predictive performance and capability to perform large-scale analysis, respectively. Moreover, GIS-based tools are easy to use in any GIS software achieving high performance. Shallow landslides are a considerable geo-hazard within the CTNP owing to severe recurrent rainfalls along with the peculiar geomorphological and land-use setting [47,48]. On the other hand, CTNP is increasingly popular: every year millions of tourists visit this wonderful terraced coastal landscape to enjoy a beach holiday and to practice hiking activities along a dense trail network. Landslide susceptibility evaluations are significant tools in current land-use planning. Unfortunately, most of the available landslide zoning maps only examine susceptibility to landslide detachment. Usually, hazard zoning maps do not take into account transit and invasion analyses. These shortcomings have a great impact on final hazard models and on urban planning. Susceptibility mapping including detachment and runout can represent an important tool for urban planning administrations and policymaker. Furthermore, a complete analysis of these geological hazards is fundamental in risk assessment and management, since the landslides runout mainly considers the areas located in the lowland zones. In the detachment susceptibility analyses, these areas are always underestimated and classified with low to extremely low landslide susceptibility.

2. Study Area

The name Cinque Terre (literally Five Lands) originates from the five historic coastal hamlets of Monterosso al Mare, Vernazza, Corniglia, Manarola and Riomaggiore (Figure 1). These hamlets are located in the easternmost Liguria region (north-western Italy), along a coastal stretch which represents a worldwide known tourist attraction. Because of its environmental, cultural, and historical heritage, the Cinque Terre area was listed as a World Heritage Site by UNESCO in 1997 and was declared National Park in 1999.

2.1. Geological and Geomorphological Setting

From a geological point of view, CTNP belongs to a NW-SE oriented portion of the Northern Apennine, an orogenic chain formed during the Tertiary [50]. This sector is made up of five tectonics units belonging to the Tuscan, Sub-Ligurian and Ligurian domains [51]. These units are involved in a wide overturned antiform SW-verging fold that is bounded to NE by a major normal fault (La Spezia fault). Structural features correlated to the Apennine orogeny, and especially to the Plio-Quaternary tectonic uplift [52,53], strongly control the current topography. From top to bottom, the five overlapping tectonics units include: Gottero Unit, Ottone Unit, Marra Unit, Canetolo Unit and Tuscan Nappe [50]. The Gottero Unit, which crops out in the westernmost part of the National Park, is composed of ophiolitic and sedimentary formations; the first ones mainly include serpentinites and gabbros, while the second ones consist of turbidite sequence (Mt. Gottero Sandstones Fm., Late Cretaceous). The Ottone, Canetolo and Marra units crop out along thin and NW-SE-trending strips of territory that occupy the western and the central sectors of the CTNP. The Ottone Unit is prevalently made up of pelitic rocks (i.e., Monte Veri Complex Fm.), being composed by dark grey marls with calcareous and arenaceous intercalations. The Marra and Canetolo Units mainly consist of argillites and claystones with intercalations of calcarenites, siltstones, sandstones and silty-sandy turbiditic rocks (e.g., Canetolo shales and limestones Fm., Torrente Pignone marls Fm.). The predominant geological formation is a flysch complex consisting of sandstones and clayey siltstones known as Macigno Fm. (Tuscan Nappe, Upper Oligocene). This formation crops out along almost the entire coast between Monterosso al Mare and the eastern border of the CTNP. Finally, the Tuscan Nappe also includes small outcrops of limestones, only visible in the proximity of the easternmost border of the CTNP (Maiolica and Scaglia Toscana limestones).

The geomorphological setting of CTNP is characterized by small coastal catchments, ranging from 1–3 km² (e.g., Fegina, Pastanelli, Groppo and Riomaggiore stream basins) to around 6 km² (i.e., Vernazza stream basin), separated from the inner Vara Valley by the main watershed that runs parallel and very close to the coastline (Figure 1b). Such basins are characterized by steep slopes which are typically mantled by thin soil covers (on average 1–2 m thick) made up of eluvial-colluvial deposits [47]. Over the past centuries, these slope materials have been extensively human–reworked to obtain cultivable areas. CTNP is frequently affected by rainfall-induced shallow landslides (e.g., debris slides and debris avalanches), often evolving into flows, with the failure surfaces usually set at the soil–bedrock interface [6,54]. Landforms and deposits related to mass movements are extremely widespread, since gravity is one of the main geomorphologic agents [55]. The high steepness of slopes represents an important predisposing factor for landslides along with favoring areal and channelled erosion due to running water [55].

According to the Unified Soil Classification System (U.S.C.S.), most eluvial-colluvial soils consist of silty gravels, clayey silty gravels and clayey gravels [6]. The thickness of slope cover depends on land-use and geological settings. Very thin soil covers (<0.5 m) occur on wooded slopes characterized by ophiolitic bedrock [48]. Slightly greater soil thicknesses (up to 1 m) can be found on wooded slopes with sedimentary bedrock [48]. In terraced areas, due to the anthropogenic reworking of detrital covers, soil thickness can reach several meters [6]. The presence of high elevations in proximity to the coastline does not permit a well-developed hydrographic network. Catchments are generally characterized by narrow and deep-cut valleys with short streams showing an ephemeral hydrological regime, whereas the coast is prevalently shaped by plunging rocky cliffs [48].

The morphological features of the CTNP deeply control the local climate setting. The combined action of southerly exposition and orographic effect enhances the features of the Mediterranean climate, here characterized by particularly warm dry summers and temperate winters. However, due to the humidifying action of the sea on the air masses coming from the southern quadrants along with the effect of the Alpine–Apennine chain, relatively abundant rainfalls affect the autumn and winter seasons [48]. The mean annual temperature ranges between 14 and 15 °C; the hottest months are July and August, with mean temperatures of 22–24 °C; the coldest month is January, with a mean temperature of 7–8 °C. On average, the Cinque Terre area receives 900–1100 mm of rainfall annually [48,54]. Extraordinary rainfall events can occur in late summer and early autumn, especially in the ongoing climate change context [48,55,56]. During the rainfall event of 25 October 2011, a cumulative rainfall of 382 mm and maximum rainfall intensities of about 90 mm/h, 195 mm/3 h and 350 mm/6 h were recorded at Monterosso rain gauge [6,57].

2.2. Agricultural Terraces

Most of the Liguria region is characterized by a rugged morphology and by a lack of wide flat areas suitable for cultivation, both along the coast and in the hinterland. To obtain well-suited areas for agricultural harvesting, slope soil covers have been extensively reworked over the centuries to build terraces sustained by dry-stone walls. In the study area, this impressive work of slope modification resulted in an exceptional example of human integration with the natural landscape [58]. Indeed, wide portions of the CTNP area, from the shoreline up to 400–500 m a.s.l., have been modified through slope terracing, producing an unusual man-made coastal landscape which is perfectly integrated with the local geomorphological setting (Figure 2). Terraced slopes cover an area of about 15 km², corresponding to about 41% of the entire CTNP area; however, approximately 67% of the terraced areas have been abandoned [55]. Terranova [59] estimated that the maximum linear extension of dry-stone walls could be approximately 6000 km. It is widely documented that agricultural terraces, if properly maintained, have a positive role on soil conservation by contributing to reduce runoff velocity and soil erosion [61 and references therein]. On the other hand, a lot of studies demonstrated that when these traditional systems are abandoned, land degradation issues (i.e., gully erosion, terrace failure, mass movement, piping, hydrological connectivity) arise [60,61,62]. The effects of terrace abandonment may even result in increased risk conditions when slopes are located in proximity to urbanized areas. In the study area, abandoned terraced slopes are particularly susceptible to shallow landsliding, especially due to intense rainfall [6,18,48]. During high-intensity rainfall, materials mobilized by mass transport and landslide phenomena have a considerable impact on the solid discharge, flow and energy of streams. These dynamics, which are quite common in the entire Liguria region, often cause flow-like movements, debris floods and flash floods, which are increasingly affecting, in particular, the coastal settlements [11,48,63,64]. Within the CTNP, these phenomena represent a serious threat to human settlements, inhabitants and trail users, as dramatically observed after the intense rainstorm that hit the Monterosso and Vernazza areas on 25 October 2011 [48,65,66,67,68], which triggered hundreds of shallow landslides as well as destructive debris floods [6,47,56].

3. Materials and Methods

To achieve the objective of this study a three-step methodology was applied (Figure 3). The first step consisted of modelling and mapping the Landslide Detachment Susceptibility (LDS). After selecting input data (i.e., landslide inventory and predisposing factors) and identifying any collinearity problems, ensembled ML models were applied. The ML techniques used in this work are: (i) Artificial Neural Network (ANN) [69], (ii) Generalized Boosting Model (GBM) [70] and, (iii) Maximum Entropy (MaxEnt) [71].

The second step was aimed at modelling Landslide Runout Susceptibility (LRS). The “reach angle” approach [72] was adopted and implemented on a Geographical Information System (GIS), allowing to identify Maximum Invasion (MI) points of landslides and to establish their degree of susceptibility. Eventually, as a third step, information coming from the previous analyses were used to generate the Landslide Detachment, Transit and Runout Susceptibility (LDTRS) map.

3.1. Landslides Inventory

CTNP has been historically affected by recurrent landslide events, mostly due to intense rainfalls [64]. A landslide inventory map of CTNP was recently compiled by Raso et al. [73] by means of field surveys and visual interpretation of satellite images and aerial photos. This map was produced at 1:45,000 scale and it includes landslides with a size greater than 25 m². This inventory consists of more than 450 landslides classified according to Cruden and Varnes [74] and Hungr [75]. From this inventory, a total of 281 rainfall-induced shallow landslides, mostly debris slides, debris avalanches and debris flows were selected (Figure 4).

To perform the LDS analysis, as the algorithm works with punctual features, the selected landslides were represented by points. The highest elevation point of each landslide was chosen as representative, according to the approach adopted by the Italian I.F.F.I. project (Inventario dei Fenomeni Franosi in Italia—Landslide Inventory in Italy), where each landslide could be statistically identified by mean of its P.I.F.F. (Punto Identificativo del Fenomeno Franoso—Landslide Phenomenon Identification Point). By convention, P.I.F.F.s correspond to the highest elevation points on the landslide crown [76]; hence, they can be regarded as landslides source points.

3.2. Predisposing Factors

Based on the previous knowledge of the study area, the following Predisposing Factors (PFs) were selected for the LDS assessment: (i) slope angle, (ii) lope aspect, (iii) planform curvature, (iv) profile curvature, (v) distance from roads, (vi) distance to streams, (vii) degree of abandonment of agricultural terraces, (viii) land-use and, (ix) geo-lithological features.

All morphometric variables were derived from a 10 × 10 m Digital Terrain Model (DTM) produced by the Liguria Region Administration in 2017 (data downloaded from https://geoportal.regione.liguria.it/ accessed on 11 February 2021). It should be noted that this analysis has been focused only on environmental variables, without taking into consideration triggering factors (e.g., intense rainfalls) [48]. Slope angle is one of the most significant PFs in landslide occurrence and its values are notoriously inversely connected to slope stability and the soil cover thickness [5,77]. Aspect plays a fundamental role in controlling soil moisture along slopes through solar radiation and rainfall [78] as well as influencing many factors regulating slope stability such as hydrological processes, weathering and vegetation cover [79]. Curvature provides information on slope shape: planform curvature relates to the convergence and divergence of flow across a surface, while profile curvature describes the direction of maximum slope, indicating the acceleration and/or deceleration of the flow across the surface. Distance from roads and rivers represents a significant factor controlling the occurrence of slope mass movements [80]. Specifically, roads represent morphological discontinuities of the slope profile while streams are preferential ways where landslides can increase their mobility. As different studies demonstrated, the degree of abandonment of agricultural terraced slopes represents an important landslide susceptibility factor in the study area, strongly influencing the magnitude of rainfall-induced mass movements [6,54,81,82,83,84,85]. In this study, special care has been paid to the recognition and mapping of both agricultural terraces and non-terraced areas. In particular, using data from previous studies [12,55], four classes were distinguished: 1. cultivated terraced areas; 2. abandoned terraces with poor vegetation cover (initial state of abandonment); 3. abandoned terraces with dense vegetation cover (advanced state of abandonment); 4. non-terraced areas (including either urban areas, outcropping rocks and woods). Land-use, because of its temporal and spatial variations, can greatly affect landslide occurrence as PFs [15]. In particular, it controls landslide hazard and influences the development of the spatial pattern of elements at risk [86,87]. The land-use map has been obtained from the Corine Land Cover project (CLC), released in 2018 by the European Environment Agency (available from https://geoportal.regione.liguria.it/progetti/uso-del-suolo-e-tipi-forestali.html accessed on 11 February 2021). The geo-lithological map (1:20,000 scale) was re-elaborated by using bibliographic information [50] and archive data available on the Liguria Region Administration geoportal (https://geoportal.regione.liguria.it/catalogo/mappe.html accessed on 11 February 2021).

3.3. Landslide Detachment Susceptibility (LDS) Assessment

To assess the detachment susceptibility for shallow landslides, ML approaches were applied. Specifically, Artificial Neural Network (ANN) [88,89], Generalized Boosting Model (GBM) [90] and Maximum Entropy (MaxEnt) [71,91] techniques were used because of their reliable predictive performance [24,85,92,93,94,95]. The goodness of the predictive accuracy of the susceptibility models was estimated using a K-Fold cross-validation approach [22]. This procedure consists in randomly partitioning the total dataset into “k” parts of an equal number. Afterwards, at each step, one dataset is used for model testing while the combined remaining “k-1” datasets are used for model training. This approach avoids both overfitting problems and asymmetrical sampling (i.e., bias) of the training dataset, which is typical of methods involving the division of the dataset into two parts (i.e., training and testing dataset). The whole procedure was iteratively repeated for each model and, lastly, the average predictive accuracy was specified [96]. Precisely, each algorithm was implemented 50 times, splitting 80% for training and the residual part for testing. To detect multicollinearity problems, a Variance Inflation Factor (VIF) was applied. VIF measures how much the behavior (i.e., variance) of an independent variable is influenced, or inflated, by its interaction/correlation with the other independent variables: a high VIF indicates that the associated independent variables are highly collinear with the other variables in the model. There is no well-defined and recognized threshold in the literature; the general rule of thumb is that VIFs exceeding 4 necessitate further investigation, while VIFs exceeding 5 or 10 show serious multicollinearity requiring correction [97]. The predictive performance was evaluated with spatial cross-validation estimates of the Area Under the Receiver Operating Curve (AUROC), which is independent of a specific decision threshold [98]; such a curve is one of the most common statistics to assess model performance. The AUROC curve plots all the possible true positive rates (i.e., sensitivity) to the corresponding false-positive rates (1—specificity). AUROC values close to 50% indicate no discrimination, while an AUROC close to 100% indicates perfect discrimination between binary prediction classes. After the completion of the stand-alone methods, the Ensemble Modelling (EM) technique was run (Figure 3). EM allows one to create multiple models to predict an outcome, either using many different modelling algorithms or using different training datasets. The main reason for using EM is to reduce the generalization error of prediction. Since the stand-alone models are independent, the prediction error of the model decreases when the EM approach is used. Even though EM has multiple base models, it operates as a single model. Both modelling and statistical analysis were conducted entirely in R [99], principally using the “Biomod2” package [100]. Four ensemble techniques, namely mean of probabilities (PM), median of probabilities (PME), binary mean or committee averaging (CA) and weighted mean of probabilities (PMW), were used. For the highest performance value, the weighted mean was chosen (AUCWMean = 0.9), whereas to estimate the predictive capability of each method, the AUROC curve was adopted. More precisely, for AUROC values lower than 0.5, outcomes are poor, values in the range 0.5–0.7 are sufficient, fair in the range 0.7–0.9, and excellent for values greater than 0.9 [101,102]. To identify the potential landslide source areas, the continuous LDS map was classified into six susceptibility classes, namely from “very low” to “very high” susceptibility, using the Natural Breaks classification methods. The Natural Breaks criterion sets the limits between two classes in correspondence with discontinuities or “jumps” in the frequency distribution [103]. Such a method aims to define the best data distribution in significant classes, taking advantage of the intervals naturally present in the distributions of data, trying to minimize the variance within each class and maximizing that one between the classes themselves. Furthermore, the intersections between landslide points and the LDS map have been carried out as validation procedure. Such a procedure was performed to find out the percentage of the landslide areas and landslide source points, falling into different susceptibility classes, identified by the model. In order to obtain an outcome comparable with the existing landslide hazard zoning maps provided by the local authorities, six landslide susceptibility classes have been chosen.

Another significant consideration is related to the importance of variables, which can support the understanding of the environmental predictors’ contribution to the model outputs. The adopted criterion is based on the randomization of a single variable of input data. Therefore, a model prediction with a ‘randomized’ dataset was made. Subsequently, a simple correlation was computed (Pearson’s correlation by default) between reference predictions and the ‘randomized’ one. The returned score is 1-cor (pred. reference, pred. randomized). A good correlation score between two predictors means that one of the two variables has a low influence and thus it is considered not important by the model. On the contrary, a low correlation score indicates that both variables have a high influence on the model. It should be noted that this method does not account for interplay between variables but every factor has to be considered individually and should be evaluated more as an information tool for each model independently [100]. In this work, to overcome this issue, the variable importance was divided by the highest variable importance value so that values ranged between 0 and 100.

3.4. Landslide Runout Susceptibility (LDR) Assessment

The first experiences on runout assessment were made by Heim [104] and Terzaghi [105]. The main objective of these studies was to understand the characteristics of landslides that occurred in the past to predict the possible future invasion areas, therefore preventing economic losses and fatalities [106]. Calculating the runout of mass movements and understanding their rapidity are basic elements of any landslide hazard or risk assessment [106]. In this work, LRS assessment was performed adopting the geometrical method of the “reach angle” [107], also known as “fahrböschung” [104], which is one of the most used empirical approaches to evaluate landslide mobility. The main parameters used in this approach are shown in Figure 5. Runout distance (L) is defined as the horizontal projection of the line (straight line) that connects the highest point of the landslide source area with the landslide deposit toe; slope height (H) represents the difference in elevation between the landslide crown and the farthest edge of the displaced mass, while θ is the slope angle before the failure occurrence [30]. The model assumes that the landslide will arrest at the intersection point of the straight line with the topography. The intersection of the straight line with the topographic surface provides both H and L [9], from which the reach angle (α) can be determined. Corominas [30] demonstrated that the reach angle depends on the landslide volume, the type of movement and topographic constraints along the path. Moreover, runout behavior of flow-like landslides can be controlled by several factors, such as land-use, soil properties, local morphology and water content. Therefore, it is not easy to determine the runout path and the travel distance for future debris flows [27].

In the first phase of LRS assessment, values of H and L were carefully measured, in an open-source GIS environment. This was carried out both for landslides triggered along or upstream of channels and for failures on open/planar slopes. Subsequently, the reach angle was calculated through the following formula Equation (1):

$$\alpha = \arctan \frac{H}{L}$$

(1)

The reference values of the reach angle for the two different slope types (i.e., channelled and open/planar) were obtained through the definition of the modal values. During the second phase, using the LDS map as a base map, different topographic sections were drawn. These sections identify the hypothetical/potential path that a detached landslide could travel. The sections were drawn along the main creeks, in case of channelled slopes, and along the maximum slope line, in case of open/planar ones. Topographic sections were automatically obtained using the Profile tool plug-in [109] in open-source QGIS. In the third phase, MI points were identified for different potential detachment points positioned in the upper part of the slope (Figure 6a). The straight-line equation Equation (2) was used to define the MI point corresponding to a certain potential detachment point. Such a straight-line is described by the following equation:

$$y=x \mathrm{tg}\alpha +c$$

(2)

where:

α = reach angle;

x = progressive distance;

c = highest detachment point altitude.

Starting from the highest elevations, different detachment points (C_n) were considered along the topographic profile. Each point was positioned within a certain susceptibility class derived by the LDS modeling. Subsequently, every point of the straight-line was iteratively compared with the local topography obtained through the Digital Elevation Model (DEM). In this way, the MI point (Y’_n) is identified by the intersection between the straight line and the local topographic profile (Figure 6b). The susceptibility assigned to the MI point is the same as its potential detachment point.

3.5. Landslide Detachment, Transit and Runout Susceptibility (LDTRS) Assessment

To produce the LDTRS map, both LDS and LRS information were used. Specifically, the implementation of the LDTRS map was carried out through areas enveloping procedure aimed at (i) harmonizing areas with different LDS and (ii) including information derived from the LRS assessment (i.e., areas showing the same susceptibility degree to both detachment and runout were enclosed in a single area). To ensure higher safety conditions, a conservative approach was adopted during the areas enveloping procedure. Hence, priority was always given to the highest susceptibility classes (Figure 7). In this way, the areas located downstream of a portion with a higher susceptibility are included in these ones, even if they are characterized by a lower susceptibility.

Figure 8 summarizes the different operations carried out during the areas enveloping procedure. LDTRS areas were drawn in order of decreasing susceptibility class (Figure 8b–f).

Lastly, an attempt was made to obtain an indication on the goodness and effectiveness of the LDTRS model. The test area of the Vernazza catchment, where data from the catastrophic event of 25 October 2011, were available [12], was selected. The adopted testing methodology consists of a quantitative comparison of the deposition areas of shallow landslides (mostly debris slides and debris avalanches), detected through field surveys and aerial imagery analysis, with the respective degree of susceptibility resulting from the output LDS and LDTRS maps. As debris flows propagated into the drainage network, they were not included in this analysis. In fact, in this case, the deposits were swept away along the channels and uncertainties occurred in identifying the boundary between debris flows and debris flood deposits [14].

4. Results

4.1. Landslide Detachment Susceptibility (LDS) Assessment

In landslide susceptibility modelling, the predictive capability of the predisposing factors should be considered and factors with non-predictive value should not be included in the models, improving the performance of the prediction. In this study, VIF was examined to estimate the possible correlations between each variable. There are no well-defined and recognized VIF thresholds reported in literature (see Section 3.3). The general rule of thumb is that VIFs exceeding 4 need more analysis, while VIFs exceeding 5 or 10 denote serious multicollinearity issue. In this study, VIF values of each predisposing factors are smaller than 4. The highest VIF is 1.94, corresponding to slope angle, while the smallest VIF is 1.02, which is associated with slope aspect. No multicollinearity problems were found among the predisposing variables considered in the analyses and the outcomes have been summarized in

Table 1. Multicollinearity diagnosis indices for the selected variables and their importance for the Weighted mean (Wmean) ensemble model.

Predisposing Factors	Variance Inflation Factor	Variable Importance (Wmean)	Values between 0 and 100 (%)
Slope angle	1.94	0.22	58
Slope aspect	1.02	0.24	63
Profile curvature	1.35	0.04	10
Planform curvature	1.54	0.12	32
Degree of abandonment of agricultural terraces	1.39	0.38	100
Distance to road	1.18	0.02	5
Distance to river	1.28	0.01	3
Geology	1.36	0.04	10
Land Use	1.14	0.06	16

.

Another meaningful outcome, which allowed to understand the contribution of PFs to the model outputs, is related to the significance of the selected variables. The results of the analysis show that slope angle, slope aspect and the degree of abandonment of agricultural terraces are the most important factors affecting landslide models (Table 1).

The most important PFs for LDS assessment are related to agricultural terraced slopes that were abandoned and covered by Mediterranean scrub, unlike the areas where agricultural terraces are still cultivated or they have been recovered in recent years. Other relevant environmental variables for LDS assessment are some morphological factors such as slope angle and slope aspect (Table 1).

The LDS map that resulted from the calculation of the weighted mean probability of the different models was reclassified into six classes describing the different susceptibility to landslide detachment (Figure 9a).

As can be seen from Figure 10, the distribution of landslides is characterized by an increasing trend, with the highest number of landslides placed in the high and very high susceptibility classes (ca. 63%). On the contrary, the areal distribution shows a decreasing trend, with the majority of areas characterized by low and very low susceptibility values (ca. 58%) and a small area characterized by higher values (ca. 14%). The opposite direction of the two trends reveals that most of the landslides fall within the highest susceptibility classes, although the range recognized as susceptible has a limited extension.

4.2. Landslide Runout Susceptibility (LRS) Assessment

As explained in Section 3, the first step of LRS assessment was to determine the reach angle of the inventoried landslides. As reported in Figure 11, reach angle values calculated for channelled landslides were generally lower than those related to landslides that occurred on open/planar slopes. In particular, the most frequent reach angle values for channelled and open/planar slopes landslides are in the range 26°–30° and 36°–40°, respectively. The choice of the reference values was made through the definition of the modal value, namely the higher frequency identified in the two respective classes, adopting a conservative approach. Therefore, values equal to 28° and 36° for landslides potentially occurring on channelled and open/planar slopes were chosen, respectively.

Topographical sections were drawn as reported in Section 3.4. In this way, 453 topographic sections were drawn (Figure 9b). This procedure allowed to identify 1853 MI points; subsequently, an LRS class has been assigned to each of them (see Section 3.4). The number of identified MI points is higher than the topographic sections drawn since, as explained in Section 3.4, along each section it was possible to identify up to five MI points.

4.3. Landslide Detachment, Transit and Runout Susceptibility (LDTRS) Assessment

The final product is a detailed LDTRS map resulting from the integration of the procedures used for LDS and LRS assessment (Figure 9a,b). In this way, besides susceptibility to shallow landslides detachment, both transit and runout susceptibility were taken into account (Figure 9c,f). The analysis of LDTRS map (Figure 9c) highlights differences in the spatial distribution of susceptibility. Generally, terraced coastal slopes show high and very high susceptibility. Regarding the inland portions of the territory, high susceptibility values can be observed within the Vernazza catchment and in the areas behind the Monterosso al Mare hamlet. Interestingly, in the Vernazza catchment, there is a marked difference in susceptibility values between the western and the eastern side of the basin. In this case, the greatest landslide occurrence and susceptibility values were found on the right side of the catchment. On the other hand, the lowest landslide susceptibility classes are located at a higher altitude. Furthermore, most of the stream channels show high susceptibility values, highlighting their crucial role during debris flows runout.

Eventually, indication on the goodness of the LDTRS model was obtained for the test area of Vernazza through the comparison between deposition areas of shallow landslides that occurred on 25 October 2011, and respective degree of susceptibility in the LDS and LDTRS maps. The obtained results are summarized in Figure 12. It can be noticed that the percentage of deposits falling in areas with high and very high susceptibility increases from about 48% for LDS to about 82% for LDTRS. At the same time, the percentage of deposits falling in areas with very low to medium susceptibility decreases from about 52% for LDS to about 18% for LDTRS.

5. Discussion

Landslide susceptibility maps are important tools in contemporary land-use planning. Regrettably, most of the current landslide zoning maps only consider susceptibility to landslide detachment, neglecting landslide transit and invasion. Consequences of incorrect urban and land-use planning can be disastrous, especially for areas where flow-like movements are likely to propagate [110]. In fact, these phenomena are frequently triggered in hilly and mountainous areas and, after travelling long distances, can impact buildings and infrastructure located in lowland areas, threatening people’s safety. In Italy, the well-known landslide disasters of Sarno (Campania Region, southern Italy) [10,77,111] and Giampilieri (Sicilia Region, southern Italy) [112,113], which occurred in 1998 and 2009, respectively, are impressive examples of loss of lives and damage experienced by urbanized areas and infrastructure due to long-distance flow-like movement propagation (up to some km from detachment). As regards the Liguria region, the recent cases of Cinque Terre in 2011 [6,12] and of Leivi (central-eastern Liguria) in 2014 [114,115,116] are still impressed on local people’s memory since both caused fatalities and severe damage. These natural disasters undoubtedly suggest the importance of improving existing urban planning tools, namely landslide susceptibility maps, including runout analyses. In this work, a combined approach was used to obtain a thematic map considering both aspects. Outcomes of the first part of this work have been focused on the LDS assessment. Such results revealed that the main causes of slope instability are: (i) abandonment of agricultural terraces (Wmean parameters importance = 0.38), (ii) slope aspect (0.24) and (iii) slope angle (0.22). The high LDS of terraced areas can be chiefly interpreted as the outcomes of two conditions. Firstly, the flat portion of agricultural terraces hampers runoff avoiding soil losses, and, secondly, it increases the water infiltration [6,62]. During rainfall, the combination of these two hydrological functions leads to the formation of a groundwater table at the interface between layers (i.e., slope deposit and underlying bedrock) characterized by different permeability. The aforementioned conditions strongly influence slope stability [6,60]. The lack of maintenance of agricultural terraces and artificial drainage networks caused by the progressive abandonment of agricultural practices becomes one of the main factors in increasing LDS [6,54]. Generally, the highest amount of rainfall-induced slope instabilities occurred on abandoned agricultural lands, while a small portion of them affected cultivated terraced slopes. Moreover, slope failures affecting cultivated zones were characterized by a lower magnitude than those occurring on abandoned ones [6,84]. Furthermore, Cevasco et al. [6] reported that terraces abandoned for a short time showed the highest landslide susceptibility.

Slope angle is notoriously one of the most important predisposing factors, as gravity is the engine of landslide movements. Namely, the steeper the ground surface, the greater the magnitude of driving forces and the expected velocity of landslide movements. Indeed, more than 70% of the CTNP presents slope angles greater than 21°. Slope aspect plays a fundamental role in controlling soil moisture through solar radiation and rainfalls; in fact, the highest number of landslides are located on slopes facing SE-SW quadrants (72.3%). This could be explained with the presence of highly humid winds coming from the sea. Such winds increase humidity in south-faced slopes and, when heavy/intense rainfall occurs, the difference in water content of south-facing soil covers can play an important role in influencing triggering landslides [117,118]. Furthermore, north-facing slopes suffer from more severe weathering, since in the Northern Hemisphere these slopes preserve higher moisture than south-facing slopes. Therefore, it is expected that north-facing slopes have higher landslides occurrence compared to south-facing slopes. However, it should be noted that south-facing slopes may be subjected to greater landslide activity due to more frequent wetting–drying cycles compared to north-facing slopes [119].

Despite the complex geological and geomorphological setting of the study area, ML methods implemented in this work provided LDS maps of good quality, especially after the implementation of the Ensemble Method that boosted the performance of LDS assessment (AUROC = 0.9). This statistical approach allowed to consider a great amount of available input data. In particular, the easy implementation of a great number of predisposing factors over a large area made ML techniques a more reliable solution for LDS mapping. During the subsequent phase of this study, the reach angle identification was performed. The obtained results showed different modal values of the reach angle for channelled and open/planar slopes (28° and 36°, respectively). This difference, which is in agreement with several works reported in the literature [30,120], confirms that, on average, channelled landslides travel longer distances than landslides occurring on open/planar slopes. Moreover, it is worth noting that the obtained modal reach angle values are higher than those obtained in other areas [9,39]. This could be related to the fact that the reach angle depends on different factors (e.g., landslide volume, type of movement, land use, soil properties, water content, topographic constraints) that may differ from one setting to another. Given the high number of landslides selected in this study, it was not possible to identify with accuracy, among these factors, those responsible for the differences observed. Using the LDS map as a basis, the application of the geometric method of the “reach angle” allowed mapping of MI points and their susceptibility. Through the combination of the procedures used for LDS and LRS assessment, the LDTRS map was obtained. The comparison between LDS and LDTRS maps highlights how certain portions of the study area show an opposite degree of susceptibility. This especially occurs in flat areas located at the base of slopes that, despite a low LDS, could be affected by the flow-like landslides runout (Figure 8). Information in this sense were obtained comparing the propagation areas of shallow landslides occurred on 25 October 2011, at Vernazza with the degree of susceptibility resulting from the LDS and LDTRS maps (Figure 13). The results of this comparison were very encouraging, indicating a higher performance of LDTRS map with respect to LDS map (percentages of landslide deposits falling in areas with high and very high susceptibility of 82% and 48%, respectively). The availability of future events data could provide further information about the reliability and effectiveness of the LDTRS map obtained in this study.

6. Conclusions

In mountainous and hilly areas, rainfall-induced shallow landslides have always been important worldwide threats for settlements located at the base of steep slopes. During detachment phases, such landslides are often shallow, but subsequently they can evolve in higher-magnitude, dangerous mass flows with strong impacts on the social fabric. In this paper, a methodology aimed at producing a shallow landslide susceptibility map taking into account both detachment and runout was presented. The proposed methodology was applied to the Cinque Terre National Park, which is an area characterized by increasing shallow landslide risk. The novelty of this work consists of a processing chain that considers both the detachment and runout susceptibility assessment. Moreover, the use of input data simply derived from DEMs allows for achieving a good level of accuracy at such a scale of analysis and predictive efficiency also in case of a lack of exhaustive information. The results obtained in this work indicate that the susceptibility map taking into account both detachment and runout performs better than a simple detachment susceptibility map. Therefore, the proposed methodology may be particularly useful for areas where rainfall-induced shallow landslides constitute a threat to human lives, buildings and infrastructures. In particular, the obtained map will be helpful for urban and land planning, as well as for decision makers and stakeholders, to predict areas where rainfall-induced shallow landslides are likely to occur and travel in the future and to identify areas where hazard mitigation measures are required. Hence, this tool represents a step forward for a more proper territorial planning.