Earthquake-triggered landslide susceptibility in Italy by means of Artificial Neural Network

The use of Artificial Neural Network (ANN) approaches has gained a significant role over the last decade in the field of predicting the distribution of effects triggered by natural forcing, this being particularly relevant for the development of adequate risk mitigation strategies. Among the most critical features of these approaches, there are the accurate geolocation of the available data as well as their numerosity and spatial distribution. The use of an ANN has never been tested at a national scale in Italy, especially in estimating earthquake-triggered landslides susceptibility. The CEDIT catalogue, the most up-to-date national inventory of earthquake-induced ground effects, was adopted to evaluate the efficiency of an ANN to explain the distribution of landslides over the Italian territory. An ex-post evaluation of the ANN-based susceptibility model was also performed, using a sub-dataset of historical data with lower geolocation precision. The ANN training highly performed in terms of spatial prediction, by partitioning the Italian landscape into slope units. The obtained results returned a distribution of potentially unstable slope units with maximum concentrations primarily distributed in the central Apennines and secondarily in the southern and northern Apennines. Moreover, the Alpine sector clearly appeared to be divided into two areas, a western one with relatively low susceptibility to earthquake-triggered landslides and the eastern sector with higher susceptibility. Our work clearly demonstrates that if funds for risk mitigation were allocated only on the basis of rainfall-induced landslide distribution, large areas highly susceptible to earthquake-triggered landslides would be completely ignored by mitigation plans.


