Landslide Susceptibility Mapping by Comparing GIS-Based Bivariate Methods: A Focus on the Geomorphological Implication of the Statistical Results

Landslide susceptibility is one of the main topics of geomorphological risk studies. Unfortunately, many of these studies applied an exclusively statistical approach with little coherence with the geomorphodynamic models, resulting in susceptibility maps that are difficult to read. Even if many different models have been developed, those based on statistical techniques applied to slope units (SUs) are among the most promising. SU segmentation divides terrain into homogenous domains and approximates the morphodynamic response of the slope to landslides. This paper presents a landslide susceptibility (LS) analysis at the catchment scale for a key area based on the comparison of two GIS-based bivariate statistical methods using the landslide index (LI) approach. A new simple and reproducible method for delineating SUs is defined with an original GIS-based terrain segmentation based on hydrography. For the first time, the morphometric slope index (MSI) was tested as a predisposing factor for landslides. Beyond the purely statistic values, the susceptibility maps obtained have strong geomorphological significance and highlight the areas with the greatest propensity to landslides. We demonstrate the efficiency of the SU segmentation method and the potential of the proposed statistical methods to perform landslide susceptibility mapping (LSM).


Introduction
Landslide susceptibility (LS) can be defined as the probability that a landslide will occur in a given area based on terrain conditions without any temporal consideration [1]. It follows the principle that future slope failures are more likely to occur under the same conditions that have triggered past and present slope failures [2]. Certain terrain conditions can be considered as predisposing factors for the development of landslides. Triggering factors are temporary conditions that can directly cause landslides (e.g., earthquakes, rapid snow melt, and heavy rain). Landslide susceptibility maps are produced by considering conditioning factors only, whereas landslide hazard maps consider both predisposing and triggering factors [2].
Many models for landslide susceptibility mapping (LSM) have been developed using different methods, scales, and evaluation criteria [3], from knowledge-driven to processand statistical-based models [4]. The most common and effective approaches in terms of the readability and usability of the output maps are based on more or less complex mathematical and statistical techniques, including (amongst others) logistic regression, neural network analysis, data overlay, index-based, and weight of evidence analyses, as well as machine learning (see [3] for a complete and up-to-date review). They require the a priori definition of the landslide inventory and predisposing factors and a rigorous procedure to calibrate and validate the models [5,6]. However, the complexity of the Remote Sens. 2021, 13, 4280 2 of 29 elaboration methodology often prevents or at least makes their reproduction and application to different areas difficult. The most used mapping unit for LSM is the grid [3], which is in a matrix form, can be easily obtained in GIS, and is simple to handle for data processing. However, grid cells are not associated with geological-geomorphological environments, as they cannot represent the physiographic conditions of terrain. A more representative segmentation method comprises the unique condition units (UCUs) [7], which are homogenous domains with morphodynamically constrained spatial limits [8]. They maximize the internal homogeneity and the external heterogeneity of the selected parameters and therefore better approximate the morphodynamic response of the slope to the occurrence of landslides [9]. A particularly powerful segmentation method is based on geomorphological features, such as slope units (SUs), which are obtained by splitting the two halves of subcatchments and considering the slope gradient and aspect [10]. Some studies have demonstrated the efficiency of this kind of terrain segmentation, which can even outperform grid-based models [11][12][13]. One of the main strengths of SU segmentation is that the obtained landslide susceptibility maps are more readable and directly linked to the terrain structure [14].
The topic of landslide susceptibility has been the object of numerous scientific studies, many of which have an exclusively statistical approach focused on comparing a number of different techniques on the same dataset without a discussion of the meaning of the predictive relationships or of their coherence with morphodynamic models. It is necessary to translate the models into forecast maps that are useful for administrators as an effective tool for risk management.
The aim of this study was to obtain a landslide susceptibility map at the catchment scale for a key area (populated, with active slope failures), based on simple statistical methods and SUs segmentation, characterized by a strong geomorphological significance. The entire elaboration procedure was based on a multidisciplinary approach in which LiDAR-derived digital terrain model (DTM) analysis, digital orthophoto interpretation, and GIS processing were combined with direct field surveys and past geomorphological map analysis.
We constructed a new simple and reproducible method for delineating SUs with an original GIS-based terrain segmentation method based on hydrographic features. The SUs were chosen to be large enough to contain the entire alimentation area of a landslide or a portion of that but small enough to be the reference slope area of a single landform.
We performed LS analysis by comparing two GIS-based statistical methods. We chose to apply bivariate statistical methods based on the landslide index approach [15], originally developed for grid-cell analysis, by adapting the equations. These two methods are simple and intuitive, and the advantage of using bivariate statistics is that they immediately reveal the relation of the predisposing variables with the propensity for landslides.
We selected the Piomba Stream basin, in Central Italy, as the test area. This is a strongly anthropized area, and the hillslopes are particularly prone to slope failures due to the geomorphological and geological features. This area has been studied for many years as an important site for the interaction of both human pressure and slope morphometry on erosion processes (e.g., [16][17][18][19]). Recent studies have focused on the role of precipitation on landslides in this area and evidenced that rainfall is the main triggering factor [20,21]. However, no studies have applied LS modeling.
The predisposing factors were selected after an intense and prolonged field survey in order to choose the most effective factors based on field evidence, picking the parameters most commonly used in the literature for evaluating LS [22,23]. We considered only one parameter for each terrain characteristic, avoiding redundancies, which could have generate overlapping effects and interrelationships that are difficult to discriminate [3]. Moreover, we tested the use of the morphometric slope index (MSI) [24] as a predisposing morphometric factor for landslides and the possibility of considering MSI as an alternative to the slope gradient and other slope morphometric variables. It is a single effective parameter for slope morphometry and is the result of the geomorphological dynamics, as it highlights all those processes that have changed the morphological structure of the terrain over time. It reflects the history of the slope and its modeling.
The entire LS analysis was framed within an experimental design that combined statistical methods, predisposing factors, and the areal threshold of landslides. The main outcomes of this research demonstrate the efficiency of the SU segmentation method and the potential of the proposed statistical methods for determining the LSM, producing easy-to-read susceptibility maps. Beyond the purely statistical values, the proposed models adhere to the geomorphological significance of the terrain and highlight the areas with the greatest landslide propensity.

