Combining multi-typologies landslide susceptibility maps: a case study for the Visso area (central Italy)

ABSTRACT The research proposes a simple but geomorphologically adequate method to produce a combined landslide susceptibility map. In fact, in a logic of real use, offering type-specific landslide susceptibility maps to land use planners and administration could be not a successful solution. On the other hand, the simple grouping of more types of landslides could be misleading for model calibration considering that the relationships between slope failures and geo-environmental predictors should be conveyed by the abundance of each type of landslide resulting not specific and diagnostic for each typology. In this test, after having produced independent models for flow, slide and complex landslide by exploiting MARS (Multivariate Adaptive Regression Splines) and a set of type-specific geo-environmental variables, a combined landslide susceptibility map was obtained by combining the scores of the three source maps. The combined map was finally validated with a new unknown archive, showing very good performances.


Introduction
One of the main issues of researchers in landslide hazard mapping, and specifically, landslide susceptibility evaluation study is to offer an easy-to-apply/read map for the final user (land use planners and administration, Bufalini et al., 2021;Martinello et al., 2021;Smith et al., 2011).This reflects not only in the adequate partitioning of the study area, which should define mapping units geomorphologically suitable both for the modelling procedures and for obtaining final landslide susceptibility maps (Martinello et al., 2021(Martinello et al., , 2022a)).In fact, considering the different typologies of landslides (Hungr et al. 2014) and the different influence of each factor in landsliding, specific modelling procedures should be analysed, and the related final maps produced.At the same time, the management of several landslide susceptibility maps, typespecific, could be not easy for local administration.For this reason, efforts should be addressed to obtain combined landslide susceptibility maps, grouping different typologies of phenomena according to geomorphologically adequate criteria.In this sense, a frequently adopted approach is to produce an a priori grouped model which means aggregating landslides of more than one typology for calibrating the model.However, generally, these studies do not focus their attention on different behaviours and characteristics of grouped phenomena or, even more, on different the relationships between the geo-environmental predictors and each landslide type.
In this study, an integrated landslide susceptibility map combining flow, slide and complex types is presented for the Visso (Marche, Italy -Main Map, figures a and b) area, which was obtained by combining the final score of the three source independent susceptibility maps.The original landslide inventory was obtained by the reviewing the existing CARG (Geological and Geomorphological Mapping of Italy, Regione Marche 2014) landslide dataset for the 'Visso' map.In facts, mapped single phenomena were often grouped inside large polygons, with a low resolution both in terms of the spatial pattern and typological characterization.For this reason, by exploiting topography maps and orthophotos, the original polygons were split according to the movement typology thus producing separated slides, flows, and complex landslides archives (198, 91, and 51 cases, respectively).By exploiting the LCL_SLU (Martinello et al., 2021(Martinello et al., , 2022a) ) as mapping units and Multivariate Adaptive Regression Splines (MARS) statistical method, three specific susceptibility models were prepared by regressing each landslide inventory on a specific set of physical-environmental predictors.Furthermore, multicollinearity and variable importance analyses were carried out to verify their relevance and influence in landslide susceptibility assessment.After attesting the sharing of some of the most important variables for the three models, a combined landslide susceptibility map was obtained by adding the three scores of the three single maps.The final integrated map (whose combined score was scaled to a range between 0 and 1) was finally submitted to double validation: the first with respect to the modified CARG archive (putting together slide, flow and complex archives); the second, by exploiting a new inventory obtained by remote checking on Google Earth.The high performance of the model suggests a geomorphologically suitable final map.