Introduction
The concept of landslide susceptibility defines the expectation of where landslides may occur in a given landscape, thus providing information on the spatial component of the landslide hazard definition (Varnes and the IAEG Commission on Landslides and Other Mass-Movements, 1984;Guzzetti et al. 1999;Lombardo et al. 2020a). The numerical expression of a landslide susceptibility corresponds to the probability of landslide occurrences within a given mapping unit (Lima et al. 2017;Broeckx et al. 2018;. As reported by Bird and Bommer (2004), earthquake-triggered landslides often cause the largest damage related to earthquakes and affect transportation routes inhibiting recovery and safety operations during emergency phases (Martino et al. 2019;Mantovani et al. 2019). In detail, earthquake-triggered landslides caused more than 50% of the total worldwide losses due to landslides in the last decades (Petley 2012).
To preemptively reduce the risk associated with these processes, predictive models have been proposed to estimate the distribution of earthquake-triggered ground effects scenarios (Sassa 1996;Jibson et al. 2000;Prestininzi and Romeo 2000;Romeo 2000;Wasowski and Del Gaudio 2000;Del Gaudio et al. 2003;Jibson 2007; Hsieh and Lee 2011; Lombardo and Tanyaş 2020, among others) representative of a uniform hazard distribution or seismic shaking scenarios. Out of several available options proposed to preemptively estimate earthquake-triggered effects, the proposed approaches essentially boil down to two types: physically based approaches (Van Westen et al. 2006) and statistically based ones (Guzzetti et al. 2005). The first type of approach implies that slope stability analyses are performed to quantify safety factors (Martino 2016) and/or the expected seismically induced displacements of the landslide masses and include Newmark's method (Del Gaudio and Wasowski 2004). The recently proposed PARSIFAL (Probabilistic Approach to pRo-vide Scenarios of earthquake-Induced slope FAiLures) approach allows obtaining probability maps of expected Newmark's displacements at regional scales (Esposito et al. 2016;Martino et al. 2018Martino et al. , 2019. As regards the statistically based counterpart, the general framework is quite similar when data-driven (statistical and machine learning) models are used. Typically, a mapping unit is chosen between grid-cells and slope-units and a dichotomous status expressing the absence or presence of landslides (or 0/1) is assigned. In a subsequent step, the binary status is fitted to a set of predictors chosen to represent predisposing factors of slope instability and the outcome of the modeling procedure is a probability (Amato et al. 2019). However, the algorithmic architecture one chooses to implement has notable repercussions on the performance each model provides. For instance, simple bivariate statistical models provide quite straightforward interpretation of the functional relations existing between factors and landslides (e.g. Weight of Evidence; Bonham-Carter 1989;Van Westen 2002;Martino et al. 2019). But, this is achieved at the expense of statistical rigour (the model does not assume any underlying probability distribution nor the interaction among explanatory variables) and performances. For instance, multivariate statistical routines assume that landslides are distributed over space according to the Bernoulli probability distribution . And, they allow to model linear relations (in the case of Generalised Linear Models; Ayalew and Yamagishi 2005;Castro Camilo et al. 2017) or a combination of linear and nonlinear relations (in the case of Generalised Additive Models; Brenning 2008; Goetz et al. 2011) between predisposing factors and landslides occurrences. These models offer excellent performance while keeping a clear interpretability at each step and for each model component (Lombardo et al. 2014;Frattini et al. 2010). Ultimately, machine learning methods provide equally and often even higher performance than the other two approaches mentioned above, this time though at the expense of the interpretability of each step. The reason behind this characteristic is due to the fact that machine learning algorithms are often based on the combination of highly nonlinear functions which are difficult to be individually and multivariately traced as the model evolves converging to the best solution (Liu et al. 2014;Zhou et al. 2016Zhou et al. , 2018. Because of the high performance provided, the number of machine learning applications n landslide science rose in recent years (Marjanović et al. 2011;Huang et al. 2017;Zhu et al. 2017;Amato et al. 2021). And, Artificial Neural Networks (ANNs) (e.g. Ermini et al. 2005;Gomez and Kavzoglu 2005), including their more recent convolutional extensions (e.g. Wang et al. 2019) have demonstrated to be a valid tool for landslide susceptibility assessment. Neural networks (Hassoun et al. 1995) are characterised by the possibility of modelling the relationship between independent and dependent variables in a complex non-linear way and are by nature prone to overparameterization of the model itself. These aspects lead to the advantages that ANN have the ability to model complex relations when these are not known a priori and in the fact that, thanks to overparameterization, they are very little sensitive to problems of collinearity (De Veaux and Ungar 1994) between independent variables (i.e. the internal dependence among variables that may cause that one variable explains the variability of another, or more than, one). Thus, ANNs approaches are particularly suitable for big data.
Most of the works that adopt ANN approaches to analyse earthquake-triggered landslide susceptibility perform analysis at a regional scale (Song et al. 2012;Umar et al. 2014;Zhou and Fang 2015) using input landslide inventories that are limited in time and space to single earthquakes (Tanya et al. 2017;Shrestha and Kang 2019;Tanya and Lombardo 2020). A classical approach to perform earthquake-triggered landslide susceptibility analysis by means of ANN techniques usually sees, as a first step, inventorying landslides triggered by a selected large earthquake. Afterwards, a susceptibility model, driven by the inventoried landslides, is produced for the region hit by the selected earthquake (Song et al. 2012;Umar et al. 2014;Zhou and Fang 2015;Shrestha and Kang 2019). As a result, the susceptibility analysis ends to have a relatively local impact, being hardly exportable to other seismic regions, even within the same country. This is mainly because the model is earthquake-specific. To generalize the model and make it more exportable, landslides triggered by other relevant earthquakes, even from surrounding regions, should be included in the starting inventory.
In this study, we performed a susceptibility analysis of earthquake-triggered landslides byimplementing an Artificial Neural Network approach for the whole Italian territory. Slope units have been adopted as mapping units (Alvioli et al. 2020;Titti et al. 2022a;Loche et al. 2022a). The input landslides inventory used to train the network has been accessed via the Italian Catalogue of Earthquake-Induced Ground Failures (CEDIT; Fortunato et al. 2012;Martino et al. 2014;Caprari et al. 2018), which collects ground effects caused by earthquakes occurred over the whole Italian territory from the XII century to present days. As far as the authors know, this represents the first study dealing with earthquake-triggered landslide (EQtLs) susceptibility for the whole Italian territory. To clearly evaluate and present the achieved results, we quantified the ANN classification performances through commonly adopted metrics, and we generated the first Italian EQtLS susceptibility map. Besides, we investigated the importance of the predictors performing a Permutation Feature Importance (PFI) of single predictors and explored how the classification performance varies by selecting all the possible combinations among predictors groups (i.e. terrain, seismic, geothematic, hydrological, and anthropic predictors). Ultimately, we checked the obtained susceptibility map for every Italian administrative region, by using an additional landslide dataset, which was not included during the ANN training phase. The resulting percentage of unstable territory for every Italian region has been computed to highlight priorities in land management practices at more local scales other than the national one.

Italian morphotectonic settings
Italy is the European country most affected by landslides (Herrera et al. 2018), with over 620,000 landslides recorded in the framework of the IFFI dataset, the most complete and detailed landslides inventory existing in Italy (Trigila et al. 2013). The main triggering factors for landslides in Italy are intense rainfalls and earthquakes. And in recent years, anthropogenic factors such as road cuts have been assumed to also play an increasing role. From 1999 until December 2019, almost two billion Euros were dedicated by Italian institutions to interventions for landslide risk mitigation, which great majority was realised within or close to areas classified with a high landslide hazard (link here). This highlights the high relevance that the classification of the territory assumes to establish remediation funding priorities. Italy is also characterised by active geodynamics related to the geological evolution of the two major mountain chains, i.e. the Alps in the north and the Apennines throughout the peninsula, as testified by the distribution of earthquakes and volcanic activity. More specifically, the Alps' chain shows a double-verging growth, involving the exhumation of metamorphic rocks. Conversely, the Apennines chain consists of a single-east-verging belt, mostly characterised by thinskinned tectonics. As a consequence, earthquakes show prevalent compressional focal mechanisms at the fronts of the two chains and extensional mechanisms along the Apennines backbone (Carminati et al. 2010). The highest magnitude seismic events, with peak ground acceleration (PGA) values higher than 0.225 g and a return time of 475 years, are expected in the central-southern Apennines, Calabria region (on the southwest of the Italian peninsula), in the southeastern part of Sicily island and in the north-eastern sector of the Alps chain. Medium to low seismic acceleration values (PGA up to 0.225 g) are expected with a return time of 475 years along the entire Alpine Arch, along the entire western Italian coast and the peri-Adriatic regions (eastern Italian coast). Ultimately, ttablehe Sardinia island is the only sector with very low seismic hazard (link here).

CEDIT catalogue
The EQtLs susceptibility model we built in this study is based on data collected in the CEDIT catalogue (Prestininzi and Romeo 2000;Fortunato et al. 2012;Martino et al. 2014Martino et al. , 2019Martino et al. , 2020aCaprari et al. 2018), which is available online at (link here). The CEDIT collects all the ground-failure effects produced by earthquakes that occurred in Italy taking information from historical documents and literature, especially in the case of earthquakes that occurred before 1980, while ground effects induced by more recent events have been surveyed directly on the field by the CERI (Research Centre for the Geological Risks of Sapienza University of Roma) working group (see Martino et al. 2017, for more details on the standard cataloguing procedure). Among the many information contained in the CEDIT, the error estimation assigned to each ground effect location according to the following scheme (Martino et al. 2014) was found particularly useful for this work: • Class 5: high-quality location from historical documents or GPS measurement with any or negligible associated error; • Class 4: location coordinates with an average error of 1 km; • Class 3: location coordinates with an average error of 3 km; • Class 2: location coordinates with an average error of 10 km; • Class 1: location coordinates with an average error of 30 km.
The EQtLs were extracted from the CEDIT catalogue and were split into two different subsets ( Fig. 1): 1. An "Input dataset" containing 1545 landslides, all belonging to the georeferencing class 5. These were induced by the earthquakes that occurred in Italy from 1908 to 2018. 2. A "Check dataset" containing 465 landslides with georeferencing classes ranging from 1 to 4, induced by all the earthquakes contained in the CEDIT catalogue, and 54 landslides belonging to the georeferencing class 5, induced by earthquakes that occurred in Italy before 1908.
With the aim to provide the susceptibility analysis with a high-quality geolocalised dataset, we used the Input dataset for the training and cross-validation-test cycles of the neural network, whereas we used the Check dataset to perform an a-posteriori and independent verification of the EQtLs susceptibility map of Italy. The spatial distribution of the two here considered datasets (i.e. earthquake-induced landslides for Input and Check) are shown in Fig. 2.

Mapping unit
A mapping unit in landslide science is considered to be a geographical object upon which the landscape is partitioned (Carrara 1988). The vast majority of the landslide Fig. 1 Bar chart showing the distribution of the type of EQtLs for Input dataset (a) and for Check dataset (b), together with the georeferencing class distribution for the Check dataset only (b). "Generic Landslides" means that the landslide type is not known Fig. 2 a Spatial distribution of EQtLs belonging to the Input dataset and b of EQtLs belonging to the Check dataset, coloured on the basis of their georeferencing class susceptibility literature is based on regular mapping units shaped as a squared (e.g. Jibson et al. 2000;Steger et al. 2020) or hexagonal (Avolio et al. 2013;Lupiano et al. 2018) lattice. However, when it comes to statistically based applications, the way these units are used is generally flawed for a few reasons. These have been extensively described in Reichenbach et al. (2018) and we direct the reader to this article for more details. Most of these issues do not affect a valid alternative represented by a slope unit (SU) partition which is mapping units bounded by ridges catchment/sub-catchment divides and streamlines (Carrara et al. 1991(Carrara et al. , 1995. SUs are able to capture the variability of the landscape associated with the failure process, by maximising homogeneity of slope steepness and aspect within a single unit and heterogeneity of the same between adjacent SUs (Alvioli et al. 2016;Titti et al. 2021). On the other hand, a SU choice requires an additional step which corresponds to the aggregations or upscaling of properties that are represented over space with a much higher resolution. Oftentimes, mean and standard deviation values are extracted for numerical properties at the scale of the single SU (e.g. Guzzetti et al. 2006). But, these could also be expressed via different summary distribution metrics, e.g. such as quantiles (Amato et al. 2019). Similarly, it is not standardised the way categorical properties such as lithology or land use are aggregated at the SU scale and different approaches can be found in the literature (e.g. Castro Camilo et al. 2017;Schlögel et al. 2018).
In this study, we select a SU partition of the Italian territory. In addition to the above-mentioned reasons, for such a large study area, choosing a small regular lattice would have inevitably produced several tens of millions of grid-cells. The alternative of seeking a reasonable size of the dataset would have instead produced grid-cells which would have been individually very coarse, neglecting any geomorphological representation of the landscape under study. The SUs we used were made available by Alvioli et al. (2020) at the following address (link here). In Fig. 3, we summarised the distribution of all the SUs' planimetric areas, ranging from approximately 0.1 to 10 km 2 .
To support the analyses in this study, we assigned unstable conditions' labels to SU intersected by landslides contained in the Input dataset while we assigned a stable label to all the others.

Predictor variables
To support the modeling protocol at a scale comparable to the whole Italian territory, we selected a broad set of predictors aimed at expressing properties known or assumed to influence landslide occurrences. InTable 1 in the Appendix, we provide a general overview of these predictors by grouping them into geological, seismic, anthropic, terrain, and hydrological macro-classes. Before describing these properties, we stress to the reader that in order to aggregate each predictor at the SU scale, we used the mean and standard deviation criterion for continuous properties as well as the dominant class for categorical properties. And we recall that, if collinearity among predictors, which may negatively affect the model by inflating the variance estimates (McElroy and Jach 2019), exists, ANNs are able to handle it by spreading the estimated weights over the collinear variables to take into account the different noise levels, taking actually advantage in terms of predictive performance. Therefore, we have chosen to keep the whole predictor set which consists of 167 layers (more details will be provided in "Artificial Neural Network").

Geothematic predictors
We considered three geo-thematic properties, detailed below: 1. Landforms are specific geomorphic features on the earth's surface which encompass both large-scale terrains such as plains or mountain ranges and small-scale characteristics such as single hills or valleys (Jacek 1997). The work of Guisan et al. (1999) first and Jenness (2006) later has pioneered the automatic extraction of such features from DEMs. More recently, Jasiewicz and Stepinski (2013) have implemented an efficient automatic classification tool for landforms named geomorphon, which returns 10 terrain The 95% confidence interval is calculated as the difference between the 97.5 and the 2.5 percentiles of the SU area morphologies in 10 classes: (1) flat, (2) summit, (3) ridge, (4) shoulder, (5) spur, (6) slope, (7) hollow, (8) footslope, (9) valley, (10) depression. In this study, we used geormophon to initially calculate the ten landforms, and in a subsequent step, we have aggregated this information at the SU scale by assigning to a given mapping unit the class with the largest planimetric extent. 2. Similarly, we have assigned to each SU the predominant lithological type. This geological information was retrieved from the Geological Map of Italy at 1:500,000 scale. This map was based on 1:100,000 and 1:50,000 national geological cartography or geological maps (Tacchia et al. 2005). Overall, after the aggregation step, 21 lithology classes have been assigned to SUs across the whole Italian territory. Table 2 in the Appendix offers a description of each class. 3. The predominant soil type was assigned to each SU on the basis of the European Soil Map compiled by the European Commission-Joint Research Centre (Finke and Montanarella 2001). In this map, soil type classes are classified according to the World Reference Base (WRB) system, which consists of a two-levels terminology. The first level defines the Reference Soil Groups whereas the second level is nested within the first and consists of a set of principal and supplementary qualifiers (for more details, see link here). In this study, SUs have been classified on the basis of 91 soil type classes, which have been used for modeling purposes and mainly belong to the Reference Soil Groups reported in Table 3 in the Appendix.
Concerning the categorical predictors, slope units have been labelled with "1" in correspondence with the predominant classes of geomorphon, soil type and lithology, and "0" for all the other classes.

Seismic predictors: distance to seismogenic features
Seismic information has been considered in the form of Euclidean distance to the nearest active fault and the Euclidean distance to the nearest seismogenic source. Specifically, for each SU, the mean distance value and its standard deviation have been computed. Data required to produce these predictors have been accessed from the Italian Database of Individual Seismogenic Sources (DISS) catalogue (link here), which continuously updates the study on the Italian seismogenic sources. An Individual Seismogenic Source is obtained by parameterising the geometry and kinematics of large active faults considered capable of generating earthquakes with a magnitude (M w ) greater than 5.5 (Basili et al. 2008;DISS-Working-Group 2018). This corresponds to an active fault that has accumulated some displacement in the recent past and can be considered very likely to produce a new offset in the near future (link here).

Terrain predictors
Concerning Terrain predictors, we used the 20m DEM released by the Italian Institute for Environmental Research in 2013 (link here). And, for each slope unit, we calculated the mean value and the standard deviation of the following derivatives: • Elevation (e.g. Ayalew and Yamagishi 2005) can be considered a proxy for climate-related characteristics (e.g. ground temperature or even the precipitation itself when high ridges play the role of meteorological barriers). And, its standard deviation per slope unit mimics the signal of surface roughness. • Eastness and Northness, these are computed as the sine and cosine of the aspect expressed in radians, respectively ). These are two linear components of the nonlinear slope exposition signal, a common proxy for strata attitude and localised dry/wet soil conditions. • Slope gradient (Zevenbergen and Thorne 1987) expresses the potential gravitational forces acting over a given slope. • General, Longitudinal, and Tangential Curvatures (Evans 1980;Wood 1996); Planar and Profile Curvatures (Heerdegen and Beran 1982). Plan and profile curvatures carry the signal of the potential soil availability, and potential small-scale hydraulic and gravitational forces (Ohlmacher 2007). Conversely, cross-sectional curvature measures the curvature perpendicular to the downslope direction. As a result, it detects small-scale features such as channels. Longitudinal curvature plays a similar role but parallel to the downslope direction (Patel and Sotiropoulos 1997). • Topographic Positioning Index (TPI, De Reu et al. 2013) measures the difference between the elevation of a focal cell and the average elevation within a predetermined radius. • Topographic Roughness Index (TRI, Riley et al. 1999) expresses rough terrain conditions. • Topographic Wetness Index (TWI, Beven and Kirkby 1979) expresses the terrain's tendency to retain water at a given location, as a function of local slope steepness and upslope contributing areas. Therefore, it conveys information related to potential high pore pressure conditions distributed over the landscape or the presence of open floodplains. • The area of each slope unit (A SU ) controls the availability of potential material to fail.

Anthropic predictors: distance to roads
An ideal situation to inform any predictive model of the potentially destabilising effect of road cuts would be to collect the exact location and height of the cut. However, such information is available only for the location component and no height characteristics can be accessed for the whole Italian road network. For this reason, we opt to compute the Euclidean distance from roads at buffers equal to 5, 10, 50, and 100 m. Subsequently, a series of statistical metrics of the distances to roads have been calculated for each SU, namely the mean, maximum, and minimum distance of the unit from the closest road and the portion of the territory extending within certain distance ranges. Therefore, the following statistics have been calculated for every slope unit using the Zonal Statistics Plugin, in QGis 3.10.4 (Graser 2016).
• Count: the count of the number of pixels at a < 100-m distance; • Sum: the sum of the pixel distance values; • Mean: the mean distance; • Min: the minimum distance; • Max: the maximum distance; • Range: the range (max-min) of distance; • Majority: the most represented distance within a slope unit; • Count < 5: the count of the number of pixels at a < 5 m distance.

Hydrological predictors: distance to watercourses
The Euclidean distance from watercourses has been computed similarly to the road network case. This time though, we extracted ten equally spaced (100 m wide) buffer zones from 0 up to 1000 m from each streamline. The same summary statistics calculated for the distance from the road network have been computed also for the hydrological network with respect to each slope unit.

Artificial Neural Network
The used ANN architecture was previously proposed in Amato et al. (2021) and has been here optimised to perform a binary classification between stable and unstable slope units. Stable slope units are those SUs with no EQtLs while unstable SUs contain at least one landslide of the Input dataset. The ANN training is performed on balanced class datasets (Titti et al. 2022b). The used network is a "shallow" ANN whose architecture is a two-layer fully connected feed-forward network. For the hidden layer, a sigmoid activation function has been considered. The output layer is a "softmax layer", in which the outputs are normalised into probabilities proportional to the exponentials of the input values (Fig. 4). The network is trained by scaled conjugate gradient backpropagation. To limit any overfitting effect an "early stopping by validation" training criterion has been adopted. The classification process associates a probability value, from 0 to 1, to each slope unit to be susceptible to EQtLs. In order to be correctly trained to distinguish between stable and unstable slope units, the ANN needed to learn from samples of both classes. We set a fixed number of samples per class (equal to the number of all the slope units with landslides). Therefore, the Input dataset counted 644 positives (i.e. slope units with landslides) and an equal number of negatives (i.e. slope units without landslides), the latter chosen randomly from the larger number available. The Input dataset was then split as follows: 70% of samples were used to train the network, 15% was used for validation and 15% as the test dataset. The training dataset is used to optimise the weights and the bias assigned to each node of the ANN. After each step of the iterative training, the ANN classification is applied also on the validation dataset and the classification performances on the two datasets are monitored. As the classification performance continues to improve on the training dataset but worsens on the validation dataset, the training process is early stopped and overfitting of the model is avoided. Finally, the test dataset is a completely independent dataset used to test the reproducibility of performances obtained on the first two sets. In order to build a statistically significant distribution of the classification results and performance metrics, we replicated the training procedure 100 times. To ensure the maximum statistical independence, for each of the 100 replicates, the training, validation, and test datasets are recreated from scratch as described before. Furthermore, the initial values of ANN weights and biases are randomly changed. Fixed the ANN architecture, some of the operating network hyperparameters, and in particular, the number of nodes in the hidden layer, have been tuned to achieve the best and more reliable performances. In the "tuning" tests, the ANN performance was calculated as True Positive Rate (TPR, or Recall). TPR is the ratio between the number of true positives (i.e. those samples correctly predicted by the model as belonging to the given class) and the sum of true positives and false negatives (i.e. those samples that the model predicted as not belonging to a given class while they belong). A number varying from 1 to 6 nodes in the hidden layer has been tested. It resulted in a TPR increase as the number of nodes increased. The number of nodes was finally set to 4 as being the smallest number of nodes, which still produced a significant increase in performance. At the end of each of the 100 training replicates, the ANN was run on all the SUs, covering the whole national territory. The mean of the probability values output from the 100 classification replicates, as well as their standard deviation, was calculated and was used to plot the Earthquake-induced Landslide Susceptibility Map of Italy.

Performance assessment: validation routines
In the case of binary classification, a probability threshold is set to associate a particular sample to one of the two possible classes. The most appropriate way to investigate the discriminatory capabilities of a binary classifier for each possible value of the discrimination threshold between 0 and 1 is commonly the receiver operating characteristic (ROC; Rah-mati et al. 2019) plot. ROC plots report the TPR on the y-axis and the false positive rate (FPR or fall-out) on the x-axis. FPR is defined as the ratio between false positives and all the negatives, namely false positives + true negatives. False positives are samples classified as belonging to the class of interest while they were not, whereas true negatives are those samples correctly predicted by the model as not belonging to the class of interest. The Area Under the Curve (AUC) is strictly linked to the shape of the ROC curve, and it is a good proxy of the overall capability of a model to distinguish between two classes, regardless of what classification threshold is chosen. AUC assumes values between 0 and 1, gradually increasing with the classification capabilities of the model. If AUC is 1, the model is perfectly able to distinguish between positive class and negative class (Hosmer and Lemeshow 2000).
In this study, a probability threshold of 0.5 has been chosen to classify each slope unit as stable or unstable. This choice is confirmed by examining the point of the average ROC corresponding to a threshold value of 0.5 (as also reported in Fig. 5a). This point is the closest one to a TPR equal to 1 and an FPR equal to zero. Once the threshold value has been chosen, it is possible to further investigate the obtained discrimination capabilities by the means, for instance, of a Confusion Plot Lombardo et al. 2020b). Conversely to ROC (and AUC), the Confusion Plot is a threshold-dependent method to evaluate classification performance. It shows, fixed the classification threshold, the value of TPR as a function of the corresponding TNR for each of the 100 ANN replicates, allowing to evaluate the statistical robustness of the classification to different training, validation, and test datasets. In model performance evaluation, TNR stands for True Negative Rate and is the ratio between the number of true negatives and the sum of true negatives and false positives. In this study, TNR refers to the success rate in classifying slope units as belonging to the "stable" class, and TPR refers to the "unstable" one. Besides Confusion Plot and ROC (plus AUC), which are considered good indicators of the general performance of a model and commonly adopted in the scientific literature , the performance obtained by the network in this study has been evaluated also representing the importance assumed by each predictor during the classification by performing a Feature Importance analysis. This procedure highlights those predictors that gave a major contribution to the success of the susceptibility analysis. To make this, the Permutation Feature Importance (PFI) was adopted. The method is based on the assumption that a random variation of the value of an important predictor has a negative impact on the performance of the model greater than that of the random variation of a less important predictor (Putin et al. 2016). The permutation allows the random variation of the predictor while preserving the natural distribution of the values of the predictor itself (Gao et al. 2020). In the current study, the PFI was applied to each of the predictors. The model reduction, i.e. the PFI score of a predictor, was calculated as the ratio between the TPR of the non-permuted model and the TPR of the permuted model. FPI scores were evaluated for each of the 100 ANN replicates thus allowing the evaluation of a statistical distribution of the predictors' importance. Also, we grouped the 167 predictors into 5 groups (Road, Hydro, Geo, Terrain, Seismic; see Table 1 for more details) and we investigated how the network performance varies by running the classification 20 times with each of all the possible different combinations of the five groups. Finally, the susceptibility map was verified by means of a comparison with the Check dataset, and for each Italian administrative region, an additional check TPR was calculated, as well as the percentage of territory classified as unstable. Table 4 shows the average values and standard deviation of the TPR, TNR, and AUC general performance indicators obtained through the 100 ANN replicates. The results are reported for the three types of datasets we considered, namely training, validation, and test. Furthermore, we also report the values obtained for the dataset composed of the sum of the three subsets (All). The results in the table show high performances for the three indicators considered. The average values for the three datasets are also comparable, demonstrating that the approach followed is able to limit any evident overfitting effect and the consequent loss of generality in the slope unit classification phase. Very limited values of the standard deviations also demonstrate the robustness of the method, which is able to obtain comparable  performances regardless of the specific datasets used in each of the 100 replicates.

Results
Considering the comparability of the performances obtained on the three training, validation, and test datasets, for the following results, it was considered appropriate to report those obtained on the overall dataset composed by the three. Figure 5 a shows (in grey) the ROC obtained for each of the 100 replicates of the ANN training. The average ROC is shown in red. In this study, AUC = 0.91 has been reached on average, with a standard deviation smaller than 0.02 (see Fig. 5a and Table 4). Besides the best classification threshold that resulted in being about 0.5, in Fig. 5a, the TP and FP rates related to other eight different thresholds (from 0.1 to 0.9) are indicated by the means of black circles in order to allow a deeper and direct interpretation of the EQtLs susceptibility probability value that the model associates to each SU. As an example, by choosing a threshold value of 0.8 a very low FPR (about 0.06) is obtained, meaning that only a very limited fraction of the stable SU would be wrongly classified as unstable. As a result, those SUs that have been classified with a probability higher than 0.8 to be susceptible to EQtLs, are statistically very significantly likely to have actually experienced landslides/be true positives. Figure 5 b shows, for the 100 ANN replicates, the values of the TPR parameter according to the TNR parameter. Mean and standard deviation ranges are also reported for both TPR and TNR. On average the classification has a very similar success rate for both classes (about 0.84) with a small standard deviation (0.02), demonstrating that the classification is carried out with the same accuracy for both classes. The low value of the standard deviation and the absence of correlation between the values of TPR and TNR also make it possible to assert that the results obtained are robust with regard to the statistical representativeness of the samples considered and the absence of bias introduced.

Susceptibility mapping
After every training replicate, the ANN was applied to all the slope units of Italy and 100 susceptibility values for each SU have been generated. The mean susceptibility of each SU, and its standard deviation, after 100 replicates has been considered to produce the EQtLs susceptibility map of Italy (Fig. 6a).
In the EQtLs susceptibility map of Italy, flat lowland areas have been taken out from the classification and resulted in grey colour. Orange to red areas represent moderately to highly susceptible slope units (probability > 0.5), while green to blue areas have been classified as stable. Susceptible areas are frequent in the northeastern part of Italy and along an NW-SE-oriented longitudinal belt that corresponds to the Apennine mountain chain. In particular, red areas are located in correspondence with the epicentral area of historically strong earthquakes and a moderate density of unstable slope units is present in the Calabria region, the most southern region of the Italian peninsula. Conversely, most of the western side of the peninsula and of the alpine region, in the north, are lowly susceptible to being affected by EQtLs. Also, the south-east and the two main Italian islands, Sicily and Sardinia, are widely blue-coloured. The standard deviation of the resulting classification (Fig. 6b), associated with the mean susceptibility of every SU, is very low (< 0.1) in correspondence with the high susceptibility SUs in central Italy and in the north-east, as well as for most of the highly stable areas. In general, the standard deviation of the susceptibility is low (0.1-0.18) for the overall Italian territory. Higher values are present in limited spotted locations and more concentrated in the Calabria region.
A point of novelty of this study is represented by the comparison of the landslide Susceptibility map of Italy with an EQtLs dataset that was not used to train the network.
The overlapping between checking landslides and unstable SUs has been evaluated to verify the correctness of the susceptibility map. In order to make the checking process reliable, a radius sized as the associated error has been taken into account around the less precisely georeferenced landslides. When more than half of the area of the resulting circle overlapped with unstable slope units, that landslide was considered a true positive (TP). Conversely, when the overlap was limited to less than half of the circle area, landslides were considered false negatives (FN). When some parts of the uncertainty circles included areas with no classification (e.g. lowlands or sea), only the portion overlapping with classified slope units was considered. Consequently, the checking TPR (C-TPR) has been calculated for every Italian region. On the basis of the susceptibility map, also the regional percentage of unstable territory has been computed. As a result, in most of the Italian regions, the number of TP was higher than FN, although not all the regions counted the same number of landslides from the checking dataset. In cases of regions with only a few (< 10) checking landslides, C-TPR generally reached very small values. Conversely, Friuli, Veneto, Emilia-Romagna, Tuscany, Abruzzo, Molise, Campania, and Basilicata show very good performances (C-TPR ≥ 70%) and a high number of checking landslides (> 14), which makes the resulting statistics quite reliable. In these regions, the percentage of unstable territory varies from around 20-40% to more than 60%. Contextually, Lombardy, Latium, Sicily, and Calabria show low to very low C-TPR despite the good number of checking samples. In Calabria, 36% of the regional extent has been classified as unstable, while in the other three regions, the unstable territory is < 20%. Nevertheless, considering the low C-TPR values achieved, these percentages might have been probably underestimated.

Predictors' importance
PFI provided an interesting analysis of the importance that the single predictor had in order to achieve the final classification.
In Fig. 8, it can be seen that the ANN mainly relies on five predictors, all listed among Geothematic and Seismic ones. In particular, soil type (code 122), distance from seismogenic sources (123), lithology (31), distance from active faults (125), and geomorphon (10) have the highest PFI score, respectively. The first terrain predictor in order of importance is represented by the mean tangential curvature of a slope unit (code 144). Its importance, however, varied significantly among the 100 replicates. Following, all the other predictors, such as other terrain predictors and the road-related ones, account for very little contribution to the classification and the associated PFI standard deviation is small. On the basis of the EQtLs susceptibility map, the a posteriori distribution of the classes of GEO predictors among the unstable slope units has been analysed at the national level in Fig. 9.
To make the chart clearer, only soil types with unstable slope units higher than 10% have been reported. Concerning soil types, slope units mainly covered by Dystric Cambisol resulted highly susceptible to EQtLs, and 75% of them have been classified as unstable, although they are not numerous (< 5000 in the whole national territory). In the WRB system, "Dystric" indicates a soil with base saturation of less than 50% at a given depth, and Dystric Cambisol is located in small parts of central and south Apennine, in seismically very active areas, which have been historically hit by strong earthquakes. Further, more than 60% of the slope units composed of Rendzic Leptosol have been classified as unstable. Rendzic Leptosol is described as very shallow soils immediately overlying highly calcareous material and is quite frequent in Italy, particularly in central and south Apennine as well as in Friuli and Veneto regions. Almost 50% of slope units characterised by the main presence of Chromi-calcaric Luvisol have been classified as unstable. Chromi-calcaric Luvisol is defined by the WRB system as a reddish calcareous with a marked textural differentiation whose surface horizon was depleted of clay, which accumulated more in-depth. And, according to the pedological map of Europe, it is very rare in Italy. Finally, almost 40% of slope units were dominated by Lithic Leptosol, a very shallow soil with continuous hard rock within 10 cm from the soil surface (Table 3), resulted susceptible to EQtLs. In Italy, its occurrence is limited to central Apennine, between Latium and Abruzzo, and in Sicily island. Concerning lithology, 75% of slope units are mainly constituted by chaotic sedimentary complexes, and 50% of those composed by marls have been classified as unstable. The first lithology is composed of turbiditic sandstones and clays, locally with evaporites and limestones. It is mainly spread in central Italy, along the eastern side of the Apennine chain. The second type of rock is spread in central Italy and the north-west. Successively, 25-40% of arenaceous and limestone slope units resulted susceptible to EQtLs. Arenaceous formations crop out all over the Italian territory, mainly in mountain areas, while Limestones are spread in seismic regions of central Italy, such as Umbria and Abruzzo, in the southern part of the Alps, and along the coasts of south Italy. Finally, metamorphic rocks, mainly granitoid gneiss, whose almost 20% of slope units are considered unstable, crop out only in the northern part of the Alps and in small parts of Calabria and Sicily regions. Concerning the slope morphology, valley, and concave slope units interestingly resulted to be relatively more unstable than slope units located in other parts of the slope. In detail, the 25-35% of hollow, valley, and depression slope units have been classified as unstable against the 15-20% of summit, ridge, spur, and slope classes. Finally, slope units which are linked with flat areas, such as flat, shoulder and footslope, are generally stable.
In this paragraph, the analysis of how the classification performance changes varying the combination of groups of predictors used by the ANN is also provided. Predictors have been grouped as Terrain, Seismic, Geo (i.e. Geothematic), Hydrological and Roads (Anthropic) as described in "Material and methods." All possible combinations made up of a variable number of groups have been taken into account (one group at a time up to all five groups together). For each of the possible combinations among these groups, the ANN has been run 20 times, and the related AUCs have been calculated. Figure 10 shows the box plot of the AUC values distribution among the 20 replicates and for all the possible combinations of predictor groups.
Combinations are ordered by the medians of the AUC distributions. The background colour varies according to the quartiles of the distribution of the median AUCs calculated over the 20 replicas per combination. The median quartiles are at AUC values of 0.84, 0.88 and 0.89. Lower performances (AUC < 0.84) are generally achieved with only one or two groups, or with 3-groups combinations that contain Hydrology and Roads but not Geo. Good performances (0.84 < AUC < 0.89) are achieved with all the 2-groups combinations that include Geo. In this regard, Geo+Seismic performs the best. Also, combinations with three or four groups achieve good performances. Finally, those combinations of predictors groups whose AUC is entirely included in the dark red band and greater than the third quartile (AUC > 0.89) can be considered the best-performing ones. Among these, two 3-groups' combinations are listed.
1. Geo + Seismic + Road 2. Geo + Seismic + Terrain whereas the three 4-groups' combinations are: 1. Geo + Seismic + Terrain + Hydrology 2. Geo + Seismic + Hydrology + Road 3. Geo + Seismic + Terrain + Road and only one combination exists with all five groups. From the analysis of the best-performing combinations, it is clear as Geo and Seismic predictors must be both considered in order to achieve a median AUC higher than 0.89, and that at least another group is also needed. The importance of Geo (i.e. lithology, soil type and geomorphon of slope units) and Seismic (i.e. distances from active faults and seismogenic sources) predictors was previously indicated also by the PFI analysis. Nevertheless, what and how many predictor groups are needed besides Geo and Seismic was not straightforward. Related to this, on the basis of the interquartile range and the median of AUC values, Geo+Seismic+Terrain+Road and Geo+Seismic+Terrain seem to perform slightly better than all the other combinations. Figure 11 represents a heatmap of the mean AUC value obtained by adding one of the five groups of predictors to each of all their possible combinations.
Each row contains one of the possible combinations and is sorted from top to bottom by the increasing number of groups. In each column, one of the five groups is present. The mean AUC obtained after 20 ANN replicates considering the combination in the row and the addition of the group in the column is reported in each cell of the heatmap. "Null" row and column respectively indicate that none of the possible combinations has been considered and that no groups have been added. Figure 10 confirms what has been previously seen in Fig. 10 that the higher the number of groups within a combination, the higher the performance. Nevertheless, not all the groups have the same effect. When the classification has been carried out taking only one group at a time (first row on the top), Terrain and Geo performed the best, with mean AUC = 0.84, and significantly better than Seismic (mean AUC = 0.79) although some of the Seismic features resulted among the most important in the full model PFI analysis. Nevertheless, Terrain + Seismic reaches AUC > 0.9 only when Geo is added while, conversely, Geo + Seismic reaches AUC > 0.9 also with Roads appearing that, when combined with Geo, Seismic provides a bigger contribution than Terrain. This led to infer that Terrain and Geo groups might bring partially overlapping information and that those brought by Seismic better combine with Geo than with Terrain features. In general, when the Geo group is added to whatever combination (second column from the left in Fig. 11), the mean AUC reaches 0.9 in seven cases and it never goes below 0.8. This means that the Geo predictors have a high importance for the ANN and their presence ensures very good performances, whatever another group is added to the combination. Similarly, Seismic predictors allow reaching mean AUC ≥ 0.9 when added to six different combinations. Further, when they are present, performance decreases below 0.8 only in one case and the combination Geo + Seismic achieves AUC = 0.89. Conversely to Geo and Seismic, only four combinations that include Hydrological predictors allow achieving a mean AUC of at least 0.9 and, in all these cases, Geo is present. Also, two combinations that include Hydrological predictors do not reach AUC = 0.8. Finally, five combinations containing Roads and five combinations containing Terrain reach a mean AUC = 0.9. This means that the probability to reach very good performance by a combination that contains Road-related predictors is the same as a combination that includes terrain predictors. However, when Roads is added to Geo and Seismic, AUC arrives at 0.91 and vice-versa, when Terrain is added to Geo and Seismic, AUC averagely arrives at 0.92. Adding Roads to Geo + Seismic + Terrain brings a contribution lower than 0.01 while, adding Hydrology, performance decreases to 0.90. Besides, 0.92 is the highest mean AUC reached by the ANN and is due to the main contribution of Geo and, secondly, Seismic information. Terrain predictors would have a much higher importance, when grouped, than that resulting from the single predictors' analysis. But its information might be partially provided also by Geo predictors and, when combined with other groups, it accounts for slope units variability less than Geo group, ending to provide only + 0.03 to the combination Geo + Seismic. From this analysis, the key role of Geo and Seismic predictors is confirmed and emphasised. Also, a significant contribution of Terrain has been proven. At the same time, the non-significance of distance to rivers as a predictor for EQtLs susceptibility resulted and a not ignorable contribution to improving the classification performance given by the presence of roads. Finally, concerning what can be selected as the most performing combination among all the possible and tested ones, it should be noted that the differences between the mean AUC values for the three best median AUC combinations, which are Geo + Seismic + Terrain, Geo + Seismic + Terrain + Road, and Geo+Seismic+Terrain+Road+Hydro, are not statistically significant (p-value = 0.86 with one-way ANOVA test).

Discussions
The sections below are meant to provide the reader with an overview of the strengths and potential weaknesses of the modeling protocol we implemented, these discussed both from the data as well as the modeling strategy perspectives.

ANN performance overview
The ANN performance was very good. In detail, after 100 replicates mean AUC was 0.91 and the associated standard deviation was 0.01. Considering that both positive and negative samples (i.e. slope units with and without landslides) within training, test, and validation datasets changed at every replicate, the very low standard deviation is an excellent result, which demonstrates the solid stability of the network. Also, the ability to distinguish between the two classes was high: averagely, TPR, namely the ability to correctly classify unstable slope units, was 0.85 while TNR, proficiency in The dark bars represent the regional estimated percentage of unstable territory according to the susceptibility map Fig. 8 Resulting scheme of the Permutation Feature Importance analysis. Predictors codes are provided in Table 1 classifying stable slope units, was 0.84. Both metrics show a standard deviation lower than 0.02 after 100 replicates confirming the robustness of the classification. In particular, the classification error plot shows low standard deviations especially for those SUs classified as extremely stable (mean susceptibility < 0.25) or unstable (mean susceptibility > 0.75), giving rise to high reliability of the final susceptibility model. These outputs fulfil the aim of the work to perform a Fig. 9 Distribution of slope units among the three geothemathic variables classes: bars refer to the percentage of unstable slope units out of the total number of slope units, per class; diamonds indicate the total number of slope units, per class. a Refer to soil type classes. To make the chart clearer, only soil types with unstable slope units higher than 10% have been reported. b Refer to lithology classes and c to geomorphon classes robust susceptibility analysis of earthquake-triggered landslides at the national scale which, being trained on landslides distributed over more than one century and over the whole Italian territory, could serve as a basis to prioritise funds for remedial interventions at national to regional levels.

EQtLs susceptibility patterns
The EQtLs susceptibility map of Italy obtained by the means of the neural network approach was compared with a landslide distribution map of Italy derived from the IFFI inventory, which reports landslides whose genesis is linked to rainfall, earthquake, snowmelt, and anthropic effects.
The comparison reveals an interesting output which regards the main distribution of earthquake-induced landslide all along a more internal portion of the Apennine Chain backbone (Fig. 12). As a result, the eastern coastal zone is less predisposed to landslide triggering due to earthquakes.
On the north, along the Alps Chain, the highest susceptibility zone corresponds to the eastern area, namely parts of Veneto and Friuli regions, where seismogenic sources are more concentrated.
The average dimensions of the chosen mapping unit, i.e. 0.7-km 2 slope units, provide a detailed level of spatial resolution to the susceptibility map, which represents an accurate model when observed at the regional scale and clearly identifies what are the more susceptible areas with respect to the more stable ones. In those regions where the a-posteriori model check reaches high performance (C-TPR > 70% in Fig. 7), such as Veneto, Tuscany, Friuli, Abruzzo, Emilia-Romagna, Molise, Campania, and Basilicata regions, the produced EQtLs susceptibility map can be taken as a reliable instrument to drive the decision makers toward appropriate funding management, i.e. in order to provide priority lists of local interventions. The comparison between the overall landslides density map and the here presented EQtLs susceptibility map clearly indicates that areas highly susceptible to earthquake-triggered landslides might be not taken into account in the frame of landslide mitigation National funds since not necessarily exposed to a high generic landslide hazard, e.g. rainfall-induced. Contextually, in these areas, the likely dedicated funds for earthquake risk mitigation might tend to be used primarily for building reinforcement, keeping on ignoring the significant slope stability matter demonstrated with this study.

Predictors' role
From predictors importance analyses, Geothematic and Seismic resulted in being the most important predictors, as largely supported in the literature by other susceptibility  Tanyaş et al. 2019;Loche et al. 2022b). From the PFI analysis, it resulted that soil type, distance from seismogenic sources, lithology, distance from active faults, and geomorphon are the most important predictors for the network and, as consequence, for the good result of the classification between stable and unstable slope units. As one can expect from an application on earthquaketriggered landslides, distance from seismogenic sources and distance from active faults played a key role in the classification, demonstrating to well represent the slope units variability due to seismic predictors. In this regard, the seismogenic source represents the portion of a fault that is more likely to enucleate a M w > 5.5 earthquake (Basili et al. 2008;DISS-Working-Group 2018). Nevertheless, landslides can occur even along the tip portions of a fault, after M w < 5.5 seismic events or in correspondence with secondary segments; therefore, distances from both source and fault line have been considered. Taking into account seismogenic sources and active faults, the susceptibility analysis presented in this study resulted in an inclusive model that is not bound to some specific seismic events and can be applied to the whole National territory accounting for more local variability than that provided by the national PGA. Soil type resulted to be the most important predictor from the PFI analysis. Statistically speaking, this may be partially due to the fact that slope units have been characterised on the basis of 91 different soil types, giving rise to a high, detailed, pedological variability. This represents an impressive quantity of data for the ANN to take useful information from in order to perform the classification. Related to this, lithology and geomorphon, which only count 21 and 10 classes respectively, might have provided lesser, albeit meaningful, information, that resulted as third and fifth more important predictors, respectively, among all the considered ones. In percentage, the most unstable soil categories resulted in poorly developed pedotypes, generally thin, and derived from the alteration of rocky or highly calcareous bedrocks. This result is in line with what resulted from the analysis of the instability percentage per lithology class, which shows as calcareous and arenaceous formations, besides Fig. 11 Heatmap of the mean AUC values after 20 replicates for all the combinations of predictors groups. Combinations are obtained by adding the group in the column to the combination in the row. Combinations in rows are sorted from top to bottom by the increasing number of groups. In each column, one of the five groups is present. "Null" row and "null" column indicate that none of the possible combinations has been considered and that no groups have been added, respectively clayey and marly lithologies, are largely present within the unstable slope units. These results reflect the high abundance of disrupted failures that affect rock masses during an earthquake, like rock fall, which also represents the most numerous landslide type in the Input dataset (Fig. 2). It is particularly relevant that slope units with prevalent Chromi-calcaric Luvisol, although they are very rare in Italy, are in percentage rather unstable. They are located in the Veneto region, within a restricted area quite close to a few landslides that occurred as a consequence of an earthquake that struck the region in 1936 and had an epicentre 60 km away. In this case, a clustering effect cannot be ruled out: it may have been the earthquake, with consequent EQtLs, to occur in correspondence with areas with Chromi-calcaric Luvisol rather than the presence of this soil favouring the trigger of EQtLs.
Concerning lithology, 75% of slope units are mainly constituted by chaotic sedimentary complexes, and 50% of those composed by marls have been classified as unstable. The first lithology is composed of turbiditic sandstones and clays, locally with evaporites and limestones. It is mainly spread in central Italy, along the eastern side of the Apennine chain. The second type of rock is spread in central Italy and the north-west. Successively, 25-40% of arenaceous and limestone slope units resulted susceptible to EQtLs. Arenaceous formations crop out all over the Italian territory, mainly in mountain areas, while Limestones are spread in seismic regions of central Italy, such as Umbria and Abruzzo, in the southern part of the Alps, and along the coasts of south Italy. Moreover, an exploratory analysis of the data highlights chaotic sedimentary complexes as the lithotype is mostly prone to slope failures. This may be due to the unstructured bond between constituents, which under seismic stresses tend to generate mass movements more than in other more cohesive/massive lithologies.
Concerning geomorphon, its importance may be linked to several aspects. Intuitively, slope morphology is strictly correlated with the probability of landslide occurrence. Related to this, geomorphon classes could properly represent the most important morphological features of a slope, accounting for the contribution given by the terrain predictors to the ANN, which did not result as important as one could expect. In particular, our analysis highlights that hollow, valley, and depression slope units have been classified as more frequently unstable than summit and ridge areas. This aspect can be partially linked with the presence of minor scarps along the slope (due to strata bedding, faults and fractures, terraces, and other natural reasons) and with the presence of roads in the lower part of the valley sides, as it occurs in many mountainous regions of Italy. Martino et al. (2019) and Tanyaş et al. (2022) found out that the presence of road cuts at the bottom of deeply incised V-shaped valleys played a conditioning role in the spatial distribution of EQtLs triggered in 2016 in Central Italy higher than road cuts located elsewhere. In this regard, from the analysis of the predictor group combination, the contribution of the distance from roads to the classification performance resulted to be not negligible. And cases of EQtLs mostly occurred along slopes modified by road cuts are documented in the literature (Keefer et al. 2006;Delgado et al. 2015). Therefore, despite the role of roadcuts in favouring EQtLs occurrence is probably more appreciable in applications in small study areas, this study also contributes to the ongoing debate drawing attention to the potential importance of this predictor.

Opposing arguments: Validation routine through the Check data
The Italian EQtLs susceptibility map, albeit resulting from a robust iterative training-validation-test procedure, shows heterogeneous results between different regions from the comparison with the Check EQtLs dataset, which was not used to train the ANN. This can be partially due to a low number of check landslides in regions like Piedimont, Aosta Valley, Liguria, and Apulia (Fig. 7). Nevertheless, Calabria and Sicily show low C-TPR in spite of a high number of checking landslides. A reason for this could be that, as shown in zoom 1 in Fig. 12, the Input landslides from which the ANN was trained are concentrated in the area of the Strait among these two regions, while the checking EQtLs are more spread over the regional territory. This may have led to a too-low-generalised training of the network. Further, the exact location of seismogenic sources in Calabria is an argument of debate in the scientific community. Considering the importance assumed by the Seismic predictors in our analysis, their potential location's uncertainty can significantly affect the accuracy of our model.
Ultimately, an ideal model should holistically include all the factors which are known to control the slope failure mechanisms. Therefore, one should theoretically build an analogous ANN by featuring soil moisture characteristics or at least a proxy of this property, usually expressed in terms of precipitation. However, our training dataset spans multiple earthquakes, even historical ones. For this reason, we did not include any information related to soil moisture and its preparatory effect for slope failures under seismic shaking. It is our intention to fulfil this requirement for a subset of the data we considered here. Specifically, we are already Fig. 12 a Map of landslide density derived by the IFFI inventory. Each pixel is 5 × 5 km. b Earthquake-triggered susceptibility map of Italy produced in this study. c Map of active faults in Italy. Zooms of the maps in a and in b are shown in the upper and the lower bands, respectively. Circles represent the landslides of the Input dataset used in the ANN training process. Triangles represent the landslides of the checking dataset used for an ex-post evaluation of the susceptibility map ◂ working on the earthquake sequence that affected central Italy in 2016. This case would offer a perfect stage to include soil moisture in the analyses.

Concluding remarks
This study represents the first example of susceptibility analysis to Earthquake-triggered Landslides built at the national scale in Italy by using an ANN approach. At this aim, we exploited the CEDIT catalogue, and we implemented an ANN at the Slope Unit scale featuring predictors that take into account predisposing factors of morphological, lithological, and seismotectonic characteristics of the Italian territory. To train the ANN, a sub-dataset made of 1545 EQtLs, most accurate in terms of geolocalisation, and related to strong earthquakes from 1908 to date, was extracted from CEDIT.
The ANN was optimised and trained for the classification of Slope Units in terms of susceptibility to earthquake-triggered landslides, based on 167 predictors. The performances of the ANN have been evaluated by carrying out 100 training on independent datasets to assess its robustness. These produced a remarkable performance distinguishing between the stable and unstable SU classes. We also tested our model via an external dataset composed of 465 EQtLs, which were not featured in the model-building phase due to a lower geolocation accuracy. This expost verification confirmed the overall suitability of our ANN analytical protocol.
The analysis shows that a large portion of the Italian national territory is highly prone to earthquake-triggered landslides. This is especially the case throughout the Apennine arc, where high susceptibility values are associated with more than 50% of Abruzzo, Marche, Molise, and Umbria and ∼ 25% of Tuscany, Emilia-Romagna, Campania, and Basilicata regions. Furthermore, high susceptibility values are associated with approximately 25% of the territory of the Friuli region, in the eastern sector of the Alpine arc. Conversely, north-western regions, Sicily, Sardinia, and most of Latium and Apulia appear to be quite stable with minor percentages of the territory characterised by susceptible slopes.
The 0.7-km 2 SU mean size is extremely detailed, even for a regional or provincial scale assessment. For instance, our susceptibility map would be a perfect addition to Civil Protection in order to design emergency plans at the provincial level. Or it could be used at the level of each region to plan long-term risk mitigation actions and allocate funds. Our work demonstrates that if funds were allocated only based on rainfall-induced landslide distribution, large areas highly susceptible to EQtLs would be completely ignored by mitigation plans.
Nevertheless, the resolution of these mapping units is still far from the requirements for planning purposes, or for seismic microzonation studies, at a municipal scale.
As regards future improvements we envision for this study, the main extension consists of scaling down the model to a greater spatial resolution focusing on a specific datarich sub-region, of particular interest. Here, we envision to run also a suite of physically based models, besides the ANN-base procedure and compare the two outputs in order to investigate potential differences, both in terms of strength and weaknesses, between each modeling routine. Table 2 List of the predictors assigned to each slope unit. Codes reported in "Predictors code and description" have been used to represent the results of the permutation feature importance. Predictors have been grouped as indicated in "Group" to perform the combination-group analysis. Terrain characteristics were calculated from a 20-m Digital Elevation Model (DEM)

Group Code Description
Geothematic predictors Geomorphon (10 classes) 10 Class of "geomorphon" which covers most of the area of the Slope Unit (SU) Lithology (21 classes) 31 Lithology covering most of the SU area; it is taken from the geological map of Italy at a scale of 1:50.000 or 1:100.000 Soil type (91 classes) 31 Type of soil that covers most of the area of the SU; it is taken from the ecopedological map of Italy at a scale of 1:250.000 Seismic predictors Distance to seismic features 123 The average distance of a SU from the nearest seismogenic source 124 The standard deviation of the distance to the nearest seismogenic source 125 The average distance of a SU from the nearest active fault line (capable or not) 126 The standard deviation of the distance of a SU from the nearest active fault line (capable or not)

Group Code Description
Anthropic predictors Distance to roads 127 Count of the pixels of a SU covered by any buffer of distance from roads. The buffer ranges are 10, 50, and 100 m 128 Sum of the buffer values of the pixels of a SU covered by any buffer of distance from roads 129 The buffer of distance from a road that takes up most of the area of the SU 130 The maximum value of the buffer of distance from a road within a SU 131 The average value of the buffer of distance from a road within a SU 132 The minimum value of the buffer of distance from a road within a SU 133 Range (max-min) of distance values from roads included in a SU 134 Count of the pixels of a SU that fall within 5 m of distance from roads Hydrological predictors Distance to watercourses 135 Count of the pixels of a SU covered by any buffer of distance from rivers. The buffer ranges are 10, 50, and 100 m 136 Sum of the buffer values of the pixels of a SU covered by any buffer of distance from rivers 137 The buffer of distance from a river that takes up most of the area of the SU 138 The maximum value of the buffer of distance from a river within a SU 139 The average value of the buffer of distance from a river within a SU 140 The minimum value of the buffer of distance from a river within a SU 141 Range (max-min) of distance values from rivers included in a SU 142 Count of the pixels of a SU that fall within 5 m of distance from rivers Terrain predictors Curvature 143 The average Tangential Curvature of a SU 144 The standard deviation of the Tangential Curvature of a SU 145 The average Profile Curvature of a SU 146 The standard deviation of the Profile Curvature of a SU 147 The average Plan Curvature of a SU 148 The standard deviation of the Plan Curvature of a SU 149 The average Longitudinal Curvature of a SU 150 The standard deviation of the Longitudinal Curvature of a SU 151 The average General Curvature of a SU 152 The standard deviation of the General Curvature of a SU Elevation 153 The average Elevation of a SU 154 The standard deviation of the Elevation of a SU Exposure 155 The average Exposure of a SU from north to south 156 The standard deviation of the Exposure of a SU, from north to south 157 The average Exposure of a SU from east to west 158 The standard deviation of the Exposure of a SU, from east to west Slope 159 The average Slope of a SU

160
The standard deviation of the Slope of a SU Topographic Wetness Index 161 The average Topographic Wetness Index of a SU 162 The standard deviation of the Topographic Wetness Index of a SU Topographic Position Index 163 The average Topographic Position Index of a SU 164 The standard deviation of the Topographic Position Index of a SU Topographic Ruggedness Index 165 The average Topographic Ruggedness Index of a SU 166 The standard deviation of the Topographic Ruggedness Index of a SU Area 167 SU areal extent  Table 4 Description of the reference soil groups compared in the classes of the categorical predictor "Soil type"

Soil type WRB code Description
Cambisols CM Soils with little or no profile differentiation-Moderately developed. Cambisols are developed in medium and finetextured materials derived from a wide range of rocks, mostly in alluvial, colluvial, and aeolian deposits Regosols RG Soils with little or no profile differentiation-No significant profile development. Regosols are developed in unconsolidated materials. Regosols are extensive in eroding lands, in particular in arid and semi-arid areas and in mountain regions Luvisols LV Soils with clay-enriched subsoil-High-activity clays, high base status. The main characteristic is an argic horizon, a subsurface zone with higher clay content than the material above it. This typically arises as clay is washed downward by water and accumulates at greater depth Andosols AN Soils distinguished by Fe/Al chemistry-Allophanes or Al-humus complexes. Andosols are generally quite young soils found in volcanic areas formed in volcanic tephra. Andosols are usually defined as soils containing high proportions of glass and amorphous colloidal materials, including allophane, imogolite and ferrihydrite Calcisols CL Accumulation of moderately soluble salts or non-saline substances-Accumulation of secondary carbonates. Calcisols are developed in mostly alluvial, colluvial, and aeolian deposits of base-rich weathering material. They are found on level to hilly land in arid and semi-arid regions. The natural vegetation is sparse and dominated by xerophytic shrubs and trees and/or ephemeral grasses Leptosols LP Soils with limitations to root growth-Thin or with many coarse fragments. Leptosols are very shallow soils over hard rock or a deeper soil that is extremely gravelly and/or stony. Leptosols can be found on hard rocks or where erosion has kept pace with soil formation or removed the top of the soil. The very shallow, less than 10 cm deep, Lithic Leptosols in mountain regions are the most extensive Leptosols on Earth Fluvisols FL Soils with little or no profile differentiation-Stratified fluviatile, marine, and lacustrine sediments. Fluvisols are found on alluvial plains, river fans, valleys, and tidal marshes on all continents and in all climate zones. Under natural conditions, periodical flooding is fairly common. The soils have a clear evidence of stratification. Soil horizons are weakly developed, but a distinct topsoil horizon may be present