Study Area
The Piomba Stream is located in the hilly coastal area of the Central-Eastern Apennines (Abruzzo Region, Italy; Figure 1). This W-E-elongated basin flows from Mt. Giove (747 m a.s.l.) to the Adriatic Sea, near the town of Silvi (PE). The drainage pattern is mostly trellislike and subordinately dendritic, with few small tributary channels. The main tributary is Fosso del Gallo in the hydrographic left side. as an alternative to the slope gradient and other slope morphometric variables. It is a single effective parameter for slope morphometry and is the result of the geomorphological dynamics, as it highlights all those processes that have changed the morphological structure of the terrain over time. It reflects the history of the slope and its modeling. The entire LS analysis was framed within an experimental design that combined statistical methods, predisposing factors, and the areal threshold of landslides. The main outcomes of this research demonstrate the efficiency of the SU segmentation method and the potential of the proposed statistical methods for determining the LSM, producing easy-to-read susceptibility maps. Beyond the purely statistical values, the proposed models adhere to the geomorphological significance of the terrain and highlight the areas with the greatest landslide propensity.

Study Area
The Piomba Stream is located in the hilly coastal area of the Central-Eastern Apennines (Abruzzo Region, Italy; Figure 1). This W-E-elongated basin flows from Mt. Giove (747 m a.s.l.) to the Adriatic Sea, near the town of Silvi (PE). The drainage pattern is mostly trellis-like and subordinately dendritic, with few small tributary channels. The main tributary is Fosso del Gallo in the hydrographic left side. The area is characterized by a Mediterranean climate with a sub-temperate littoral regime [25]. Mean annual temperatures are generally between 12.5 and 15.5 °C, and the mean annual rainfall is between 600 and 800 mm. Precipitation increases toward the mountains with a trend not proportional to the elevation. The maximum rainfall is in The area is characterized by a Mediterranean climate with a sub-temperate littoral regime [25]. Mean annual temperatures are generally between 12.5 and 15.5 • C, and the mean annual rainfall is between 600 and 800 mm. Precipitation increases toward the mountains with a trend not proportional to the elevation. The maximum rainfall is in autumn, with a secondary maximum in spring. Summers are rather dry, especially near the coast and on the hills with lower elevation.
The study area belongs to the outermost eastern sector of the Apennine chain, which features a NE-verging thrust belt geometry, dissected by several younger normal faults. The piedmont areas are mainly affected by extensional basins to the west and buried thrusts Remote Sens. 2021, 13, 4280 4 of 29 overlain by homocline plateau, mesa, and cuesta landscapes to the east [26]. The sediment outcrops are mainly foredeep siliciclastic sequences belonging to the Plio-Pleistocene succession, which become progressively younger moving toward the Adriatic Sea and are arranged according to a homocline slightly dipping towards the NE. They are mainly composed by clayey/sandy-clayey sediments, often interbedded with clastic deposits (sands and conglomerates with lenticular geometry, sometimes very thick) at various stratigraphic heights. In the areas near the coast, the hills' summits are covered by thick deposits of sands, gravels, and conglomerates of fluvial-deltaic and coastal environments [27,28]. At the highest elevations of the foothills, clayey marl and marly clay alternate at sandy and clayey levels. They belong to three main geological formations [29]: the Laga Formation (Messinian-Lower Pliocene p.p.), which includes the older and innermost turbiditic sequence (the westernmost part of the Piomba basin); the Cellino Formation (lower Pliocene), which includes a younger and more external foredeep turbiditic sequence (the central part of the Piomba basin); and the Mutignano Formation (upper Pliocene-lower Pleistocene), the post-orogenic sequence deposited within the Periadriatic Basin (the easternmost part of the Piomba basin) ( Figure 2). The homoclinal valleys are engraved along the clayey levels with marked asymmetry. The rivers, which flow according to the regional topographical gradient from W to E, have a cuesta morphology [17]. autumn, with a secondary maximum in spring. Summers are rather dry, especially near the coast and on the hills with lower elevation.
The study area belongs to the outermost eastern sector of the Apennine chain, which features a NE-verging thrust belt geometry, dissected by several younger normal faults. The piedmont areas are mainly affected by extensional basins to the west and buried thrusts overlain by homocline plateau, mesa, and cuesta landscapes to the east [26]. The sediment outcrops are mainly foredeep siliciclastic sequences belonging to the Plio-Pleistocene succession, which become progressively younger moving toward the Adriatic Sea and are arranged according to a homocline slightly dipping towards the NE. They are mainly composed by clayey/sandy-clayey sediments, often interbedded with clastic deposits (sands and conglomerates with lenticular geometry, sometimes very thick) at various stratigraphic heights. In the areas near the coast, the hills' summits are covered by thick deposits of sands, gravels, and conglomerates of fluvial-deltaic and coastal environments [27,28]. At the highest elevations of the foothills, clayey marl and marly clay alternate at sandy and clayey levels. They belong to three main geological formations [29]: the Laga Formation (Messinian-Lower Pliocene p.p.), which includes the older and innermost turbiditic sequence (the westernmost part of the Piomba basin); the Cellino Formation (lower Pliocene), which includes a younger and more external foredeep turbiditic sequence (the central part of the Piomba basin); and the Mutignano Formation (upper Pliocene-lower Pleistocene), the post-orogenic sequence deposited within the Periadriatic Basin (the easternmost part of the Piomba basin) ( Figure 2). The homoclinal valleys are engraved along the clayey levels with marked asymmetry. The rivers, which flow according to the regional topographical gradient from W to E, have a cuesta morphology [17].  The geomorphological evolution of the area was strongly influenced by the outcropping lithology, the tectonic structure, and the regional uplift, and by climatic-and sea-level changes [27]. This produced the relative sea-level changes and local base level changes (negative in the long term in comparison to today's sea level). The present landscape was shaped by the Pleistocene-Holocene geomorphological dynamics. During the Pleistocene, the rapid uplift and the alternation in climatic phases produced the deepening of the hydrographic systems and the dismantling of the summit conglomerate, which led to the outcrop of the clayey substrate [26]. Fluvial and slope processes were recorded in a complex sequence of continental deposits (post-orogenic) overlying the bedrock [30]. In the Holocene, fluvial, gravitational, and anthropic processes allowed the landscape to be shaped especially through mass movements (slides, flows, and falls), badlands, gully and rill erosion, slope remodeling, and de-and reforestation due to agricultural practices [16,17,19,20]. The mass movements and landslides are mostly active and recent, and the main triggering factor is high precipitation, the effect of which has combined with the general drainage dissection, terrain morphology (particularly slope gradient and concavity [21]), and land use [20].