Geological, geomorphological and climatic setting
The Umbria-Marche Apennines constitutes the external sector of the central Apennines, a fold-and-thrust belt, NE verging, formed by the convergence between the Corsica-Sardinia and the Adria continental margins.This sector of the Apennines is characterized by arcuated contractional structures (Neogene folds and thrusts faults) involving a Mesozoic-Tertiary sedimentary sequence (Figure 1).
The study area consists of a ∼ 30 km 2 sector of the Sibillini Massif, that characterizes the southern portion of the Umbria-Marche Apennines.The geological framework is characterized by the Umbria-Marche succession (Calamita & Deiana, 1988;Pierantoni et al., 2013) composed by a multilayer of pelagic and emipelagic formations, with alternating calcareous, marly calcareous, calcareous-siliceous and siliceous rocks.
The physical landscape is characterized by a predominantly mountainous territory, with reliefs that exceed 1800-2000m a.s.l.(2051m a.s.l. at Mount Bicco).The slopes are often steep, with angles that exceed even 70°, while the valley floors are narrow and strongly incised.
In these conditions the continental deposits (Pleistocene-Holocene in age) are reduced both in extension and in thickness and are mainly made up of slope deposits, sometimes cemented, debris and alluvial fans, and fluvial deposits (Gentili et al., 2017).Gravitational morphogenesis is also very important and represented by rockfalls, slides (both rotational and translational), and flows (mainly debris flows); very frequent are the complex phenomena usually made by slides that evolve downstream in rapid flows (Aringoli et al., 2010;Buccolini et al., 2020;Materazzi et al., 2021) (Figure 2).
From a climatic point of view, the study area is characterized by an Apennine-Adriatic regime with rainfall almost uniformly distributed throughout the year; the highest values are recorded late in autumn and during spring while the lowest ones, in July and January (Gentilucci et al., 2020(Gentilucci et al., , 2021)).The total annual precipitation often exceeds 1500 mm with 1000 mm isohyet encompassing the whole mountain area.A snowy contribution from November to April is also present; the quantity and persistence of snow increase with the altitude.The study area shows average temperatures of almost 11°C, with a significant increase of 0.5 C in the last 30 years on average, due to climate change (Gentilucci et al., 2019).

Landslide inventory dataset and diagnostic area
For this study, the data from the landslide inventory of the Marche Region in the year 2000 were initially collected and homogenized, later also used for the editing of the Regional Geological Map at 1: 10,000 scale (Marche Region, 2014).Subsequently, to verify their evolution over time, integration and re-measurement of the landslide areas were carried out using a colour orthophoto of 2020 at the same scale.For the classification of landslides, according to Cruden and Varnes (1996), three main categories (slides, flows, and complex) and four states of activities (new activation, inactive, reactivated, and partly reactivated), were assumed; among these, only one seismic-induced landslide as a consequence of the 2016 central Italy earthquake seismic sequence has been recorded (Aringoli et al., 2021;Farabollini et al., 2018;Romeo et al., 2017;).Fall and topple as well as deep-seated landslide were here excluded as literature generally employs a direct and/ or physical-based approach (e.g.Agnesi et al., 2015;Cafiso et al., 2021;Cafiso & Cappadonia, 2019;Cappadonia et al., 2019Cappadonia et al., , 2021;;Pappalardo et al., 2021).
According to the inventory review, 351 phenomena (Figure 3) were included, whose summary characteristics are given in Table 1.
The data collected show that the most numerous phenomena fall into the 'slides' category although complex landslides are those that, referring to the single phenomenon, have the greater areal extension.The most elevated landslide density is, instead, that of flows which represent, among other things, almost all the new phenomena (90%) that took place between 2000 and 2020.
The recognized landslides were mapped by using a Landslide Identification Point (LIP), which corresponds to the highest point along the crown of the landslide area here assumed as diagnostic in potentially marking unstable slope conditions (Cama et al., 2015;Lombardo et al., 2014Lombardo et al., , 2016;;Rotigliano et al., 2011Rotigliano et al., , 2018Rotigliano et al., , 2019)).

