This section will briefly introduce the study area and the landslide inventory we used. Subsequently, we will describe the spatial partition and the covariate set, together with the model we selected, referring to the articles where an extensive mathematical formulation is provided.

## 2.1. Study area and landslide inventory

The study area where we tested our modeling strategy is located within the *Cinque Terre* National Park, Italy. Such a name derives from the five epochal villages of *Monterosso al Mare*, *Vernazza*, *Corniglia*, *Manarola* and *Riomaggiore*. These hamlets are placed along a coastal stretch, which represents a worldwide known tourist attraction. Because of its environmental, cultural, and historical heritage, the *Cinque Terre* territory was listed as a World Heritage Site by UNESCO in 1997 and was declared National Park in 1999. Being most of the *Liguria* region characterized by a rugged morphology and by the absence of large flat areas adapted for cultivation, slope soil covers, over the century, were intensively reworked to build terraces sustained by dry-stone walls (Brandolini, 2017). Terraced slopes cover a wide portion of the study area (about 15 km2); even if, about 67% of the terraced areas were abandoned (Raso et al., 2020). Agricultural terraces, if correctly preserved, assume a positive role on soil conservation by reducing runoff velocity and soil erosion (Moreno-de-las-Heras et al., 2019). Conversely, when the terraces are abandoned, land degradation issues (i.e., gully erosion, terrace failure, mass movement, piping, hydrological connectivity) arise (Tarolli et al., 2014; Arnáez et al., 2015; Di Napoli et al., 2021). In addition, when intense rainfall occurs, abandoned terraced are particularly prone to shallow mass movements (Brandolini et al., 2018; Cevasco et al., 2015). High-intensity precipitation mobilizes huge amounts of materials having a considerable impact on the solid discharge, flow and energy of streams. These dynamics often cause flow-like movements, debris floods and flash floods, which are increasingly affecting, in particular, the coastal settlements (Zingaro et al., 2019; Brandolini et al., 2012).

These phenomena represent a serious threat to human settlements, inhabitants and trail users, as dramatically the area experienced after the intense rainstorm that hit the *Monterosso* and *Vernazza* areas on 25 October 2011 (Cevasco et al., 2015; Rinaldi et al., 2016), which triggered hundreds of shallow landslides as well as destructive debris floods (Cevasco et al., 2014). Figure 1 shows an overview of the study area and the landslides triggered by the convective storm occurred on October 25th 2011. During that day, up to 382 mm of rain were discharged in few hours, as recorded at the weather station of *Monterosso* (Cevasco et al., 2015).

Landslide inventory was redacted from detailed field surveys and analysis of high-resolution aerial images (ground resolution from 3 to 50 cm, according with the altitude) taken a few days after the catastrophic event. The photointerpretation analysis, executed from georeferenced orthophotos provided by *Liguria* Regional Administration, was validated by an intensive field surveys carried out from November 2011 to March 2012 (Cevasco et al., 2013). The mapped mass movements were classified according to Cruden and Varnes (1996) based on the prevalent type of movement and material. The main landslide types recognized can be associated to debris flow, debris avalanche and debris slide, involving colluvial and anthropically reworked deposits overlaying the fractured bedrock. The average landslide density was about 65 landslides/km2 and landslide areal extent ranged between a few tens up to a few thousands of square meters. Lastly, the event inventory featured 695 mass movements.

## 2.2. Mapping Units