Landslide Susceptibility Analysis
The data necessary to complete the landslide susceptibility analysis were obtained by means of a multidisciplinary approach that included geomorphological field survey, remote sensing (air photo analysis), DTM processing, and institutional open data. We used color aerial photos, regional topographic maps (Carta Tecnica Regionale, CTR, 2007) and DTMs provided by the Cartographic Office of the Abruzzo Region (Table 1) to derive the morphological and hydrological variables (i.e., the geometry of SUs, landslide inventory, land use, and drainage density). This improved the feasibility and repeatability of the analysis and will lead to future improvements with multitemporal analysis for the investigation of the trends in landslide susceptibility. The aerial photos were obtained with LiDAR flights, which were orthorectified and georeferenced with a pixel resolution of 0.2 m. The CTRs were produced using photogrammetrical techniques from the digital orthophotos at a 1:10,000 scale. The DTMs were LiDAR-based with a 10 m cell-size resolution. The geological and geomorphological variables were obtained through the field survey with the support of the Hydrogeological Setting Plan of Abruzzo Region (Piano per l'Assetto Idrogeologico, PAI) [31] and the Italian Landslides Inventory (Inventario dei Fenomeni Franosi in Italia, IFFI) [32] maps. The input data were collected in a geodatabase built in the GIS environment and organized in two categories: landslide inventory and predisposing factors ( parameter processing and morphometrical analyses were performed in ArcGIS (ArcMap ® 10.1, ESRI, Redlands, CA, USA). The data analysis was performed via bivariate statistics by comparing the results of the two landslide indexes. The SUs were considered as mapping units of the LS models.

Landslide Inventory
The landslide inventory was compiled through an integrated approach combining field and remote investigations. In the preliminary phase, pre-compiled regional maps were analyzed: PAI [31] and IFFI [32]. Afterward, the remote observation from aerial photographs allowed us to map the surface distribution of the landforms with GIS software. Finally, all the landforms were observed and measured in the field for ground-truthing.
A detailed geomorphological field survey was performed between spring 2017 and autumn 2019. In particular, from January to April 2017, heavy rainfall activated numerous shallow mass movements that were investigated in the field survey. The final product was the Geomorphological Map of Torrente Piomba catchment at a 1:10,000 scale (resized in a schematic map in Figure 3), from which the landslide inventory map and some predisposing factors maps were derived. The landslide inventory includes active/dormant features, inventoried in the GIS-based geodatabase, which are classified according to the Varnes [33] classification of landslides.

Predisposing Factors
The predisposing factors for landslides were chosen in relation to the terrain morphometry, geology and geomorphology, hydrology, vegetation, and anthropic features. The selection was based on field evidence, expert knowledge of the study area, a methodological approach from previous studies, and local data availability [22,23,34]. Seven controlling factors were considered from those that showed the greatest influence on landslides: bedrock lithology (L), land use/land cover (LULC), deposits thickness (T), slope aspect (A), drainage density (D), slope gradient (S), and morphometric slope index (MSI). Each factor was subdivided into different classes that represent their entire variability, according to the accuracy and scale of the data source. For the categorical factors (L and LULC), the majority class within the SU was considered. The numerical factors (D, S, and MSI) were calculated from vector layers (polygons, polylines) applying the appropriate formula. The slope aspect (A) was derived from the DTM, and the median was calculated in order to consider the most prevalent value in each SU. The deposit thickness (T) was extrapolated from morphometric variables and then subdivided into qualitative classes. Afterward, A and T were considered categorical variables.
Bedrock lithology (L). The bedrock lithology was selected as it is recognized to be one of the most important factors in LSM because it influences the slope stability, the size and the type of landslides, and subsequently the susceptibility degree [35,36]. In our study, the bedrock lithology map was derived from the field-based geomorphological map. We followed a lithotechnical approach, moving from a lithological basis and grouping the geological units into four classes of rock types ( Figure 4a): fluvial deposit, sandstone and conglomerate, clay, and alternating sand and clay. Each SU was assigned the prevailing class in terms of areal distribution.      Land use/land cover (LULC). The importance of LULC in slope stability is recognized in particular for areas subjected by a strong anthropic impact [16,37,38], such as the piedmont areas of Central Italy. The variability in vegetation cover can directly influence the soil infiltration capacity and consequently the slope stability [39]. Moreover, natural or artificial changes in LULC can weaken and destabilize the slopes. In our study, the LULC was derived from the Corine Land Cover project-Level 3 [40] and implemented through the analysis of orthorectified aerial photos. Nine types of LULC are relevant in the study area ( Figure 4b): vineyards, olive groves, complex cultivation patterns, land principally occupied by agriculture with significant areas of natural vegetation, annual crops associated with permanent crops, broad-leaved forest, non-irrigated arable land, transitional woodland shrub, and sparsely vegetated areas.
Deposits' thickness (T). The thickness of superficial deposits (i.e., mainly continental colluvial and slope deposits in the study area) was selected as it is recognized to be one of the strongest factors influencing surface processes determining infiltration rates, overland flows, and water storage potential [41,42]. These thickness data were derived from direct observations (e.g., boreholes) or from indirect estimation. In the study area, direct observations were not available in sufficient amounts to derive an extensive deposit thickness map to be able to predict their depth with an accuracy compatible with the requirements of LSM [43]. The available data were concentrated into small, urbanized areas, and could not be extrapolated for the entire basin. In these cases, an indirect geostatistical analysis can be applied to estimate the thickness through DEM-derived morphometric variables [43,44]. Therefore, we extrapolated T by using the equation developed by Sciarra et al. [45] for an area with similar bio-geo-climatic conditions. This is a linear pixel-based relation that links T to the local slope (tan β): The few data from direct measurements were used to check the values of T. The depths were validated by comparing the direct measurement data with the values of the estimation via the formula and making sure that the measured value was within the range of the estimated value. The most prevalent value of T was extracted in each SU as a median value via the zonal statistics. We divided T into two classes (Figure 4c) to discriminate between absent/very shallow and deep deposits (≤2 m and >2 m, respectively).
Slope aspect (A). The role of the slope aspect in slope stability is debated since different studies have reported opposite outcomes ( [46] and references therein). Some results suggested that A is an important predisposing factor for superficial landslides in clayey deposits and is correlated with other factors such as bedrock structure [46,47]. As already stated, we derived A from the DTM and calculated its median value in each SU using the Zonal Statistics tool. Then, we grouped the eight main directions with an angle step of 45 • (N, NE, E, SE, S, SW, W, and NW) (Figure 4d).
Drainage density (D). The drainage density is a traditional hydrological parameter that was developed by Horton [48] and represents the degree of fluvial dissection. It is strictly influenced by the slope gradient and the lithology. We selected D because it influences the slope stability and highlights hydrological critical conditions or critical fluvial erosion [49][50][51]. D is expressed as the length of the drainage network per unit area and calculated as the ratio between the total stream length and the basin area. The value of D can be influenced by the mapping scale and strictly depends on the channel head's source area [50]. We derived the drainage network as a vector layer from the Regional Topographic Maps of the Abruzzo Region at a 1:5000 scale (CTR, 2007) at the same scale as the other predisposing factors and verified and implemented this through aerial photo analysis. The hydrographic layers are LiDAR-based and accurate enough with respect to the terrain morphometry. Then, we calculated D for each SU. This factor was classified into five classes using the natural break algorithm (Figure 4e), with each corresponding to a certain range of D values (in this case, the ranges are 0, 0.0001-1.2, 1.2-2.5, 2.5-5, and 5-9.7 km −1 ).
Slope gradient (S). The slope gradient is the main morphometric variable that affects the slope stability at different scales (e.g. [3,52,53], amongst many others). It controls not only slope-gravity-induced processes but also hydrological processes, runoff, soil erosion, weathering, vegetation cover, and human activities, which in turn influence slope stability. We calculated S using the Add Surface Information tool, which attributed the average percentage of the slope to each SU and then transformed it into degrees. This factor was divided into four classes (0 • -10 • , 11 • -15 • , 16 • -20 • , and >20 • ; Figure 4f).
Morphometric Slope Index (MSI). The morphometric slope index (MSI) was recently introduced as a unique reference index to study the influence of pre-incision slope morphometry on the evolution of drainage networks and mass movements within Italian badlands and small clayey basins [19,24,54,55]. It groups the main (linear and areal) morphometric features of a morphological unit with the formula: where L is the length, R C is the circularity ratio, A 3D is the three-dimensional surface area, and A 2D is the plane surface area. In this study, we used MSI to characterize all morphometric features of the SUs and provide a more complete alternative to the use of S. It considers the slope's gradient (expressed through the relationship between the 3D and 2D area), shape, and size (length and width). The single parameters were calculated as vector (polygon) characteristics via specific tools in ArcGIS (Minimum Bounding Geometry for L, Interpolate Shape/Add Surface Information for A 3D , Calculate Geometry for A 2D , and the perimeter p used for R C = 4πA 2D /p 2 ). MSI was subdivided into eight classes (Figure 4g) by modifying the natural break algorithm ranges (11-35, 35-50, 50-90, 90-170, 170-230, 230-300, 300-400, and >400 m).

Mapping Units
In this study, we used the SUs as mapping units to assess landslide susceptibility. The SUs were derived using the Arc Hydro toolbox in ArcGIS, with little or no handling by the operator.
As a base layer, the LiDAR-based DTM with a 10 m resolution was used from the Abruzzo Region geodatabase (Figure 1). The step-by-step procedure to delineate the SUs is shown in Figure 5. The first phase was the terrain preprocessing, in which the DTM was corrected to clearly identify the drainage cells. The second phase was the terrain processing, in which the sub-basins were delineated and then divided into two by the main streams, in order to obtain the two halves of the catchment. To carry out this operation, the main drainage lines were extended until they intersected the watershed divides using the Longest Flow Path for Catchments tool. Subsequently, this longest flow path was used to cut the hydrographic unit into two SUs. Then, the aspect was derived via its median value in each SU using the Zonal Statistics tool. Using an iterative procedure, the smaller areas were aggregated to the larger adjacent ones, and the neighboring units in the same aspect range were joined. In this way, only the SUs greater than 10,000 m 2 with a uniform aspect were maintained.

Landslide Susceptibility Methods
To perform the LSM of the study basin, we chose to test two GIS-based bivariate statistical methods and to compare their results to establish which is the most powerful method for predicting landslides. The first is the landslide susceptibility index (LSI) developed by Romeo et al. [56] according to Lee and Min [57]; the second is a simplified version of LSI named the landslide index (LI) developed by Sciarra et al. [45]. The advantage of these methods is to provide an immediate measure of the role played by each factor and related classes on landslide susceptibility [56]. Both methods were designed to work on raster data and have been modified and adapted so that they can also be used on vector data. They measure a percentage-weighed sum of the predisposing factors ranging from 0 (least influence) to 100 (most influence). The complete workflow of the statistical analysis is shown in Figure 6. Remote Sens. 2021, 13, x FOR PEER REVIEW 13 of 33 Landslide Susceptibility Index (LSI) Originally, the LSI was created for raster-based LSM and was calculated by summing the ratio between the number of cells in which landslides occur and the number of cells in which landslides did not occur for each class of a factor. We transformed this calculation for the SU-based LSM into the ratio between the landslide area and the non-landslide area in the SUs for each class of a factor.
LSI is defined as the ratio between a and b, expressed in the following equation: The parameter a is the ratio between the landslide surface of each class (A(fi)) and the total landslide surface in the study area (A (F)), while parameter b is the ratio between the non-landslide area of each class (A(si)) and the total non-landslide area in the study area (A(S)). LSI was calculated for all the classes of each factor and then normalized to compare the data: where LSInor represents the value of LSI for each class (i) of a factor divided by the maximum LSI value in the same factor and multiplied by 100. It quantifies the impact of the classes taken individually. The total LSI was obtained from the combination of the LSInor for each SU: Landslide Susceptibility Index (LSI) Originally, the LSI was created for raster-based LSM and was calculated by summing the ratio between the number of cells in which landslides occur and the number of cells in which landslides did not occur for each class of a factor. We transformed this calculation for the SU-based LSM into the ratio between the landslide area and the non-landslide area in the SUs for each class of a factor.
LSI is defined as the ratio between a and b, expressed in the following equation: The parameter a is the ratio between the landslide surface of each class (A(f i )) and the total landslide surface in the study area (A (F)), while parameter b is the ratio between the non-landslide area of each class (A(s i )) and the total non-landslide area in the study area (A(S)). LSI was calculated for all the classes of each factor and then normalized to compare the data: LSI nor,i = LSI i /LSI i(max) × 100.
where LSInor represents the value of LSI for each class (i) of a factor divided by the maximum LSI value in the same factor and multiplied by 100. It quantifies the impact of the classes taken individually. The total LSI was obtained from the combination of the LSInor for each SU: LSI tot = ∑ (LSI nor,i × w j ) w j = LSI nor,i(min) /∑ LSI nor,i(min) . (5) LSItot is given by the sum of the LSInor of each factor multiplied by the weight of the factor (wj). The weight of a factor is expressed as the ratio between the value of the minimum LSInor of the classes and the sum of the minimum LSInor of all the factors. LSItot expresses the propensity for failure due to the weighted combination of the considered predisposing factors.
LSItot was assigned to each SU and the corresponding values w3re classified according to the Jenks' method (natural breaks). Five susceptibility classes were created: very high, high, moderate, low, and very low.
Landslide Index (LI) The LI was created to work in raster-based LSM, in the same manner as LSI. In this case, however, we applied the equation to the SU-based LSM as it is because LI is calculated as the ratio between the landslide area of each class (i) and the total surface of the class: The values obtained for the classes of each factor were normalized to 100: where LI i represents the value of the LI for each class, which is divided by the maximum LI value within the same factor and multiplied by 100. The total landslide index was obtained by the sum of the LI nor of each factor divided by the total number of factors (N): This expresses the propensity for instability due to the mathematical combination of predisposing factors. Furthermore, in this case, the SUs were classified and mapped according to the natural breaks method into five levels of susceptibility: very high, high, moderate, low, and very low.

Training and Validation
The training and validation phases were performed using the success rate and prediction rate curves to test the model accuracy and prediction skills.
For the training phase, the whole area of the basin of the T. Piomba without the basin of the Fosso del Gallo was considered as the test area (Figure 1), from which the susceptibility model was developed to determine the performance of the model. The calibration of the models was constructing by creating the ROC curves as plot charts, in which the cumulative percentage of the SU area for each susceptibility level is reported on the x-axis and the cumulative percentage of landslide areas for each susceptibility level is reported on the y-axis. To evaluate the reliability of the statistical method, the area under the curve (AUC) was calculated. The AUC can range between 0 (invalid) and 1 (high reliability). Generally, the model is considered valid if AUC ≥ 0.7 [58]. The calculation of AUC allows us to compare different models.
For the validation phase, the Fosso del Gallo sub-basin was chosen as the validation area ( Figure 1). This allowed us to evaluate whether the models can work even in areas other than those in which it was created and thus to test its capacity to predict future landslides. Once the model was obtained, it was applied to the validation area according to the same susceptibility levels, and its predictive performance was measured through the calculation of the AUC of the ROC curves.

Experimental Design
The research was organized in a 2 × 3 × 2 experimental design (Table 2) combining the statistical methods (2×), the combination of predisposing factors (3×), and the areal threshold of landslides (2×). Both statistical methods, LSI and LI (2×), were applied with different combinations of factors (3×), where one considered all the seven factors and the other two considered only one morphometric factor (alternatively S or MSI) together with the other six. This was carried out to test the hypothesis that MSI alone is sufficient to describe the morphometric characteristics of the SU, without the need to consider other factors, such as S. Furthermore, for each application, two models (2×) were employed according to the areal threshold of landslides: one considering all the mapped landslides and the other considering the SUs with less than 2% of the area covered by landslides as free of landslides [59]. A total of 12 landslide susceptibility maps were generated and then compared through the training and validation procedure to establish which achieved the best performance in terms of reliability and prediction skills.

Landslides
In the landslide inventory map (Figure 7), a total of 265 landslides were reported, covering an area of 15.5 km 2 , corresponding to 14.8% of the basin (105 km 2 ). Many of the landslides occurred on cultivated sites, and their morphological evidence was modified or removed by farming practices. Therefore, the pre-existing PAI and IFFI inventories allowed us to identify the correct location of some of the mapped landforms whose boundaries were correctly traced during the field activity and the air-photo analysis. In addition, more recent landslides were mapped during the field activity and the air-photo analysis. As observed after the rainfall events, the main triggers of mass wasting in the study area were heavy rainfalls, which caused the activation or the reactivation of most of the landslides.
The majority of the landslides were of the flow type (223 out of 265). In the areas in which they occurred, the substrate was mainly clayey, favoring their (re-)activation in the rainy periods due to water infiltration. There were 24 rotational landslides, leading to a number of counterslopes mostly set at the foot of the concave sliding surfaces. All the 16 complex landslides had a rotational component, and most of them were located near the source area of T. Piomba. A single landslide with a translational component was located on the right bank of the Fosso del Gallo. There was also a rockfall on the left bank of the Fosso del Gallo below the town of Atri, involving the summit conglomerates.
In test area, the number of SUs with landslides amounted to 259; while consid an areal threshold of 2% of landslides, the total was 251. In the validation area, the nu of SUs with landslides amounted to 40; while considering the areal threshold of 2 landslides, the total was 37. Each SU was assigned a class of factors.

Landslide Susceptibility Models
The LSI and LI statistical methods were applied to the six different situations o experimental design combining areal threshold and predisposing factors ( Table 2). A of 12 susceptibility maps were obtained: 6 maps for each statistical method. The m were coded by adding a consecutive number as a suffix to LSI and LI. The suscepti maps were arranged into five classes, and LSItot values were divided by means of na breaks each time (Figures 9 and 10, Tables A7 and A8).

Characterization of the SUs
The T. Piomba basin was divided into 613 Sus (Figure 8): 518 in the test area and 95 in the validation area. In the test area, the smallest SU was 0.012 km 2 , the largest was 1.05 km 2 (mean = 0.17 km 2 , median = 0.118 km 2 ), and most of the SUs were between 0.012 and 0.4 km 2 (96.3%). In the validation area, the smallest SU was 0.015 km 2 , the largest was 1.8 km 2 (mean = 0.2 km 2 , median = 0.117 km 2 ), and most of the SUs were between 0.015 and 0.4 km 2 (91.6%). This suggests that SU sizes were comparable in the test and validation areas. In test area, the number of SUs with landslides amounted to 259; while considering an areal threshold of 2% of landslides, the total was 251. In the validation area, the number of SUs with landslides amounted to 40; while considering the areal threshold of 2% of landslides, the total was 37. Each SU was assigned a class of factors.

Landslide Susceptibility Models
The LSI and LI statistical methods were applied to the six different situations of the experimental design combining areal threshold and predisposing factors ( Table 2). A total of 12 susceptibility maps were obtained: 6 maps for each statistical method. The models were coded by adding a consecutive number as a suffix to LSI and LI. The susceptibility maps were arranged into five classes, and LSItot values were divided by means of natural breaks each time (Figures 9 and 10, Tables A7 and A8). In test area, the number of SUs with landslides amounted to 259; while considering an areal threshold of 2% of landslides, the total was 251. In the validation area, the number of SUs with landslides amounted to 40; while considering the areal threshold of 2% of landslides, the total was 37. Each SU was assigned a class of factors.

Landslide Susceptibility Models
The LSI and LI statistical methods were applied to the six different situations of the experimental design combining areal threshold and predisposing factors ( Table 2). A total of 12 susceptibility maps were obtained: 6 maps for each statistical method. The models were coded by adding a consecutive number as a suffix to LSI and LI. The susceptibility maps were arranged into five classes, and LSI tot values were divided by means of natural breaks each time (Figures 9 and 10, Tables A7 and A8).

LSI Models
The LSI bivariate statistics described the influence of the classes of each factor (Tables A1 and A3). Most of the landslides developed in SUs with alternating sand and clay bedrock and broad-leaved forests (and subordinately olive groves), with a high deposit thickness (>2 m), northern aspect (primarily NW, and subordinately N and NE), low D (primarily 0.0001-1.2 km −1 , and in general below 5 km −1 ), medium-high slope gradient (16-20 • , and in general above 11 • ), and high MSI (>400 m). This is valid if we consider the areal threshold of 2% for the landslide surface.
The weighted equation for LSI tot was different for each case because the weights of the factors (Tables A2 and A4) assume different values according to the influence of the factor. In any case, the models produced similar results in terms of the relative importance of the factors. D, MSI, and T were the most influential factors; A and S were moderately influential; and L and LULC were the least influential. Comparing the models that considered both MSI and S (LSI_1 and LSI_2), MSI was always found to be the most influential morphometric variable. In the cases in which MSI was not included in the analysis (LSI_5 and LSI_6), the influence of S was found to be secondary. In the cases where S was not included in the analysis (LSI_3 and LSI_4), the influence of MSI was even higher.
In general, the susceptibility maps ( Figure 9, Table A7) show that the SUs with very high and high susceptibility were mainly located on the right hydrographic side of the T. Piomba. The SUs with moderate susceptibility were located on the left side. The SUs with low and very low susceptibility were mainly located on the left side and in the areas close to the mouth.

LI Models
The LI bivariate statistics describe the influence of the factor classes and not the relative importance of the factors because the equation for LI tot depends on the number of factors, assuming that they have the same weight (Tables A5 and A6). The relative influence of the classes is similar to the LSI models, as well as the distribution of the susceptibility in the SUs. There were only small differences among the models concerning the absolute values of the indices and not their relative relationship. The most influential factor classes were alternating sand and clay bedrock and broad-leaved forests (and subordinately olive groves) with a high deposit thickness, a NW aspect (and subordinately N and NE), low D (<5 km −1 ), medium-high slope gradient (>11 • ), and high MSI (>400 m). Very-highand high-susceptibility SUs were mainly located on the right hydrographic side of the T. Piomba, while moderate susceptibility SUs were on the left side, and low and very low susceptibility SUs were on the left side and close to the mouth ( Figure 10, Table A8).

Training and Validation
The training and validation procedure was performed via the ROC curves (Figure 11), as already described. The AUC values of the LSI models were between 0.65 and 0.67. Therefore, the model training accuracies were between 65% and 67%. The AUC values of the validation set were between 0.66 and 0.69, reaching a prediction accuracy between 66% and 69%. Not one of the LSI models reached the minimum threshold of the AUC value for reliability, with values slightly below 0.7 [58]. The LSI model with the highest AUC values was LSI_1, which considered all the parameters and all the landslides without the areal threshold of landslides. In general, the LI models reached the highest validity, with the best performance for models LI_7 and LI_8. Concerning the different combinations of predisposing factors, the models with higher validity were those that considered all the parameters, while the areal threshold of landslides had no influence.   The AUC values of the LI models were between 0.69 and 0.71. Therefore, the training accuracies were between 69% and 70%. The AUC values of the validation set were between 0.67 and 0.68, reaching a prediction accuracy between 67% and 68%. Almost all LI models reached the minimum threshold of the AUC value for reliability, with values around 0.7 [58]. The LI models with the highest AUC values were LI_7 and LI_8, which considered all the parameters, with and without the areal threshold of landslides.
In general, the LI models reached the highest validity, with the best performance for models LI_7 and LI_8. Concerning the different combinations of predisposing factors, the models with higher validity were those that considered all the parameters, while the areal threshold of landslides had no influence.

Discussion
In this study, we applied the LS method in a given study area (Piomba basin, Central Italy), providing a simple and straightforward way to statistically analyze the propensity for landslide based on the distribution of the predisposing factors. A multidisciplinary approach was applied that combined LiDAR-derived DTM analysis, aerial photo interpretation, GIS processing, past geomorphological maps analysis, and direct field surveys. This approach allowed us to obtain an up-to-date and reliable landslide inventory map and to select the predisposing factors with the appropriate accuracy and resolution. The use of remote sensing data provided us the possibility to validate the field survey in order to obtain ground-truth maps; moreover, the comparison of the present (field) with the past (aerial photos and DTM) terrain features highlighted the main changes in the terrain, allowing us to more efficiently detect and delimit the landforms.
The study area was segmented into SUs by developing an original GIS-based procedure starting from the terrain hydrography. The landslide index [15] approach via bivariate statistics was applied by adapting the grid-cell analysis to the vector analysis.
The outcomes of the statistical analyses can be summarized as follows: -The two statistical methods applied (LSI and LI) can be used interchangeably because both of them describe the relative influence of each class of each factor; - The method that reached higher statistical validity is the LI method (the simplified version of LSI), which resulted in a slightly higher value of AUC; -With regards to the accuracy level, the models' performance is sufficient for only a few of the factors and is almost equivalent for the two statistical methods (LI and LSI); - The best predictive models were LI_7 and LI_8 (created with the LI statistical method considering all the parameters with and without the areal threshold of landslides, respectively), reaching a prediction accuracy of 71%; - The results of the different models are substantially in agreement when detecting the factor classes that have primary importance for the development of landslides in the study area: the landslides are mainly located on slopes characterized by alternating sand and clay, occupied by broad-leaved forests, with NW aspect, gradient between 15 • and 20 • , MSI values greater than 400 m, and low D (between 0.0001 and 1.2 km −1 ); - The only significant difference observed between the two methods is the distribution of susceptibility classes in the study area; in particular, the very high susceptibility class characterizes a greater number of SUs in LSI maps than in LI maps; - The LS models created by considering the complete landslide inventory do not differ from those that consider the landslides with an areal threshold of 2%, probably because there are few landslides under the threshold of 2% and these landslides do not affect the results; -There are no differences between the models created by considering six factors (alternatively with S or MSI) rather than seven, indicating that MSI is effective in representing the morphometric characteristics of the SU: it includes a more effective quantification of the slope features related to the geomorphological processes distribution on the SUs; S characterizes the SUs only through the value of the gradient (in degrees), while MSI includes a set of morphometric parameters such as slope, shape, and size. This result confirms previous studies on the validity of MSI as morphometric factor, e.g., [50].
From a statistical point of view, no further information is needed to explain the results. However, this distribution can be better understood if the relationships between landslides and predisposing factors are analyzed in detail from a geological and geomorphological point of view. Some general considerations can be extrapolated from the local level (i.e., study case) focusing on the geomorphological implication of the statistical results.
Considering the morphometrical setting, the most influential classes of S were found in general to be above 11 • , with a slight predominance of the 16-20 • class. This indicates that intermediate and high slope gradients (>11 • ) can lead to the formation of landslides in the study area, and few landslides form on low gradient hillslopes (<11 • ). The most influential class of MSI has an MSI value >400 m, with secondary classes of 230-300 m and 301-400 m. For SUs with an equal slope gradient, these classes correspond to longer and wider hillslopes, which hold most of the landslide area. Generally, longer hillslopes are also wider, which can induce the formation of flow-type and rotational landslides on soft clay-sand bedrock, such as the Apennine hills [16,60]. The effect is directly controlled by the water infiltration capacity and terrain saturation associated with the specific hydrogeological and lithological setting. The presence of local aquifers in sand sequences confined by undelaying clay deposits [16] and the subsequent high values of neutral pressure can result in increasing flow-type landslides [61].
The slope gradient values are connected to the geo-lithological context of the study area, which is mainly constituted by soft lithologies (clay, and alternating sand and clay). In the study area, low-gradient (i.e., ≤10 • ) SUs correspond to alluvial deposits, while medium-high-gradient (i.e., >11 • ) SUs correspond almost equally to the clay, alternating sand and clay, and sandstone/conglomerate layers. The presence of thick layers of sands within the clay can allow the slope to have a higher slope gradient. Moreover, the slope can be linked to the thickness of slope deposits derived from a linear relation with the local slope [45]: T is higher when S is lower [43].
The drainage density is different within the SUs: in many SUs, the value of D is 0 km −1 , indicating the absence of streams within the SUs; a number of SUs have a D value >5 km −1 ; and most of the SUs have D values in classes ranging from 0.0001 to 5 km −1 . The primary influential class was that with lower D values (0.0001-1.2 km −1 ), but all the classes with a D value <5 km −1 had almost the same probability of being affected by landslides. Therefore, there is no distinctive class of D that influences the landslide proneness in the study area. This outcome can again be linked to the geo-lithological setting: in an area with a quite uniform distribution of lithological units (mainly constitutes by soft lithologies), the distribution of D is also uniform and reflects the landslide proneness in relation to the landslide types (flows and rotational).
Considering the LULC factor, it is homogeneous and poorly differentiated in the study area; therefore, its effect can only be local. The most influential class was found to be the broad-leaved forests. This is related to the fact that they occur only on small portions of the basin, mainly located in the upstream area. Almost all the SUs with this prevalence are affected by landslides. The secondary LULC class was found to be olive groves; in this case, the landslide clusters in the SUs also fall within this type of cultivation. The territory is principally occupied by sparsely vegetated areas, complex cultivation, and land principally occupied by agriculture, with significant areas of natural vegetation, which were found to have a low impact on the occurrence of landslides. However, in these areas, the superficial evidence of landslides can be modified by human intervention, especially by agricultural practices, which, in many cases, can temporarily obliterate the landforms [16]. Moreover, the construction and the maintenance of the country roads in some cases can affect the predisposition of a slope to small landslides or, vice versa, on the obliteration of the superficial signs of the landslides.
Regarding the slope aspect, the most relevant orientation was found to be NW, and secondarily NE and N. This parameter has often been correlated with the effect of thermal excursions and the more humid and wet condition of clayey sediments and superficial colluvial deposits with respect to the southern-facing hillslopes affected by wet/dry cycles [46]. However, in our case, the link between landslide distribution and slope aspect is more probably connected to the local bedrock structural setting, as already mentioned in previous studies [46,47]. In particular, in the main (homoclinal) valley, the cuesta structure caused the north-facing slopes (N, NE, and NW) to have a dip direction and the south-facing slopes (S, SE and SW) to have an anti-dip direction, while in the secondary (cataclinal) valleys, the E-facing and the S-facing slopes show a cross-dip direction. Generally, dip slopes have a greater propensity for landslides, while anti-dip slopes have greater stability. In this case study, the slope aspect can be considered an indirect factor of bedrock local geological structure [46].
Considering the other geomorphological processes that occur in the study area and that were included in the Geomorphological Map (Figure 3), some interesting considerations can be noted. In many cases, the SUs with a higher MSI and S (and subsequently lower T) correspond to the side slopes of badland basins. In general, they have a steep and rough morphology (which is quantified by high MSI values). This is due to the badland erosion process (i.e., linear incision inside the main gullies) [18,23,62] over drainage basins and slopes with a dense hydrographic network, and is favored by the lithology, in which the presence of thick layers of sands within the clay can allow the conservation of a high slope gradient. This determines that the corresponding SUs have also high D values (e.g., >5): in these cases, the main geomorphological processes are linear erosion (rill erosion, gully erosion, and piping) with a secondary role of the very shallow landslides (i.e., mud flows), which tend to fill the main channel [19,54].
The predisposing factors that reached higher importance in the statistical analysis are not necessarily the most important in the general geomorphological analysis, and vice versa. This is the case for D and L if we consider the LSI models. In this specific study area, they are poorly statistically significant due to the homogeneity of the terrain features, indicating that the predisposing factors are spatially associated, and their importance is dependent on local terrain features. From a geomorphological perspective, the most important predisposing factors are the geolithological features (L), to which other parameters are related. The morphometrical (MSI and S), morphological (A), and hydrographic (D) parameters are directly linked to the lithology of the terrain, which in turn is indirectly linked to other geological characteristics, such as the geological structure.

Conclusions
The results of this research study, through the creation of an up-to-date and detailed geomorphological map of the Piomba Stream basin, led to the development of an LSM that shows the areas most prone to landslides. The work was based on a multidisciplinary approach that combined direct field surveys, past geomorphological maps, LiDAR-derived DTM analysis, digital orthophotos interpretation, and GIS processing. The use of remote sensing data and applications can allow the landslide susceptibility assessment to be improved, in particular for the landslide inventory's creation and the predisposing factors, ensuring the ground truth of the field-based mapping.
Two statistical GIS-based bivariate statistical methods were applied based on the landslide index approach [15] through the simple transpositions of raster methods into vectors. An original and simple method for delineating the SUs was presented that is fully consistent with the physiography of the terrain and the conceptual image of the segmentation of the landscape into SUs [9,14]. Furthermore, the morphometric parameter MSI was successfully considered as a predisposing factor in LS analysis for the first time. The outcomes confirm the possibility of using MSI as the only morphometric factor in landslide susceptibility analysis, which was shown to be more representative of the slope's morphometry compared to S, especially in relation to the distribution of linear and areal erosion processes on the SUs.
The results of the statistics reached an adequate level of accuracy. An in-depth analysis and discussion of the results was conducted to link the outcomes of the statistical analysis to the geological and geomorphological features of the study area from a geomorphological point-of-view. The relationships between the predisposing factors and their consistency with the morphodynamic models were highlighted: the predisposing factors are strictly related with each other and with the local terrain features. We confirm that the geo-lithological factors are the main drivers for landslides in the study area, whether we consider them directly or indirectly as variables of the statistical models [3,4]. The predisposing factors used in the landslide susceptibility model of the Piomba Stream basin are spatially associated [4,58], and their importance is strictly dependent on local terrain characteristics [23,63]. These considerations can be extended to other basins of the Periadriatic foothills area, which has similar characteristics in terms of geological and geomorphological features.
This study highlighted the importance of examining the outcomes of the statistical analysis in a geomorphological perspective. The most important variable in the models is the lithology, which is directly related to the morphometry (MSI and S), morphology (A), and the hydrography (D) of the terrain and indirectly related the geological structure. However, some parameters play a priority role in the statistical analysis but a secondary role if we consider their geomorphological importance, and vice versa. This is the case with D and L, mainly due to the local conditions related to the homogeneity of the terrain features in the study area.
The observed similar statistical relevance values confirmed that both the LS methods (LI and LSI) have similar predictive skills [58] and can be used interchangeably. They provide quantitative indications, supported by geomorphological observations, of the propensity to landslide, reaching an adequate validation level. They must be further expanded and investigated, but surely can be useful if associated with some geomorphological considerations, in particular concerning the influence of one parameter on the others and the local terrain characteristics. The approach used in this study may be applied to similar morpho-geo-climatic areas and contexts, such as the Adriatic hilly area, as confirmed by previous work [45].