Mapping units and landslide conditioning factors
According to Martinello et al. (2021Martinello et al. ( , 2022a)), the LCL_SLU were used as mapping unit.This type of slope unit is obtained by intersecting the classic hydromorphological unit by means of 5000 cells contributing area threshold with the classes of the Landform Classification (LCL; Guisan et al., 1999).The LCL was obtained by exploiting the TPI-based Landform Classification (LCL; Guisan et al., 1999) SAGA tool, by fixing the inner/outer radius as 100/1000 metres.For more information about the LCL_SLU, please refer to the original research of Martinello et al. (2021) and Martinello et al. (2022a).
Each LCL_SLU was then classified as stable or unstable (positive or negative cases) depending on whether it hosts at least one LIP.
Based on the expected direct or proxied role in landsliding (Costanzo et al., 2012;Mercurio et al., 2021;Rotigliano et al., 2018Rotigliano et al., , 2019; Vargas-Cuervo  2 gives for each of the selected predictors the main characteristics and the potential proxied significance. Before starting the modelling procedures, the DEM-derived variables were submitted to multicollinearity analysis, based on the evaluation of the Variance Inflation Factor (VIF) obtained by applying the 'usdm' R-package (Naimi, 2017).All the variables passed the multicollinearity test.According to the specific characteristics of the landslide types analysed in this research, all the predictors were employed for complex type modelling, SPI predictor was excluded for the slide model, while the TWI variable was excluded for the flow model.
To assign the variables to the slope units, the continuous variables were attributed by deciles, while relative frequencies were employed for categorical ones.

Statistical model and validation tools
Between the different approaches offered by the literature, several researchers have been exploiting the  Multivariate Adaptive Regression Splines (MARS; Friedman, 1991) method to favourably determine nonlinear relationships and complex interactions in several fields of the science (e.g.Cervantes et al., 2020;Conoscenti et al., 2021).In the last ten years, the MARS method successfully was employed for landslide susceptibility evaluation (e.g.Conoscenti et al., 2015Conoscenti et al., , 2016;;Felicísimo et al., 2013;Martinello et al., 2021Martinello et al., , 2022aMartinello et al., , 2022b;;Mercurio et al., 2021;Rotigliano et al., 2018Rotigliano et al., , 2019;;Wang et al., 2015), by regressing the outcome (stable/unstable status) onto the covariates set from the controlling factor layers.
MARS is a non-parametric regression method that exploits the splitting of each independent variable into hinge functions to boost the maximum likelihoodbased adaptation skill of the logistic regression method, according to where y is the dependent variable (the outcome) predicted by the function f (x), α is the model intercept, and β i is the coefficient of the h i basis functions, given the N number of basis functions.By exploiting the 'earth' R-package (Milborrow, 2014), the MARS analysis was performed.Since this type of statistical modelling relies on a balanced assessment of status cases, a random extraction of negative in the same number of positive cases was carried out.For testing the independence of the result from the selection of negative/positive cases: (1) the extraction of negative cases was replicated one-hundred times and the models' accuracy and precision evaluated (Costanzo et al., 2014;Martinello et al., 2021;Martinello et al., 2022aMartinello et al., , 2022b)); (2) according to Chung and Fabbri (2003), random splitting procedures were employed to verify the prediction skill of the model by using 75% of the balanced archive for the training phase and the remaining 25% for the validation one.
The performances of the models were evaluated by adopting both cut-off independent and dependent metrics.In particular, the AUC value (Area Under Curve) in the ROC (Receiver Operating Characteristics -Fawcett, 2006;Goodenough et al., 1974;Lasko et al., 2005) was employed to evaluate the prediction skill of the model and to obtain the Youden index optimized cut-off (Youden, 1950).Then, by exploiting the optimized cut-off the true/false positive/negative cases (i.e.TP, TN, FP, and FN, respectively) were distinguished and confusion matrices evaluated.

Model building and validation strategies
According to the aim of this research, single models for each type of landslide were first prepared by using the three 2000 inventories.In this way, three source models and related landslide susceptibility maps were prepared and validated: reg_slide, reg_flow and reg_complex, for slide, flow and complex landslides, respectively.A combined map was then obtained by adding the three source scores and rescaling the combined value on the 0-1 range.Finally, the prediction skill of the combined map was tested in detecting both all the 2000 landslides and, according to a temporal partition, the unknown 2020 phenomena.

Results
In Figure 3 the boxplot for the three single models reg_slide, reg_flow and reg_complex, both for the validation to the blinded 25% (prediction) and the success validations are reported.In particular, based on replicates we produced through multiple random positive/ negative extraction, the distribution of the one-hundred AUC values is plot, so to estimate accuracy and precision of each of the three models.The results attest good and excellent (Hosmer & Lemeshow, 2000) AUC mean values for the prediction (0.78, 0.81 and 0.8, respectively for reg_flow, reg_slide and reg_complex).However, medium standard deviations are associated with these AUC mean values (from 0.05 to 0.09).Better performances are obtained from success models, with outstanding AUC mean values (0.95, 0.92 and 0.99, respectively for reg_flow, reg_slide and reg_complex) and very low standard deviations of replicates (from 0.01 to 0.02).
High performances are also revealed by the confusion matrices (Table 3) both for prediction and success validation.For the latter, the increase in the performances is more marked, especially for the complex type model.The sensitivity for the prediction validation is just below 0.8 while the values increase for success validations ranking from 0.82 to 0.97.Similar behaviour is reported for the specificity which attests around 0.7 for the prediction validations, increasing up to 0.97 in success validations.
Table 4 shows the 10 most important variables for the reg_flow, reg_slide and reg_complex success models.It is important to note that this type of analysis highlights only how frequently the class of the variables is employed for the construction of the model.Therefore, the analysis shows if a class of the variable is used to discriminate positive to negative cases by the model, but it is not diriment for defining if the class is positively or negatively correlated with landsliding.In this sense, the old landslide deposits class (lito12) is a very important class for all types of models considering that its overall ranks from 95 to 100 (even employed for the model construction).The max value of the elevation is a very important class of variable for reg_complex and reg_flow models.These models also share the natural grassland class of USE as important.The second most important class for reg_slide, the max value of PRF, is also the 5th for reg_flow while the Q100 of PLC, in common between reg_slide and reg_complex, is recalled around 10 times.
In figure c of the main map, the combined landslide susceptibility map for the Visso area is proposed.The four different landslide susceptibility classes are defined by employing nested Youden index cut-offs (low = 0.25, mean = 0.39, high = 0.61).For a more manageable definition of susceptibility, the ranges are linked to nominal classes described as S1 or 'low', S2 or 'moderate', S3 or 'High' and S4 or 'Very High' susceptibility levels.
The validation of the final map with respect to the 2000 inventory shows good performances with an AUC value of 0.87 (Figure 4(a)) and a very high Potential hydrological and surface hydric erosion induced disturbances (Martinello et al., 2022a) Table 3. Confusion matrix of the three single models both for prediction and success procedures.

Discussion and conclusion
Excellent performances result for the three single models (reg_flow, reg_slide, reg_flow), both in prediction validation and in the success tests.The violin plots show a larger but limited (std.dev.< 0.1) dispersion of the AUC values for the prediction analysis.On the other hand, excellent and stable AUC values were obtained in success analysis.These high performances were confirmed also by the cut-off dependent matrices.The analysis of the importance of the variables suggested lito12 ('old landslide deposits') is the more employed class for the construction of the models, for all the typologies investigated.This suggests that most of the landslides of the 2000 inventory (that correspond to the calibration archive) are a reactivation of old phenomena.The comparison shows also that other variables or classes of predictors are in common between the three models: prf (max value), ele (max value) plc (max value) and use10.
On the other hand, some classes of geo-environmental predictors are essential only for single specific models: lito5 (clayey marls) is important for flow type, lito18 ('scaglia bianca' clays) for slides while more use of dem-derived variables resulted for complex type.These behaviours encourage the grouping of the three types of landslides, suggesting also that grouping in the calibration phase could have hidden some of the relationships between a specific movement and its related geo-environmental predictors.On the contrary, by adding together the susceptibility score obtained by the single modelling procedures, the specificity of each model and, in this way, the peculiarity of landslide typology can be noticed.At the same time, the combined landslide map showed high performances in validation both with respect to the 2000 and the unknown 2020 inventories.For the latter, it is worth noting a slight lowering of the performances.According to the analysis of the variable's importance, the models are more sensitive to reactivation movements.Considering the characteristics of the two landslide inventories, this lowering could be linked to the state of activity of the phenomena: with In red, the predictors are employed by all three models: in light blue, by only two models.Legend: lito12 = old landslide deposits; lito5 = clayey marls, lito9 = calcareous units, lito19 = 'scaglia cinerea' clays, lito18 = 'scaglia bianca' clays, corine10 = natural grassland, corine8 = coniferous forest, corine7 = broad-leaved forest, lcl0 = streams, lcl5 = open slopes.respect to the 2020 year, the most of 2000 landslides are in facts still dormant.The results obtained in this application on the Visso area allow us to consider the approach employed for producing a combined landslide susceptibility map for the flow, slide and complex landslide types as geomorphological suitable.

Software
The landslides were remotely verified and remapped by using Google Earth Pro.In the first step for the production of LCL_SLU, GRASS GIS (GRASS Development Team, 2022) was utilized to obtain the SLU by using the r.watershed processing.In the second step, Saga GIS 8.0.0 (Conrad et al., 2015) was employed to produce the LCL map and to clip the SLUs with the LCL map.Saga GIS 8.0.0 was also employed to obtain the dem-derived variables and to rasterize the variables available in.shp format.Quantum GIS3.8.2 Zanzibar (QGIS.org,2022)was used to print all maps displayed in this research.The models were implemented by using the RStudio (RStudio Team, 2020) software.

Figure 1 .
Figure 1.Geological sketch of the study area.

Figure 2 .
Figure 2. Mountain portion of the Nera River basin: landslide inventory used in the present study.

Figure 3 .
Figure 3. Violin boxplots of the three single models both for prediction and success procedures.

Figure 4 .
Figure 4. Validation indices of the combined map for the validation with respect to the 2000 inventory and the 2020 one.(a) ROCcurves and AUC values and, in the bottom, confusion matrices; (b) degree of fit plot.

Table 1 .
Summary table of the landslide inventory.

Table 4 .
Most important variables for the three main models.