To model landslide intensity we chose a hierarchical structure. The high resolution mapping unit corresponds to Grid-Cells (GCs, Malamud et al., 2004). These are hierarchically combined with the coarser Slope Units (SUs, Carrara et al., 1995), at which level we computed the Latent Spatial Effect and we aggregated the intensity estimates (see (Lombardo et al., 2019a). Specifically, we selected a 20m resolution GC partition whereas we computed the SUs by using the *r.slopeunits* software (Alvioli et al., 2016). We parameterized *r.slopeunits* with a circular variance of 0.4, a minimum SU area of 12500 m2 and a flow accumulation threshold of 100000 m2. This operation returned 171 SUs. Circular variance is the parameter that *r.slopeunits* uses to control the rigid or flexible aspect criterion for the SU generation. For instance, a value close to zero indicates a very limited aspect variation allowed inside a slope unit. Taking this figuratively to the extreme, the most rigid choice would end up creating one slope unit polygon for each pixel in the aspect raster. On the contrary, a circular variance close to one would make *r.slopeunits* too generous in the SU delineation. Taking this figuratively to the extreme again, the most flexible choice would produce a single SU polygon for the whole study area because it would allow for large aspect variations. The following two parameters are dependent of the concept that catchments and half-catchments are fractal physiographic entities (La Barbera and Rosso, 1989), which means that one can apply the same concept to extract them at different geographic scales (or hydrological order; Blösch and Sivapalan, 1995). The minimum SU area is in fact the parameter that controls the convergence of the polygonal partition to a specific reference planimetric area. As for the accumulation threshold, this parameter controls the starting extent of the half-catchment delineation.

## 2.3. Covariate set

The morphometric covariates we chose to build our intensity model were derived from a 5 m Digital Elevation Model (DEM) accessed from the geo-portal of the Ligurian region. This DEM has been later resampled at 20 m resolution to match the squared lattice we defined. We computed the Euclidean distance from each GC to the nearest road or trail.

We also used the thematic properties described in (Di Napoli et al., 2021). As a result, our covariate set featured: *i*) Elevation; *ii*) Slope Steepness; *iii*) Eastness; *iv*) Northness; *v*) Planar and *vi*) Profile Curvatures; *vii*) Relative Slope Position; *viii*) Topographic Wetness Index; *ix*) Distance to road or trail; *x*) Land Use; *xi*) Terraced slope status; *xii*) Geology.

## 2.4. Landslide Intensity Modeling

By counting the distribution of slope failures per mapping unit, we can model the resulting data as a Point Process. This is possible under the assumption of spatial continuity, which is something a raster structure naturally offers. As a result, we can define a Poisson Point Process as:

$$N \left(A\right) \sim Poisson {\int }_{A}^{}\lambda \left(s\right) ds$$

1

where *N(A)* is the number of expected landslides within the study area *A*, the selected sector of the *Cinque Terre* National Park in this case, *λ* is the intensity assumed to be ≥ 0, and *s* is each of the GC within the target area. Notably, here *λ* is assumed to behave according to a Poisson probability distribution. This framework can be further extended conveniently expressing the intensity in logarithmic scale. This procedure then gives rise to what in statistics is referred to as a Log-Gaussian Cox Process (LGCP). A LGCP is particularly convenient because a Log-Gaussian structure allows one to incorporate any type of linear and nonlinear covariate effects. In our case, this leads to denoting our model as follows:

$$\text{log}\left\{\lambda \left(s\right)\right\}\sim Gaussian Process={\beta }_{0}+ {\sum _{J=1}^{J}{\beta }_{j}{x}_{j} \left(s\right)+ {f}_{Geology}+ {f}_{Land Use}+ {f}_{Terraces}+ {f}_{LSE}}_{}$$

2

where *β**0* is the global intercept, *β**j* are the fixed effects used to model continuous covariates and fGeology, fLand Use and fTerraces are the random effects for categorical properties, whereas fLSE is the random effect for the Latent Spatial Effect (LSE). We recall here that the terms fixed and random effects are synonyms of linear and nonlinear models in the context of Bayesian statistics. In other words, a single fixed effect to be exemplified here can be the Elevation, whose use as a linear covariate implies that the landslide intensity is assumed to increase or decrease as a function of the elevation with a fixed rate that cannot change at different levels of altitude. As for the random effects, these are more complex models to integrate how certain covariates contribute to the intensity estimates. Taking the Geology as an example, this is a discrete covariate whose classes contribute to the LGCP each one independently from the other. Another example of random effect can be instead visualized in the slope steepness. Here, we classified it into 20 classes in a pre-processing step, leading its original continuous information to become discrete as well. However, differently from the Geology case, each class retains some ordinal structure where all Grid-Cells contained in a given bin are always greater than all the Grid-Cells contained in the bin before and smaller than those in the subsequent bin. This is the reason why here we select a different approach for the slope steepness by introducing its effect onto the model with a random walk of the first order (see, Lombardo et al., 2018), for this structure allows for adjacent class dependence to be accounted for. Ultimately, the use of a Latent Spatial Effect is a very different tool from those explained above. In the landslide literature, Grid-Cells, catchments or any other mapping units are often modeled equally in space. In other words, mapping units that are close to each other are treated in the same way as those that are far apart. The only elements that allow the model estimates to change in space are the covariate values. However, in statistics there are solutions to inform the model about the spatial structure in the data. Specifically, adjacency matrix (see Fig. 3 of Opitz et al., 2022) can control this information which is then passed to the model as a latent covariate. The relation above corresponds to a Generalized Additive Mixed Models (GAMM, Steger et al., 2021), which we implement here in its Bayesian form via INLA (Bakka et al., 2018). We recall now two important properties of the landslide intensity. The intensity can always be converted into the most common susceptibility being the latter binary case a simpler realization of the count framework (Lombardo et al., 2019b). This can be achieved as follows:

$$Susceptibility=1- {e}^{{-\lambda }_{A}}$$

3

In addition, handling the intensity information over space is more convenient than doing the same in the susceptibility case. In fact, the susceptibility is mapping-unit dependent whereas the intensity benefits from the Poisson aggregation property across any spatial units. In this work, we use this property to aggregate *λ* values estimated for each GC contained in a given SU (see Fig. 5 in Lombardo et al., 2019a).

## 2.5. From landslide intensity to hazard and density

The landslide intensity has been shown to correlate with the total planimetric extent of landslides for each mapping unit (Lombardo et al., 2020). This contribution states that a model able to estimate landslide counts indirectly satisfies the current definition of hazard. Following the definition proposed by Corominas et al., 2014, hazard assessment aims to determine the spatial and temporal probability of slope failures occurrence, together with their mode of propagation, size and intensity. However, Lombardo and co-authors missed an important implication. In fact, if intensity and landslide sizes can be expressed one as the function of the other, this also means that one can convert landslide intensity maps into expected landslide size-related maps.

In this work, this possibility by estimating the intensity per SU and then estimating the landslide extent for each SU by multiplying the intensity for the mean landslide area was explored. We then also take a step further by dividing the estimated landslide areas for the corresponding SU size, thus returning the landslide density.

## 2.6. Performance assessment and model validation

Being our GAMM hierarchical in nature, we separately evaluate the performance at the level of the two mapping units. At the GC scale, where the data is almost binary in nature, Receiver Operating Characteristic (ROC) curves (Hosmer et al., 2013) was employed and their integral or AUC (Rahmati et al., 2019). At the SU level, the agreement between observed landslide counts and aggregated intensities via *χ**2* test and the Spearman correlation coefficient (*ρ*) were checked. In addition, the predictive performance of a fitted model was evaluated using the Leave-One-Out spatial Cross Validation (L1O-CV) method (Tanyas et al., 2019; Lombardo et al., 2020). The idea of the cross-validation is to perform several splitting of the data into training sample used for fitting the model, and into the validation sample (remaining data) employed for evaluating the predictive accuracy. In L1O procedure, each data object (in this case SU) is left out from the sample and used for validation. This means that different models are fitted (namely, 171), each one calibrates on 170 SU and alternatively predicted on the remaining one. It is important to note that this CV scheme was chosen to perturb the spatial dependence contained in the landslides distribution. In fact, if single GC were removed at random, the spatial scale at which these units act, would have not weakened the spatial structure captured via the latent spatial effect. Conversely, removing all the GC contained in each of the 171 SU partitioning the study area would have caused any residual spatial dependence to be weakened enough to test the model as a predictive tool rather than a simpler exploratory tool. A similar operation in the context of landslide modeling is extensively described in (Steger et al., 2016). In our case, we extracted all the GCs contained in each SU. Thus, being the SU different in size, a different number of GCs is extracted for each spatial cross-validation run.