Space-time landslide predictive modelling

Luigi Lombardo, Thomas Opitz, Francesca Ardizzone, Fausto Guzzetti,
Raphaël Huser

► To cite this version:

Luigi Lombardo, Thomas Opitz, Francesca Ardizzone, Fausto Guzzetti, Raphaël Huser. Space-time landslide predictive modelling. Earth-Science Reviews, 2020, 209, pp.103318. 10.1016/j.earscirev.2020.103318. hal-03164998

HAL Id: hal-03164998

https://hal.inrae.fr/hal-03164998v1

Submitted on 31 Jul 2024

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Space-time landslide predictive modelling

Luigi Lombardo<sup>a,*</sup>, Thomas Opitz<sup>b</sup>, Francesca Ardizzone<sup>c</sup>, Fausto Guzzetti<sup>c</sup>, Raphaël Huser<sup>d</sup>

<figure><img src="image_3.png" /></figure>

## ARTICLE INFO

### Keywords:

- Integrated Nested Laplace Approximation (INLA)

- Landslide hazard

- Landslide intensity

- Landslide susceptibility

- Log-Gaussian Cox Process (LGCP)

- Slope unit

- Space-time modelling

- Spatial point pattern

## ABSTRACT

Landslides are nearly ubiquitous phenomena and pose severe threats to people, properties, and the environment in many areas. Investigators have for long attempted to estimate landslide hazard in an effort to determine where, when (or how frequently), and how large (or how destructive) landslides are expected to be in an area. This information may prove useful to design landslide mitigation strategies, and to reduce landslide risk and societal and economic losses. In the geomorphology literature, most of the attempts at predicting the occurrence of populations of landslides by adopting statistical approaches are based on the empirical observation that landslides occur as a result of multiple, interacting, conditioning and triggering factors. Based on this observation, and under the assumption that at the spatial and temporal scales of our investigation individual landslides are discrete “point” events in the landscape, we propose a Bayesian modelling framework for the prediction of the spatio-temporal occurrence of landslides of the slide type caused by weather triggers. We build our modelling effort on a Log-Gaussian Cox Process (LGCP) by assuming that individual landslides in an area are the result of a point process described by an unknown intensity function. The modelling framework has two stochastic components: (i) a Poisson component, which models the observed (random) landslide count in each terrain subdivision for a given landslide “intensity”, i.e., the expected number of landslides per terrain subdivision (which may be transformed into a corresponding landslide “susceptibility”); and (ii) a Gaussian component, used to account for the spatial distribution of the local environmental conditions that influence landslide occurrence, and for the spatio-temporal distribution of “unobserved” latent environmental controls on landslide occurrence. We tested our prediction framework in the Collazzone area, Umbria, Central Italy, for which a detailed multi-temporal landslide inventory covering the period from before 1941 to 2014 is available together with lithological and bedding data. We subdivided the 79 km² area into 889 slope units (SUs). In each SU, we computed the mean and standard deviation of 16 morphometric covariates derived from a 10 m × 10 m digital elevation model. For 13 lithological and bedding attitude covariates obtained from a 1:10,000 scale geological map, we computed the proportion of each thematic class intersecting the given SU. We further counted how many of the 3,379 landslides in the multi-temporal inventory affect each SU and grouped them into six periods. We used this complex space-time information to prepare five models of increasing complexity. Our “baseline” model (Mod1) carries the spatial information only through the covariates mentioned above. It does not include any additional information about the spatial and temporal structure of the data, and it is therefore equivalent to the predominantly used landslide susceptibility model in the literature. The second model (Mod2) is analogous, but it allows for time-interval-specific regression constants. Our next two models are more complex. In particular, our third model (Mod3) also accounts for latent spatial dependencies among neighboring SUs. These are inferred for each of the six time intervals, to explain variations in the landslide intensity and susceptibility not explained by the thematic covariates. By contrast, our fourth model (Mod4) accounts for the latent temporal dependence, separately for each SU, disregarding neighboring influences. Ultimately, our most complex model (Mod5) contextually features all these relations. It contains the information carried by morphometric and thematic covariates, six time-interval-specific regression constants, and it also accounts for the latent temporal effects between consecutive slope instabilities at specific SUs as well as the latent spatial effects between adjacent SUs. We also show that the intensity is strongly related to the aggregated landslide area per SU. Because of this, our most complex model largely fulfills the definition of landslide hazard commonly accepted in the

<footnote><sup>a</sup> University of Twente, Faculty of Geo-Information Science and Earth Observation (ITC), PO Box 217, Enschede AE 7500, the Netherlands</footnote>

<footnote><sup>b</sup> INRA, Biostatistics and Spatial Processes, 84914 Avignon, France</footnote>

<footnote><sup>c</sup> Consiglio Nazionale delle Ricerche (CNR), Istituto di Ricerca per la Protezione Idrogeologica (IRPI), via Madonna Alta 126, 06128 Perugia, Italy</footnote>

<footnote><sup>d</sup> King Abdullah University of Science and Technology (KAUST), Computer, Electrical and Mathematical Sciences and Engineering (CEMSE) Division, Thuwal 23955-6900, Saudi Arabia</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: l.lombardo@utwente.nl (L. Lombardo).</footnote>
literature, at least for this study area. We quantified the spatial predictive performance of each of the five models using a 10-fold cross-validation procedure, and the temporal predictive performance using a leave-one-out cross-validation procedure. We found that Mod5 performed better than the others. We then used it to test a novel strategy to classify the model results in terms of both landslide intensity and susceptibility, which provides more information than traditional susceptibility zonations for land planning and management—hereafter we use the term “traditional” simply to refer to the majority of modelling procedures in the literature. We discuss the advantages and limitations of the new modelling framework, and its potential application in other areas, making specific and general hazard and geomorphological considerations. We also give a perspective on possible developments in landslide prediction modelling and zoning. We expect our novel approach to the spatio-temporal prediction of landslides to enhance the currently limited ability to evaluate landslide hazard and its temporal and spatial variations. We further expect it to lead to better projections of future landslides, and to improve our collective understanding of the evolution of complex landscapes dominated by mass-wasting processes under multiple geophysical and weather triggers.

# 1. Introduction

Landslides are ubiquitous in the hills, mountains, and high coasts that constellate the landmasses (Guzzetti et al., 2012), and in many areas they cause significant human, societal, economic, and environmental damage and costs (Brabb, 1989, 1991; Nadim et al., 2006; Dowling and Santi, 2014; Badoux et al., 2016; Grahn and Jaldell, 2017; Kirschbaum et al., 2009; Pereira et al., 2017; Froude and Petley, 2018; Salvati et al., 2018; Rossi et al., 2019). The reliable anticipation of landslides and their consequences is thus of primary importance.

Like for other natural hazards, the anticipation of landslides involves predicting “where” landslides can be expected (spatial prediction), “when” or how frequently they can be expected (temporal prediction), and “how many”, how large or destructive one should expect the landslides to be in an area (number, size, impact, destructiveness prediction) (Varnes, 1984; Guzzetti et al., 2005; Galli and Guzzetti, 2007; Tanyaş et al., 2018). The combined anticipation of “where”, “when” (or how frequently), and “how large” or destructive a landslide is expected to be, is called “landslide hazard” (Cruden and Fell, 1997; Hungr et al., 1999; Guzzetti, 2005; Guzzetti et al., 2005; Reichenbach et al., 2005; Fell et al., 2008; Lari et al., 2014). Differently from other natural hazards, two distinct types of predictions are possible for landslides, namely, (i) the prediction of single landslides, i.e., the anticipation of the behaviour of a single slope, or a portion of a slope, and (ii) the prediction of populations of landslides, i.e., the anticipation of the behaviour of many (tens to several tens of thousands) landslides occurring in an area, and their spatial and temporal evolution. In this work, we focus on the prediction of populations of landslides in an area, and we do not consider the anticipation of the behaviour of single slopes or individual landslides. For this purpose, we exploit a unique multi-temporal inventory of landslides that occurred over a multi-decade period in an area of Central Italy, which we use to fit and validate a set of five Bayesian statistical models constructed under the general assumption that landslides are a stochastic point process (Cox, 1965; Cox and Isham, 1980; Chiu et al., 2013).

The paper is organised as follows. We begin, in §2, by providing background information on traditional spatial and temporal landslide predictive modelling approaches, and their limitations. This is followed, in §3, by a description of the study area of Collazzone, Italy, and, in §4, of the landslide, the morphological, and the thematic data used, of our choice of the modelling mapping unit, and the pre-processing steps. Next, in §5, we describe five different spatial statistical models that we have implemented, consisting of: (i) a baseline model where the landslide spatial dimension is only carried through the explanatory variables; (ii) an improved version of the baseline model which allows for time-interval-specific regression constants; and three extensions to the second baseline model which account for (iii) spatial, (iv) temporal, and (v) spatio-temporal random effects acting at a latent level. By “latent level”, we refer to a class of properties or effects, which are present in the landslide distribution but they are not directly observable or

expressed as explanatory variables; they are thus “latent” in our statistical models (e.g., Gebregziabher and DeSantis, 2010; Lombardo et al., 2020). This is followed, in §6, by the presentation and comparison of the results for the five spatial statistical models, and the associated calibration and validation diagnostics. In §7, we discuss the results obtained, and we provide geomorphological insight on the performed statistical inference. Lastly, in §8, we summarise the lessons learnt and we outline the remaining challenges.

# 2. Prediction of landslide occurrence

In the geomorphological literature, most of the attempts at predicting the occurrence of populations of landslides in an area are based on the empirical observation that landslides are spatially and temporally discrete events that occur as a result of multiple, interacting, conditioning and triggering factors. The conditioning factors primarily influence where landslides can occur, whereas the triggering factors drive the landslides' onset, i.e., the time or period of occurrence of the slope failures. Together, the conditioning and the triggering factors control the extent of the area affected by landslides and the size distribution of the slope failures, which is linked to the landslide impact and destructiveness. Given the complexity and the variability of landslide processes, which depend on, e.g., the soil, rock, and other landscape characteristics, and on the weather or seismic triggers, and because the exact or approximate locations of the landslides are unknown before they occur, individual slope failures in a population of landslides can be considered a realisation of a stochastic process (Das et al., 2012; Lombardo et al., 2014), and are modelled accordingly.

Many approaches have been proposed to assess the landslide “susceptibility”, which refers in the geomorphological literature to the spatially-varying, time-independent likelihood of landslides occurring in an area given the local terrain conditions (Brabb, 1985; Chung and Fabbri, 1999; Guzzetti et al., 1999, 2005; Reichenbach et al., 2018). These approaches can be loosely grouped into five main categories (Guzzetti, 2005), i.e., (i) direct geomorphological mapping (Verstappen, 1983; Hansen et al., 1995; Reichenbach et al., 2005), (ii) analysis of landslide inventories (Campbell, 1973; DeGraff and Canuti, 1988; Moreiras, 2004), (iii) heuristic, index-based methods (Nilsen and Brabb, 1977; Posner and Georgakakos, 2015), (iv) deterministic, physically-based, conceptual models (Ward et al., 1981, 1982; Montgomery and Dietrich, 1994; Dietrich et al., 2001; Goetz et al., 2011; Bout et al., 2018), and (v) statistical prediction models (Carrara, 1983; Chung and Fabbri, 1999; Guzzetti et al., 1999; Van Westen et al., 2006; Lombardo et al., 2016b; Reichenbach et al., 2018). Each of these approaches has potential advantages and inherent limitations (Guzzetti, 2005; Van Westen et al., 2006; Lombardo et al., 2015).

Geomorphological mapping depends entirely on the skills and experience of the investigators. It may provide reliable results, but it is difficult to reproduce, impractical over large areas, and inadequate for quantitative hazard assessments (Guzzetti et al., 1999; Van Westen
et al., 1999). Analysis of the inventories depends on the quality and completeness of the available landslide maps (Tanyaş and Lombardo, 2020). Where an inventory is incomplete, or wrong, the susceptibility assessment will be underestimated, or biased (Guzzetti et al., 2012). Direct, heuristic/geomorphological mapping methods rely on decisions made by the investigator on how to weigh certain properties on the assumption that all landslide causes in an area are known. Moreover, they produce qualitative and subjective predictions unsuited for quantitative hazard assessments (Soeters and Van Westen, 1996; Guzzetti et al., 1999; Leoni et al., 2009). Physically- or process- based models exploit the existing understanding of the mechanical laws that control slope instability (Guzzetti et al., 1999). Their limitation lies in the inherent simplicity of the modelling equations that may not capture the complex interactions controlling the slope stability/instability conditions. Furthermore, the physically-based models require large datasets to describe the surface and subsurface mechanical and hydrological properties of the terrain, which are difficult and expensive to obtain. As a result, physically-based models are used chiefly for small or very small areas (e.g., Montgomery and Dietrich, 1994; Chakraborty and Goswami, 2016; Seyed-Kolbadi et al., 2019), albeit a few examples also exist of applications for large areas (e.g., Gorsevski et al., 2006; Raia et al., 2014; Alvioli and Baum, 2016).

Lastly, statistical approaches aim at exploiting the “functional” relationships existing between a set of instability factors, and the past and present distribution of landslides obtained typically from a landslide inventory map (Carrara, 1983; Duman et al., 2005; Guzzetti et al., 2012), or a landslide catalogue (Van Den Eeckhaut and Hervás, 2012). The large number of statistically-based approaches proposed in the literature (Reichenbach et al., 2018) almost invariably exploit classification methods, and provide probabilistic estimates suited for quantitative hazard assessments. Statistical models can be constructed using a variety of thematic and environmental variables obtained from existing maps or by processing remotely sensed images and data, in different landscape and environmental settings, covering a broad range of scales and study-area sizes. The dependent variable is obtained from different types of landslide inventory maps (Guzzetti et al., 2012) or landslide catalogues (Van Den Eeckhaut and Hervás, 2012), and is typically used in a binary structure, expressing the presence or absence of landslides in each mapping unit, where a terrain mapping unit is a regular or irregular geographical subdivision (e.g., a pixel, unique condition, slope or hydrological unit, administrative subdivision (Guzzetti et al., 1999; Guzzetti, 2005; Van Westen et al., 2006)) used to partition a study area. The fitted model is then used to assess the landslide susceptibility for each mapping unit (Guzzetti, 2005). In this framework, the uncertainty in the landslide mapping procedure (Santangelo et al., 2015a) may inevitably bias the susceptibility models and therefore, whenever incomplete landslide inventories are used, an analytical step aimed at assessing the effect of positional errors should always be included in the modelling routine (Steger et al., 2016a, 2016b).

Focusing on statistically-based susceptibility approaches, a limitation of the traditional and of most of the current models is that they predict only whether a mapping unit is expected to have (or not have) landslides, regardless of the number or size of the expected failures in each unit. In a population of landslides, the size (i.e., length, width, area, volume) of the slope failures is linked to the number of the failures. Hovius et al. (1997) and Malamud et al. (2004), among others, have shown that the area of a landslide obeys empirical probability distributions that control the average size and relative proportion of the landslides in an area. Similar empirical dependencies have been found for landslide volume (Brunetti et al., 2009), for the landslide area to volume ratio (Guzzetti et al., 2009; Larsen et al., 2010), and more recently for the landslide width to length ratio (Taylor et al., 2018).

In an attempt to overcome the inherent inability of landslide susceptibility models to predict the number of the expected landslides, Lombardo et al. (2018a) have proposed a novel framework for landslide intensity assessment. The approach treats single landslides in a

population of landslides as individual realisations from a continuous—space point process (Cox and Isham, 1980; Chiu et al., 2013). This is different from the discrete-space binary presence-absence model adopted by the traditional statistically-based susceptibility models (Carrara, 1983; Guzzetti et al., 1999; Guzzetti, 2005; Lombardo and Mai, 2018; Reichenbach et al., 2018). As a result, the approach aims at predicting the number of the expected landslides in each mapping unit adopting a Poisson distribution, whose mean is linked to the unknown landslide intensity that can be estimated from a landslide dataset. The approach was applied successfully to model populations of rainfall-induced (Lombardo et al., 2018a, 2019b) and seismically-triggered (Lombardo et al., 2019a) landslides in different morphological, geological, and climatic settings.

A second limitation of most statistically-based models is that they typically do not consider the spatial relationships between landslide occurrences in different terrain mapping units. Landslide occurrences are often assumed to be independent given the terrain conditions. In other words, the geographical location of the mapping units is not explicitly taken into account, so that adjacent, neighbouring, and distant units are considered equally by the models. Not considering the spatial dependencies (or lack thereof) among units placed in different locations in a study area may result in poorer susceptibility models and associated zonations (Reichenbach et al., 2018; Lombardo et al., 2018a). And it may bias the assessment of predictive performance if validation techniques do not take into account spatial dependence (see, e.g., Brenning, 2005; Goetz et al., 2015).

A third limitation of most statistically-based models is the fact that they consider the likelihood of landslides occurring in an area to be constant in (i.e., independent of) time (e.g., Meusburger and Alewell 2009; Cama et al., 2015). In the long run (i.e., hundreds to thousands of years), the assumption of stationarity does not hold because the landslide triggering conditions (e.g., the frequency of high or prolonged rainfall periods, of snow melt events, or of earthquakes), as well as the predisposing conditions (e.g., land use or land cover) may change with time, inevitably changing the landslide susceptibility. Recent empirical evidence shows that in some landscapes, even for short periods (i.e., tens of years), when a landslide occurs it may become an “attractor” for future landslides, with new landslides being more likely to occur inside or in the immediate vicinity of the previous landslides. Samia et al. (2017a, 2017b) called this effect of short-term temporal dependence “landslide path dependency”, and found that the spatio-temporal dependency of new landslides on old landslides disappears in their study area, in Umbria, Central Italy—the same study area of this work—approximately after 10–15 years. The evidence violates the assumption that susceptibility is constant in time, even for short periods.

In the literature, approaches to predict the temporal occurrence of landslides are equally if not more diversified. Depending on the scope, the temporal coverage, the return time of the predictions, and the extent of the study area, methods include (i) the use of empirical rainfall thresholds for the possible occurrence of landslides (Glade et al., 2000; Dai and Lee, 2001; Crosta and Frattini, 2000; Aleotti, 2004; Guzzetti et al., 2007, 2008; Saito et al., 2010; Ko and Lo, 2018; Segoni et al., 2018; Guzzetti et al., 2020), (ii) physically-based, coupled, distributed rainfall and infiltration slope stability models (Montgomery and Dietrich, 1994; Van Asch et al., 1996; Baum et al., 2008; Lanni et al., 2013; Formetta et al., 2014; Reid et al., 2015; Formetta et al., 2016; Alvioli and Baum, 2016; Bout et al., 2018), and (iii) the analysis of time series of historical landslides and landslide events (Crovelli and Coe, 2009; Rossi et al., 2010b; Witt et al., 2010). Only physically-based models (inherently) consider the spatio-temporal interactions that condition landslide occurrence, which are not considered by the threshold models or by the analysis of the historical records. However, the physically-based models are generally applicable only to small areas and for short periods of time (a few hours to a few days). Furthermore, they are not suited for the spatio-temporal modelling of landslide occurrence over large areas and for long periods, which is essential for
<figure><img src="image_4.png" /></figure>

Inventory type: H, historical. E, event. S, seasonal

(caption on next page)
land planning and management.

In this work, we construct innovative spatial statistical models that consider (i) spatial, (ii) temporal, and (iii) spatio-temporal landslide latent effects among: adjacent terrain mapping units, same mapping units but subsequent time intervals, and both conditions together, respectively. We consider this an improvement over the existing approaches to predict the spatio-temporal occurrence of populations of landslides in an area.

### 3. Study area

The area of Collazzone, Umbria, Central Italy, is well studied and represents a unique site where landslides have been mapped repeatedly over a large timespan. A description of the area can be found in Guzzetti et al. (2006a, 2006b), Ardizzone et al. (2007), Galli et al. (2008) and other references therein. Overall, our study area extends over ≈79 km², with terrain elevation between 145 m along the Tiber River flood plain NNW of Todi, and 634 m at Monte di Grutti (Fig. 1). The landscape is predominantly hilly, and lithology and the attitude of bedding planes control the morphology of the slopes, and the presence and abundance of the landslides. In the area sedimentary rocks crop out, Cretaceous to Holocene in age, including recent fluvial deposits, continental gravel, sand and clay, travertine, layered sandstone and marl, and thinly layered limestone. The climate is Mediterranean, with precipitation falling mostly in the period from October to December, and from February to May. Intense rainfall events or prolonged rainfall periods are the primary natural triggers of landslides in the area (Ardizzone et al., 2013), followed by the rapid melting of snow (Cardinali et al., 2000). Landslides are abundant and cover 17.05 km²—corresponding to a density of approximately 43 landslides per square kilometre—and range in age, type, morphology, and volume from very old, partly eroded, large and deep-seated slides and slide-earth flows, to young and shallow slides and flows (Hungr et al., 2014). Landslides are most abundant in the cultivated areas and are rare in the forested terrain, indicating a dependence of the landslides on the agricultural practices, at least in the last 20 years (Galli et al., 2008; Mergili et al., 2014).

### 4. Data compilation and pre-processing

#### 4.1. Landslide data

We obtained the landslide information from a pre-existing multi-temporal landslide inventory prepared at 1:10,000 scale through the systematic visual interpretation of five sets of aerial photographs taken in 1941, 1954, 1977, 1985, and 1997, and of five stereoscopic, panchromatic and multi-spectral satellite image pairs acquired in August 2009, March 2010, May 2010, April 2013, and April 2014 (Table 1), supplemented by field checks and surveys executed in various periods from 1998 to 2004, in May 2004, and in December 2005 (Guzzetti et al., 2006a; Ardizzone et al., 2007; Galli et al., 2008; Ardizzone et al., 2013). The multi-temporal inventory shows the location, surface area, type (Hungr et al., 2014), and estimated age (Cruden and Varnes, 1996) of 3,379 landslides, ranging in size from $A_L = 5.8 \times 10^3$ m² to $A_L = 1.5 \times 10^6$ m² (mean, $\mu = 6.9 \times 10^3$ m², standard deviation, $\sigma = 3.2 \times 10^4$ m²), for a total area covered by landslides of $A_{LT} = 17.05$ km². In the inventory, landslide age was defined as very old, old (predating 1941), or recent (from 1941 to 2014), using photo-interpretation criteria and field evidence, despite some ambiguity in the definition of the age of a landslide based on its morphological appearance (McCalpin, 1984; Guzzetti et al., 2006a).

Key to our study is the fact that in the multi-temporal inventory

landslides are separated into nineteen, irregularly spaced “temporal slices”, each shown by a different colour in Fig. 1. Landslides in the temporal slices before 2000 were detected and mapped using black and white (panchromatic) and colour (1977) aerial photographs, whereas landslides between 2009 and 2014 were detected and mapped using VHR stereoscopic satellite images (Table 1), and directly in the field. Visual interpretation of the stereoscopic satellite images was performed using different software, including: (i) ERDAS IMAGINE® and Leica Photogrammetry Suite (LPS) for block orientation of the stereo images; (ii) Stereo Analyst for ArcGIS® for image visualisation and landslide mapping; and (iii) StereoMirror™ hardware technology to obtain 3D views of the VHR satellite images.

The same interpretation criteria were adopted to identify and map the landslides on aerial photographs and the satellite images (Guzzetti et al., 2012; Ardizzone et al., 2013; Murillo-Garcia et al., 2015). In each set of aerial photographs and in each pair of satellite images, landslides that appeared “fresh” were given the date (i.e., year) of the aerial photographs, or the date (i.e., month and year) of the satellite images. The “non-fresh” landslides were attributed to the period (i.e., the “temporal slice”) between two successive sets of aerial photographs or satellite image pairs. Landslides mapped in the field after single or multiple rainfall events between 1998 and 2013 were given the date (i.e., month and year) of the field surveys. The different methods and instruments used to interpret the aerial photographs and the satellite images, and the fact that some of the landslides were mapped in the field, have conditioned the completeness and accuracy of the landslide information in the multi-temporal inventory, which is therefore not constant throughout the time slices. In general, more recent time slices showing event or seasonal inventories (E and S, respectively, in Fig. 1) have more numerous landslides of small sizes, whereas historical inventories (H, in Fig. 1) show larger landslides, and lack landslides of very small size. This is a known bias of the multi-temporal inventory used in our study (Guzzetti et al., 2006a; Ardizzone et al., 2007; Galli et al., 2008; Ardizzone et al., 2013).

#### 4.2. Mapping unit

Prediction of landslide occurrence in an area requires the preliminary selection of a suitable terrain mapping unit, i.e., a subdivision of the terrain that maximises the within-unit (internal) homogeneity and the between-unit (external) heterogeneity across distinct physical or geographical boundaries (Carrara et al., 1991; Guzzetti, 2005). In agreement with previous studies in the same (Guzzetti et al., 2006b), in similar (Carrara et al., 2003), and in different (Van Den Eeckhaut et al., 2009a; Erener and Düzgün, 2012; Castro Camilo et al., 2017; Amato et al., 2019) areas, we select the “slope unit” (SU) as the mapping unit of reference. By definition, a SU is a terrain geomorphological unit bounded by drainage and divide lines, and corresponds to a slope, a combination of adjacent slopes, or a small catchment (Carrara et al., 1991; Alvioli et al., 2016). We use a subdivision of the study area (i.e., our spatial domain) into 889 SUs ranging in size from $A_L = 6.17 \times 10^2$ m² to $A_L = 1.4 \times 10^6$ m² (mean, $\mu = 8.9 \times 10^4$ m², standard deviation, $\sigma = 1.1 \times 10^5$ m²), corresponding to an average density of one SU approximately every 0.1 km². Panel A of Fig. 2 portrays the geographical distribution of the SUs used in the study. Notably, this slope unit partition is the same one as that adopted in previous studies in the same area (see, e.g., Guzzetti et al., 2006b).

#### 4.3. Thematic data and explanatory variables

To support the modelling procedure, we used a set of morphometric

<figure><img src="image_5.png" /><figcaption>Fig. 1. Collazzone study area, Umbria, Central Italy. The map shows the multi-temporal landslide inventory, morphology, and hydrology of the study area. Coloured polygons are landslides of various ages or periods mapped through the visual interpretation of aerial photographs or satellite images of different vintages, or through field work. Individual inventories shown by different colours are of three types: E, event; S, seasonal; H, historical. See text for further explanation.</figcaption></figure>
<caption>**Table 1**</caption>

Main characteristics of the aerial photographs and the optical satellite images used to prepare the multi-temporal landslide inventory for the Collazzone study area, Umbria, Central Italy, used in this work and shown in Fig. 1. Legend: GSD, ground sampling distance.

<table><thead><tr><th>Year</th><th>Period</th><th>Type</th><th>Nominal scale</th><th>Mode</th><th>Source</th></tr></thead><tbody><tr><td>1941</td><td>Summer</td><td>Panchromatic</td><td>1:18,000</td><td>Stereo</td><td>Aerial photographs</td></tr><tr><td>1954</td><td>Spring-Summer</td><td>Panchromatic</td><td>1:33,000</td><td>Stereo</td><td>Aerial photographs</td></tr><tr><td>1977</td><td>Spring-Summer</td><td>Colour</td><td>1:13,000</td><td>Stereo</td><td>Aerial photographs</td></tr><tr><td>1985</td><td>July</td><td>Panchromatic</td><td>1:15,000</td><td>Stereo</td><td>Aerial photographs</td></tr><tr><td>1997</td><td>April</td><td>Panchromatic</td><td>1:20,000</td><td>Stereo</td><td>Aerial photographs</td></tr><tr><td>2009</td><td>12 August</td><td>Panchromatic</td><td>GSD = 0.41 m</td><td>Stereo</td><td>GeoEye-1</td></tr><tr><td>2010</td><td>March</td><td>Panchromatic</td><td>GSD = 0.50 m</td><td>Stereo</td><td>WorldView-1</td></tr><tr><td>2010</td><td>27 May</td><td>Panchromatic</td><td>GSD = 0.41 m</td><td>Stereo</td><td>GeoEye-1</td></tr><tr><td>2013</td><td>13 April</td><td>Panchromatic</td><td>GSD = 0.50 m</td><td>Stereo</td><td>GeoEye-1</td></tr><tr><td>2014</td><td>14 April</td><td>Multispectral</td><td>GSD = 2.00 m</td><td>Stereo</td><td>WorldView-2</td></tr></tbody></table>

and thematic variables, which we list in Table 2. Throughout the paper, the term covariate is used invariably for both quantitative and categorical explanatory variables. The 16 morphometric covariates were derived from a 10 m × 10 m resolution digital elevation model (DEM) obtained through the automatic interpolation of 10 m and 5 m interval contour lines taken from 1:10,000 scale topographic base maps made accessible by the Umbria Region government (link here). To represent morphometric covariates at the SU scale, we computed the two main statistical properties (mean μ, and standard deviation σ) of each covariate for each respective SU. We further obtained the thematic covariates (Table 2) from a geological map prepared at 1:10,000 scale by Guzzetti et al. (2006a), and used in previous landslide susceptibility and hazard studies in the same area (Guzzetti et al., 2006a; Ardizzone et al., 2007; Galli et al., 2008). From the geological map we extracted information on nine lithological units, and four bedding domain classes. We selected these covariates because bedrock geology and the attitude of bedding planes have been shown to control the presence (or absence) and the abundance of landslides in our study area (Guzzetti et al., 2006a; Marchesini et al., 2015). To summarise the geologic and bedding information for each SU, we computed the ratio between the relative extent of each categorical class in a given SU and the total extent of the SU. As a result, we transformed the original categorical

information into a continuous one, expressing the proportion of area coverage per class in each SU.

### 4.4. Pre-processing

We initially tested our modelling framework (see Section 5) using the 19 separate time slices shown in Fig. 1, where the original 19 landslide count distributions represented our target variable. However, these preliminary tests did not produce satisfying predictive performances (unreported results). Possible reasons for the unsatisfactory results included (i) the irregular time-span of the landslide inventories, ranging from 1 to 23 years, (ii) the highly variable number of landslides in each time slice, from 303 landslides for the 1985–1996 time slice, to 866 landslides for the 1941–1954 time slice, and (iii) the fact that some time slices have only a few landslides (Fig. 3A, dark grey). To address the issue, we aggregated the landslides both in space and time.

To aggregate the data over space, we looked into the inventories and realised that some of the landslides have a large areal extent at certain times, leading the instability to affect several SUs simultaneously. In such cases, treating these landslides as a single points (i.e., assigning one count to a single SU) would neglect the areal extent of the landslide phenomenon. Therefore, in cases where a landslide covers

<figure><img src="image_6.png" /><figcaption>**Fig. 2.** Collazzone study area, Umbria, Central Italy. (A) Geographical subdivision of the study area into 889 slope units (SUs). (B) Map showing number of adjacent SUs, and adjacency structure. (C) Adjacency graph showing links (blue lines) between adjacent SU centroids (grey dots), and adjacency matrix. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
<caption>**Table 2**</caption>

Morphometric (M), lithological (L), and bedding-related (structural, S) explanatory variables (covariates) used in the study for space-time landslide predictive modelling in the Collazzone area, Umbria, Central Italy (see Fig. 1).

<table><thead><tr><th>Type</th><th>Variable</th><th>Description</th><th>Source</th><th>Reference</th></tr></thead><tbody><tr><td>M</td><td>ELEVμ</td><td>Terrain elevation mean</td><td>10 m × 10 m DEM</td><td>1</td></tr><tr><td>M</td><td>ELEVσ</td><td>Terrain elevation st.dev.</td><td>10 m × 10 m DEM</td><td>1</td></tr><tr><td>M</td><td>SLOPEμ</td><td>Terrain slope mean</td><td>10 m × 10 m DEM</td><td>2</td></tr><tr><td>M</td><td>SLOPEσ</td><td>Terrain slope st.dev.</td><td>10 m × 10 m DEM</td><td>2</td></tr><tr><td>M</td><td>ENESμ</td><td>Eastness mean</td><td>10 m × 10 m DEM</td><td>3</td></tr><tr><td>M</td><td>ENESσ</td><td>Eastness st.dev.</td><td>10 m × 10 m DEM</td><td>3</td></tr><tr><td>M</td><td>NNESμ</td><td>Northness mean</td><td>10 m × 10 m DEM</td><td>3</td></tr><tr><td>M</td><td>NNESσ</td><td>Northness st.dev.</td><td>10 m × 10 m DEM</td><td>3</td></tr><tr><td>M</td><td>PLCRμ</td><td>Planar curvature mean</td><td>10 m × 10 m DEM</td><td>4</td></tr><tr><td>M</td><td>PLCRσ</td><td>Planar curvature st.dev.</td><td>10 m × 10 m DEM</td><td>4</td></tr><tr><td>M</td><td>PRCRμ</td><td>Profile curvature mean</td><td>10 m × 10 m DEM</td><td>4</td></tr><tr><td>M</td><td>PRCRσ</td><td>Profile curvature st.dev.</td><td>10 m × 10 m DEM</td><td>4</td></tr><tr><td>M</td><td>RSPμ</td><td>Relative slope position mean</td><td>10 m × 10 m DEM</td><td>5</td></tr><tr><td>M</td><td>RSPσ</td><td>Relative slope position st.dev.</td><td>10 m × 10 m DEM</td><td>5</td></tr><tr><td>M</td><td>TWIμ</td><td>Topographic wetness index mean</td><td>10 m × 10 m DEM</td><td>6</td></tr><tr><td>M</td><td>TWIσ</td><td>Topographic wetness index st.dev.</td><td>10 m × 10 m DEM</td><td>6</td></tr><tr><td>L</td><td>AD_R</td><td>Alluvial sediment</td><td>Lithological map, 1:10,000</td><td>7</td></tr><tr><td>L</td><td>C_R</td><td>Clay</td><td>Lithological map, 1:10,000</td><td>7</td></tr><tr><td>L</td><td>G_R</td><td>Gravel</td><td>Lithological map, 1:10,000</td><td>7</td></tr><tr><td>L</td><td>L_R</td><td>Limestone</td><td>Lithological map, 1:10,000</td><td>7</td></tr><tr><td>L</td><td>M_R</td><td>Marl</td><td>Lithological map, 1:10,000</td><td>7</td></tr><tr><td>L</td><td>S_R</td><td>Sand</td><td>Lithological map, 1:10,000</td><td>7</td></tr><tr><td>L</td><td>SGC_R</td><td>Gravel, sand, clay</td><td>Lithological map, 1:10,000</td><td>7</td></tr><tr><td>L</td><td>SS_R</td><td>Sandstone</td><td>Lithological map, 1:10,000</td><td>7</td></tr><tr><td>L</td><td>T_R</td><td>Travertine</td><td>Lithological map, 1:10,000</td><td>7</td></tr><tr><td>S</td><td>AS_R</td><td>Anaclinal slope</td><td>Bedding map, 1:10,000</td><td>7</td></tr><tr><td>S</td><td>CS_R</td><td>Cataclinal slope</td><td>Bedding map, 1:10,000</td><td>7</td></tr><tr><td>S</td><td>OS_R</td><td>Orthogonal slope</td><td>Bedding map, 1:10,000</td><td>7</td></tr><tr><td>S</td><td>US_R</td><td>Unbedded sediment</td><td>Bedding map, 1:10,000</td><td>7</td></tr></tbody></table>

References: 1, http://www.umbriageo.regione.umbria.it/pagina/distribuzione-carta-tecnica-regionale-vettoriale-1; 2, Zevenbergen and Thorne (1987); 3, Lombardo et al. (2018b); 4, Heerdegen and Beran (1982); 5, Böhner and Selige (2006); 6, Beven and Kirkby (1979); 7, Guzzetti et al. (2006a).

more than one SU, we assigned a count to multiple SUs under the constraint that the intersections between the landslide and the SUs are larger than 2% of the area of the SU (Carrara and Guzzetti, 2013). In turn, this procedure generated a moderately larger landslide count per SU than the original landslide count, as shown in Fig. 3A.

The second aggregation scheme involved the temporal dimension of

the spatially-aggregated landslide counts. Specifically, we further aggregated the 19 temporal slices in a nearly regular temporal grid with intervals of approximately 15 years, with the exception of the stand-alone snow-melt event. The procedure produced six “time intervals”, each composed of one to eight time slices, which we will refer to as T1 to T6 in the rest of the manuscript. Panel B of Fig. 3 illustrates the

<figure><img src="image_7.png" /><figcaption>**Fig. 3.** Collazzone study area, Umbria, Central Italy. Spatial and temporal aggregation schemes of landslide counts. (A): Temporal distribution of the original counts obtained by considering landslides as pure points (dark grey); and, spatial aggregation by considering the areal extent of landslides (light grey). Hence, where landslides covered more than one SU, we repeated the count over multiple SUs to respect the geomorphological and areal realisation of the instability process. (B): Spatially aggregated landslide count per SU displayed per “temporal slices” for the 19 inventories shown in Fig. 1; and, associated temporal aggregation scheme to merge the 19 time slices (single bars) into six “time intervals” (groups of bars shown by the same colour). (C): Summary stacked bar chart of the aggregated landslide counts distribution per SU for six near-regular time intervals, from T1 to T6: the x-axis reports the SUs' rank, from SUs having zero landslides (left) to SUs having up to 25 landslides (right). The y-axis shows the proportion of landslide counts across the six different time intervals, shown by different colours. White characters correspond to the counts per time interval. More information on SUs is provided in Section 4.2.</figcaption></figure>
distribution of landslide counts per SU for each time interval, and Panel C of Fig. 3 summarises our final aggregation scheme.

In addition to the covariate preparation and the spatio-temporal aggregation, the pre-processing phase involved creating the so-called adjacency matrix (Zhang et al., 2010). The ($i_1, i_2$)-th entry of the adjacency matrix is equal to one if the $i_1$-th SU shares a border with the $i_2$-th SU, and it is zero otherwise. Therefore, it is a symmetric matrix (i.e., if SU1 is adjacent to SU2, then SU2 is also adjacent to SU1), which indicates the neighbourhood structure between SUs, and it provides fundamental information to build the space and space-time models presented in this work. Maps (B) and (C) in Fig. 2 summarise the main steps involved in the calculation of the adjacency matrix, and illustrate the adjacency graph structure.

# 5. Modelling framework

## 5.1. Fundamentals of point processes

Conceptually, we identify a landslide with a single point ($s_i, t_i$), where $s_i$ is the spatial location and $t_i$ is time. In practice, we may also consider the spatial dimension of large landslides by counting them more than once. The time resolution here corresponds to the six time intervals, such that $t_i \in \{1, 2, ..., 6\}$. In our modelling approach detailed below, we do not need the exact position $s_i$ of the landslides for the purpose of estimation and mapping at the resolution of the SUs, but only the index of the corresponding SUs and the corresponding SU adjacency graph structure (Fig. 2). Our approach based on SUs is therefore robust against positional errors of landslide data, which is a major advantage to prevent potential modelling and prediction biases.

Our fundamental modelling assumption is that the space-time points ($s_i, t_i$) identifying landslide occurrences stem from an unobserved intensity function $\lambda(s, t)$ that varies over space and time. The interpretation of this intensity function is that we observe (approximately) $\lambda(s, t)$ events on average per spatio-temporal unit around $(s, t)$. For an arbitrary domain $D$ stretching over space and time, the mean number of observed events is given by the integral $\int_D \lambda(s, t) dsdt$. The natural distribution of such event counts is the Poisson distribution with intensity given by this integral. We assume that this stochastic mechanism generates the observed spatio-temporal point pattern, and we call it a Poisson point process. In a Poisson point process, the individual landslides arise independently from each other. Given two locations ($s_1, t_1$) and ($s_2, t_2$) with $\lambda(s_2, t_2) > 0$, the ratio $\lambda(s_1, t_1)/\lambda(s_2, t_2)$ indicates how much more likely it is to observe a landslide at ($s_1, t_1$) in comparison to ($s_2, t_2$).

In practice, we construct a model for the intensity function $\lambda(s, t)$ which incorporates covariate effects, and which may further capture structured variations of the space-time intensity surface of random nature, i.e., not explained by the available covariate information. Such random effects can be considered as “unobserved” effects, hence covariates during the modelling process, whose signals influence the landslide space-time pattern in the data. For instance, in case of event-based landslide studies, the precipitation or seismic triggers may not be included among the covariate set (e.g., if they are unobserved), although they influence the concentration of landslides at specific locations. Nevertheless, complex spatial models can retrieve the pattern of the unobserved trigger from the landslide distribution itself, as demonstrated by Lombardo et al. (2018a, 2019b) for storm-induced and by Lombardo et al. (2019a) for earthquake-induced landslides. As for temporal unobserved effects, complex temporal models can retrieve the influence of the “landslide path dependency” recognised by Samia et al. (2018).

When allowing for such random components in the intensity function $\lambda$, the resulting stochastic model for the observed point pattern is called a Cox point process (Cox, 1965), and we use the uppercase notation $\Lambda(s, t)$ for the intensity function to highlight its stochastic behaviour. More specifically, if we opt for the flexible and convenient choice

of additive random effects in the log-intensity $\log(\Lambda(s, t))$ possessing Gaussian process distributions, we obtain a Log-Gaussian Cox Process (Møller et al., 1998; Basu and Dassios, 2002; Diggle et al., 2013), LGCP in short.

For our dataset, we can rewrite the general structure of such models by taking into account the spatial (889 SUs) and temporal (6 intervals) resolution and by using the surface area $|A_i|$, $i = 1, ..., 889$, of SUs. We make the natural assumption that the average number of landslides observed within a SU $s_i$ is proportional to its surface area $|A_i|$, given that all other predictor components are the same. Furthermore, we model a separate spatial intensity for each of the temporal intervals $j = 1, ..., 6$, and so we write $\Lambda_j(s_i) = |A_i| \times \tilde{\Lambda}_j(s_i)$ for the integrated intensity of the $i$-th SU in the $j$-th temporal interval, where $\tilde{\Lambda}_j(s_i)$ can be interpreted as the intensity of a (rescaled) unit-area SU with the same characteristics.

We further write $N_{ij} \in \{0, 1, 2, ...\}$ for the number of landslides observed in SU $i$ during temporal interval $j$. Then, given the intensity values $\Lambda_j(s)$, the model has the following structure:

$$N_{ij} | \Lambda_j(s) \sim \text{Poisson}(|A_i| \times \tilde{\Lambda}_j(s_i)), \quad i = 1, ..., 889, \quad j = 1, ..., 6. \quad (1)$$

We will formulate five regression models (technically speaking, Poisson regressions with a log-link function) in §5.4–§5.7, with two variants of a baseline model without random effects presented in §5.4. These models integrate observed covariate effects and random effects into the log-intensity $\log(\Lambda_j(s))$ in an additive way. The models follow the general structure

$$\log(\Lambda_j(s)) = \log(|A_i|) + \text{fixed covariate effects} + \text{random effects} \quad (2)$$

where the random effect component are not accounted for in our baseline models. The term $\log(|A_i|)$ represents a fixed deterministic offset. We here adopt a discrete perspective on time where the study period is represented through 6 time intervals, not necessarily of the same length. Then, the intensity function $\Lambda_j$ is expressed for a time unit given by the length of the $j$-th interval, $j = 1, ..., 6$. We do not have to explicitly specify the length of each of the 6 intervals, which is of advantage here since there is no unique sensible partition of the study period into 6 intervals. Moreover, the first interval with no specific lower endpoint could be viewed as stretching towards minus infinity, and defining a time unit common to all 6 intervals may be awkward.

## 5.2. The Bayesian modelling paradigm

Before presenting the specific structure of the five models we tested, we first recall the general idea of fully Bayesian modelling as implemented here. We aim to simultaneously estimate the latent intensity function $\Lambda_j(s)$, and more specifically its components such as the coefficients of the fixed covariate effects and random effects, if present. Moreover, a relatively small number of so-called hyperparameters governing automatically. In Bayesian modelling, we specify so-called prior probability distributions for the components and parameters to estimate. In general, prior distributions allow incorporating expert knowledge and can help stabilise estimated models when these have a very complex structure and/or when data are very noisy. Bayes' famous Theorem then tells us how we can construct the posterior distributions (densities, expected values, etc.), i.e., parameter estimates and precise probabilistic uncertainty statements, by confronting prior information with observed data. More precisely, in Bayesian statistics, the data are typically assumed to be generated from a probability distribution with probability density $\pi(\cdot | \theta)$, which depends on some parameter vector $\theta$, which is itself distributed according to some prior distribution $\pi(\theta)$ set by the modeller. The object of interest is the posterior distribution $\pi(\theta | \text{data})$, which can be found using Bayes' Theorem as

$$\pi(\theta | \text{data}) = \frac{\pi(\text{data} | \theta)\pi(\theta)}{\pi(\text{data})} \propto \pi(\text{data} | \theta)\pi(\theta),$$
where $\pi(\text{data} | \theta)$ is called the likelihood function and $\pi(\text{data}) = \int \pi(\text{data} | \theta)\pi(\theta)d\theta$ is the marginal likelihood. The symbol $\propto$ refers to equality up to a proportionality constant, here $\pi(\text{data})$, which does not depend on the parameters $\theta$. Therefore, the Bayesian mechanism allows us to move from the specification of “data distribution given the model parameters and their prior information” to the estimation target corresponding to the “distribution of model parameters given the data and prior information”. The fitted posterior distribution $\pi(\theta | \text{data})$ can then be exploited to make inference (e.g., on model parameters and their uncertainty), derive any model-based diagnostics, and draw practical conclusions. With Bayesian hierarchical models (see, e.g., Diggle and Ribeiro, 2007; Banerjee et al., 2014), as considered here, the setting is more complex as we have an extra layer with latent parameters, $x$, to infer simultaneously, which can sometimes be very high-dimensional. In our case, the latent parameters comprise all covariate coefficients and random effects, as well as the random point process intensities $\Lambda_j(s_i)$; recall (2). Bayes' formula can also be applied in this context and we get

$$ \pi(x, \theta | \text{data}) = \frac{\pi(\text{data} | x, \theta)\pi(x | \theta)\pi(\theta)}{\iint \pi(\text{data} | x, \theta)\pi(x | \theta)\pi(\theta)dxd\theta}, \quad (3) $$

where $\pi(\cdot)$ denotes a generic density function for model components. While the joint posterior $\pi(x, \theta | \text{data})$ may be of interest, it is common to focus on marginal quantities of the form

$$ \pi(\theta_k | \text{data}) = \iint \pi(x, \theta | \text{data}) dxd\theta_{-k}, \quad (4) $$

$$ \pi(x_l | \text{data}) = \iint \pi(x, \theta | \text{data}) dx_{-l} d\theta = \int \pi(x_l | \theta, \text{data}) \pi(\theta | \text{data}) d\theta, \quad (5) $$

where $\theta_k$ denotes the k-th element of the hyperparameter vector $\theta$, and $\theta_{-k}$ is the vector $\theta$ with its k-th element removed, and so forth. Computing quantities (4) and (5) may be very difficult, especially when there are numerous latent variables, and model simplifications as well as numerical approximations are required. In what follows, we will work under the framework of Bayesian models where the latent vector $x$ has a multivariate normal distribution, and the data are conditionally independent given this latent vector $x$ and hyperparameters $\theta$. As explained in the next Section 5.3, this allows us to exploit a fast and accurate approximate Bayesian technique known as INLA to estimate these posterior quantities and make inference.

Through the specification of priors, the Bayesian paradigm allows us to resolve potential parameter identifiability issues by naturally integrating constraints in the prior distributions. For example, if some covariates have a tendency towards collinearity (e.g., if they are strongly correlated, hence providing redundant information—a common issue in landslide susceptibility and hazard studies), then the prior structure can reduce numerical instabilities and keep the model and its parameters well identifiable from the data. Therefore, the specification of priors is crucial in models where the number of unknown parameters and/or latent random variables to infer is large compared to the observed sample size.

Here, we will specify different types of prior distributions for distinct model components (i.e., fixed effects and random effects). In particular, we assume that the continuous covariates have been standardised to have mean 0 and variance 1, such that the importance of estimated coefficients can be interpreted and compared more easily, and estimated coefficients will tend to be significant if they are at least moderately large in absolute value. Since only a relatively small number of such coefficients are estimated from a quite large number of observations, these coefficients are usually well estimated in the posterior model. Therefore, it makes sense to fix a moderately informative prior distribution with fixed prior variance for such coefficients, and the exact value of the fixed prior variance has little influence on the posterior model. For example, in our models, we specify a moderately informative normal distribution with mean 0 and variance 1

independently for the global intercept and all covariate coefficients. This implies that, a priori, the probability of having an absolute covariate coefficient value $\beta$ larger than 2 is less than 5 %. However, if the data provide clear evidence that the true coefficient is higher, the final posterior estimate of the coefficient can still be much larger without any difficulty.

As for latent random effects composed of numerous random variables, a general principle resulting from the law of parsimony known as “Occam’s razor” is to avoid overly complex models by constructing prior distributions that penalise the model complexity (Simpson et al., 2016). If the stochastic behavior of such complex model components is not properly controlled by suitably chosen prior distributions, this can lead to overfitting and poor estimation and prediction performances. To penalise model complexity, a strategy is to use informative priors, which shrink complex model components towards simpler reference models, making sure that they can be reliably estimated. This general principle has been formalised recently and leads to an approach based on so-called Penalised Complexity (PC) priors (Simpson et al., 2016). We will make systematic use of PC priors in our implementation. For instance, the reference model for a latent spatial random effect (i.e., the LSE) is simply taken to be its absence, such that the prior distribution of the precision parameter of this effect is designed to avoid overly large spatial variability. In our baseline models, the random effect components are replaced by their reference distribution, i.e., they are absent. Even with principled approaches such as the PC priors, there is no general rule providing exactly the “best” choice for the prior distribution, such that some subjectivity remains. In general, if we observe rather large credible intervals of hyperparameters in the estimated model, or if numerical instabilities arise during the estimation of the model, it is recommended to study the sensitivity of estimated models to prior choices before taking the final choice of how to fix prior distributions.

### 5.3. Latent Gaussian modelling approach, and the R-INLA library

Within the last decade since its publication in 2009, the Integrated Nested Laplace Approximation (INLA) method (Rue et al., 2009) has become the most popular tool for Bayesian spatial modelling thanks to its implementation in the R-INLA library (Lindgren et al., 2015, http://www.r-inla.org/) of the statistical software R (R Core Team, 2014). The general framework of INLA is that of Bayesian latent Gaussian models, which has found widespread interest in a diverse range of applications (Lombardo et al., 2018a; Opitz et al., 2018; Krainski et al., 2018; Moraga, 2019). Essentially, the data are assumed to follow a “well-behaved” distribution function, and to be conditionally independent given some latent (multivariate) Gaussian random effects; recall the framework introduced in Section 5.2 and (3). In our context, landslide counts are modelled using the Poisson distribution, whose mean is expressed on the log scale in terms of various fixed effects and latent random effects that are correlated over space and/or time. Each of the covariate coefficients simply has a normal distribution prior, independent from the others. As for hyperparameters, the prior distributions are chosen more specifically depending on the role of the parameter. More details on each of our models are given in the following sections.

The success of INLA relies on the systematic use of random effects with sparse precision (i.e., inverse covariance) matrices within this latent Gaussian modelling framework, coupled with astute analytical and numerical approximation schemes (Illian et al., 2012; Rue et al., 2009, 2017), which provide exceptional speed-up for fitting large and complex models compared to more traditional Markov Chain Monte Carlo (MCMC) methods. For details on the theory and practice of R-INLA, we refer the interested reader to the landslide tutorial paper by Lombardo et al. (2019b) and to Bakka et al. (2018).
## 5.4. Baseline LGCP models with fixed effects only

First, we consider two “baseline” models, which we call Mod1 and Mod2, where the first one is purely spatial, in the sense that it does not include any information about the structure of time intervals, whereas the second allows for time-interval-specific regression constants. In the model Mod1, we only include the spatially-indexed covariates presented in Section 4.3 in the log-intensity:

$$ \begin{aligned} \log(\Lambda_j^{\text{Mod}1}(s_i)) & = \log(|A_i|) + \beta_0 + \sum_{k=1}^{8} \beta_{k;1}^{\text{morph}} z_k^{\text{mean}}(s_i) + \sum_{k=1}^{8} \beta_{k;2}^{\text{morph}} z_k^{\text{sd}}(s_i) + \\ & \qquad \sum_{k=1}^{13} \beta_k^{\text{them}} z_k^{\text{prop}}(s_i), \end{aligned} \quad (6) $$

where $j = 1, ..., 6$ indexes the time intervals and each $s_i$, $i = 1, ..., 889$, corresponds to a different SU. This model comprises 29 covariate coefficients $\beta$ to be estimated, here separated according to the SU-wise means and standard deviations of morphometric variables obtained from the DEM, with superscript “morph” in Eq. (6), and the 13 thematic properties, with superscript “them” in Eq. (6), expressed through SU-wise proportions (Table 2). This is a purely spatial model, which assumes that the spatial intensity is the same for each of the 6 temporal intervals (here treated as independent replicates). Moreover, this model does not account for any additional spatially-correlated nor temporally-correlated unobserved effects.

Model Mod2, with intensity $\Lambda_j^{\text{Mod2}}(s)$, has the following structure:

$$ \log(\Lambda_j^{\text{Mod}2}(s_i)) = \log(\Lambda_j^{\text{Mod}1}(s_i)) + \beta_j, \quad (7) $$

where $\beta_j$, $j = 1, ..., 6$ are additional time-interval-specific intercepts. We resolve the non-identifiability of $\beta_0$ in (6) and $\beta_j$ in (7) by imposing the sum-to-zero constraint $\Sigma_{j=1}^6 \beta_j = 0$, such that estimated $\beta_j$-coefficients ($j = 1, ..., 6$) indicate how strongly the overall number of landslides in a time interval deviates from the global average measured through $\beta_0$. As indicated in §5.2, we assign independent zero-mean Gaussian priors for each regression coefficient. However, unlike the global intercept $\beta_0$ and the covariate coefficients $\beta_{k;1}^{\text{morph}}$, $\beta_{k;2}^{\text{morph}}$, $\beta_k^{\text{them}}$ where the prior variance is set to one, the additional intercepts $\beta_j$ specific to each time interval do not have a fixed prior variance; instead, we specify a PC prior for the variance of $\beta_j$ which corresponds to a 50%-probability of being below or above 1, and we then let data decide how strongly the $\beta_j$ values should be allowed to vary between time intervals.

## 5.5. Spatial LGCP with replicated spatially structured random effects

To extend the baseline models Mod1 and Mod2, we now add a spatial random effect, with 6 replicates in time, to explain variations in the landslide intensity that cannot be explained by the observed covariates, and we write $\Lambda_j^{\text{Mod3}}(s)$ for the spatial intensity in the $j$-th temporal interval in model Mod3. Our prior assumption is that the spatial random effect should differ between SUs and between different temporal events, but that it tends to be similar for SUs sharing a boundary or being close in space; recall the adjacency graph structure in Fig. 2C.

We first write the general model formula, and we then explain how we encode this prior assumption on spatial dependence for the random effect. The model may be written as

$$ \begin{aligned} \log(\Lambda_j^{\text{Mod}3}(s_i)) = \log(\Lambda_j^{\text{Mod}2}(s_i)) + W_j^{\text{Mod}3}(s_i), & \quad i = 1, ..., 889, \\ j = 1, ..., 6, & \end{aligned} \quad (8) $$

where $\Lambda_j^{\text{Mod2}}(s)$ is the baseline intensity of model Mod2 defined in (7), and where $W_j^{\text{Mod3}}(s)$ is the latent spatial effect (LSE) with 6 replicates, one for each of the temporal intervals.

We use the notation $i_1 \sim i_2$ if the SUs $i_1$ and $i_2$ are not the same but

share a boundary, and we write $nb(i) = \{\tilde{i} | \tilde{i} \sim i\}$ for the set of all the neighbouring SUs that are adjacent to the $i$-th SU, with $|nb(i)|$ indicating their number. Since the study area is contiguous, we obtain $|nb(i)| \ge 2$ for all SUs $i = 1, ..., 889$ (Fig. 2B). The model for $W_j^{\text{Mod3}}(s_i)$ is now specified by the following two conditions:

1. The spatial fields $W_{j_1}^{\text{Mod3}}(s)$ and $W_{j_2}^{\text{Mod3}}(s)$ are independent for different times $j_1 \neq j_2$.

2. The value of the spatial random effect $W_j^{\text{Mod3}}(s_i)$ at the $i$-th SU, given all the other values $W_j^{\text{Mod3}}(s_i)$, $\tilde{i} \neq i$, follows a normal distribution whose mean value corresponds to the mean of the adjacent values, i.e.,

$$ W_j^{\text{Mod3}}(s_i)|\{W_j^{\text{Mod3}}(s_i); \tilde{i} \neq i\} \sim \text{N}\left(\frac{1}{|nb(i)|}, \sum_{i \neq i} W_j^{\text{Mod3}}(s_i), \frac{1}{|nb(i)|} \tau^{\text{Mod3}}\right), \quad (9) $$

where $\tau^{\text{Mod3}} > 0$ is a precision parameter to be estimated, which controls the dependence strength between neighbouring SUs. The parameter $\tau^{\text{Mod3}}$ of the conditional spatial distributions in (9) determines the value of the variance $1/\kappa^{\text{Mod3}}$ of the unconditional effects $W_j^{\text{Mod3}}(s_i)$; internally, R-INLA implements a parameterisation using the marginal precision $\kappa^{\text{Mod3}}$ averaged over all SUs, which may be simpler to interpret in practice.

The mechanism of this model prescribes that there is less uncertainty about the random effect value in a SU if we know the values in the adjacent SUs. This prior assumption is valid when adjacent slope units have a similar behaviour, e.g., when the sliding surface at depth affects more than one SU, or when the landslide rainfall/seismic trigger has a clear dominating spatial pattern. However, the observed data can also counteract this prior assumption if necessary, i.e., in case two adjacent slope units have strongly different behaviours. This might occur in our study area, e.g., with the common case of bedding planes dipping out of, or nearly parallel to the slope (“cataclinal” slope) in one SU, and dipping into the slope (“anACLinal” slope) in an adjacent SU across a drainage line (Marchesini et al., 2015; Santangelo et al., 2015b). The inclusion of the LSE, $W_j^{\text{Mod3}}(s)$, in Mod3 induces spatial dependence within each temporal interval, but keeps the different temporal intervals independent, be they consecutive in time or not.

A priori, we assume that the LSEs have moderate absolute values and are relatively similar between adjacent SUs, such that we construct a prior model whose complexity remains moderate when compared to the baseline Mod2 in (7). To achieve this, we use a PC prior for the standard deviation parameter $sd^{\text{Mod3}} = 1/\sqrt{\kappa^{\text{Mod3}}}$, which corresponds to an exponential distribution; in mathematical notation, we assume that

$$ \text{Pr}(1/\sqrt{\kappa^{\text{Mod3}}} > u) = \exp(-\lambda u), \quad u > 0, \quad (10) $$

where $\lambda > 0$ is a penalty rate to be defined. Here we fix $\lambda$ such that $\text{Pr}(1/\sqrt{\kappa^{\text{Mod3}}} > 1) = 0.5$, i.e., we give the standard deviation parameter a 50% chance of exceeding 1 a priori.

Model (9) has a long history (Besag, 1974) in spatial statistics since it was first studied for data observed over spatial graphs, such as the adjacency graph of slope units (recall Fig. 2). Due to its parsimonious parametrisation, its close ties to nearest-neighbor prediction (which are obvious from Eq. (9)), and advantageous numerical representations of its dependence structure, it is a robust and intuitive prior model for latent spatial effects. Many alternative formulations of prior models for latent spatial effects are possible, e.g., by using Gaussian random fields defined over continuous space with separate parameters for the range of spatial dependence and the variance (Krainski et al., 2018), but the mathematical representation and practical implementation would become more technical. Therefore, we do not compare our choice of prior model to such alternatives in this work.
## 5.6. Temporal LGCP with slope-unit-based temporal random effect

Our next model, Mod4, has a similar structure as Mod3 in (8) at first sight, but we now make an assumption about the temporal dependence of SU-based random effects within the same SU, while disregarding any direct spatial relationship between adjacent SUs. Writing $\Lambda_j^{\text{Mod4}}(s)$ for the spatial intensity in the j-th temporal interval, the model formula is given as

$$ \log(\Lambda_j^{\text{Mod4}}(s_i)) = \log(\Lambda_j^{\text{Mod2}}(s_i)) + W_i^{\text{Mod4}}(t_j), \quad i = 1, \dots, 889, \\ j = 1, \dots, 6, \qquad (11) $$

where $\Lambda_j^{\text{Mod2}}(s)$ is the baseline intensity of model Mod2 defined in (7), and where $W_i^{\text{Mod4}}(t)$ is a latent temporal effect (LTE) with 889 replicates (one for each of the SUs in our study region), defined by the following three conditions:

1. The time series $W_i^{\text{Mod4}}(t_j)$ and $W_{i_2}^{\text{Mod4}}(t_j)$ are independent when considering different SUs $i_1 \neq i_2$.

2. For fixed SU $i$, the value of the temporal random effect $W_i^{\text{Mod4}}(t_j)$ at time $t_j$ with $j > 1$, given the value $W_i^{\text{Mod4}}(t_{j-1})$ at the preceding time point $t_{j-1}$, follows a normal distribution with an autoregressive structure, i.e.,

$$ W_i^{\text{Mod4}}(t_j) | W_i^{\text{Mod4}}(t_{j-1}) \sim N^{-}(\beta^{\text{Mod4}} W_i^{\text{Mod4}}(t_{j-1}), \frac{1}{\tau^{\text{Mod4}}}), \qquad (12) $$

where the temporal autocorrelation parameter $-1 < \beta^{\text{Mod4}} < 1$ and the precision parameter $\tau^{\text{Mod4}} > 0$ have to be estimated.

3. The value at the first time point $W_i^{\text{Mod4}}(t_1)$ follows a normal distribution with mean 0 and variance $1/\kappa^{\text{Mod4}} = 1/(\tau^{\text{Mod4}}(1 - (\beta^{\text{Mod4}})^2))$, i.e., $W_i^{\text{Mod4}}(t_1) \sim N^{-}(0, 1/\kappa^{\text{Mod4}})$.

Equivalently to the conditions 2 and 3 above, we can write $W_i^{\text{Mod4}}(t_j) = \beta^{\text{Mod4}} W_i^{\text{Mod4}}(t_{j-1}) + \varepsilon_j$, where the “innovations” $\varepsilon_j \sim N^{-}(0, 1/\tau^{\text{Mod4}})$ are mutually independent and independent of $W_i^{\text{Mod4}}(t_{j-1})$. Under these conditions, the variables $W_i^{\text{Mod4}}(t_j)$ are a priori stationary in time with normal distribution $W_i^{\text{Mod4}}(t_j) \sim N^{-}(0, 1/\kappa^{\text{Mod4}})$. From an interpretation perspective, the most interesting aspect is to have direct control over the standard deviation $sd^{\text{Mod4}} = \sqrt{1/\kappa^{\text{Mod4}}}$ and the autoregression parameter $\beta^{\text{Mod4}}$. Here, we proceed as with the spatial model, and we therefore use the same PC prior as in (10), setting $\lambda$ such that $sd^{\text{Mod4}}$ has 50% chance to exceed the value 1. When specifying a prior model for $\beta^{\text{Mod4}}$, we assume that consecutive events are only weakly linked to each other, and we implement a PC prior penalising absolute values of $\beta^{\text{Mod4}}$ close to 1. Our choice of prior is such that there is a 50% chance for $\beta^{\text{Mod4}}$ to exceed 0.5 in absolute value.

## 5.7. Space-time LGCP with combined spatial and temporal structures

Finally, considering both the spatial and temporal structures in the data, we construct our most complex model, Mod5, that combines explicit assumptions for the temporal dependence of random effects between consecutive inventories in time, and for spatial dependence based on spatial adjacency relations between SUs (Fig. 2C) for contemporaneous landslides, i.e., for landslides in the same time period (Fig. 3). Similar to the spatial model Mod3 in (8), we have a single parameter $\kappa^{\text{Mod5}} > 0$ that simultaneously governs the spatial variability and dependence strength of the random effects. Additionally, the parameter $\beta^{\text{Mod5}}$ controls the strength of association between consecutive time intervals, in analogy to the preceding temporal model Mod4 in (11). Therefore, this space-time model keeps a parsimonious parameterisation with only two hyperparameters $\beta^{\text{Mod5}}$ and $\kappa^{\text{Mod5}}$ to be estimated. Writing now $\Lambda_j^{\text{Mod5}}(s)$ for the spatial intensity in the j-th temporal interval in model Mod5, we define

$$ \log(\Lambda_j^{\text{Mod5}}(s_i)) = \log(\Lambda_j^{\text{Mod2}}(s_i)) + W_i^{\text{Mod5}}(s_i, t_j), \quad i = 1, \dots, 889, \\ j = 1, \dots, 6, \qquad (13) $$

where the latent space-time effect (LSTE), $W_i^{\text{Mod5}}(s, t)$, now combines the dependence relationships of the purely spatial model Mod3 in (8) and of the purely temporal model Mod4 in (11). Specifically, we assume the following structure:

1. The LSTE obeys the following temporal dynamics:

$$ W_i^{\text{Mod5}}(s_i, t_j) = \beta^{\text{Mod5}} W_i^{\text{Mod5}}(s_i, t_{j-1}) + \varepsilon_j(s_i), \quad -1 < \beta^{\text{Mod5}} < 1, \quad j > 1, \qquad (14) $$

where the spatial “innovation fields” $\varepsilon_j(s)$ are mutually independent in time and have the same distribution as the LSE $W_i^{\text{Mod3}}$ in (9), with $\tau^{\text{Mod3}}$ replaced by $\tau_e^{\text{Mod5}}$. We write $\kappa_e^{\text{Mod5}} > 0$ for the corresponding unconditional precision parameter used internally by R-INLA, in analogy with Model Mod3 in (8).

2. The field at the first time point $W_i^{\text{Mod5}}(s, t_1)$ has the same probability distribution as the LSE $W_i^{\text{Mod3}}$ in (9), but now we denote the unconditional precision parameter by $1/\kappa^{\text{Mod5}} = 1/(\kappa_e^{\text{Mod5}}(1 - (\beta^{\text{Mod5}})^2))$. This assures that the model is stationary in time for each SU with unconditional precision parameter $\kappa^{\text{Mod5}}$ for all 6 fields $W_i^{\text{Mod5}}(s, t_j)$, $j = 1, \dots, 6$.

The prior distribution of the precision parameter $\kappa^{\text{Mod5}}$ is fixed as in the spatial model Mod3 in (8), and the prior of the temporal autocorrelation parameter $\beta^{\text{Mod5}}$ is fixed similarly to the temporal model Mod4 in (11).

In the spatio-temporal model Mod5, the data can determine if landslide counts, not well explained by the observed covariates, tend to be similar in space between nearby SUs: (i) in case of small mass movements, which separately affect different SUs; or (ii) because of single landslide bodies, whose areal extents affect more than one SU. Similarly, the temporal component of the model captures the effect, not explained by the observed covariates, that lead to landslide counts being similar through time between consecutive events for the same SU. In particular, if estimated spatial and temporal dependencies are non-negligible, then the average number of landslides in a SU for a given temporal event is often related to the average number for SUs that are located “close” in space, and for all the events “close” in time. Therefore, this model can learn about clustering in space and persistence in time in the structure of the landslide-triggering mechanism. A general understanding of this concept could be practically translated in study areas where, due to orographic conditions, the occurrence of critical precipitation amounts may always tend to arise in the same, relatively confined area in space. The model Mod5 in (13) could capture this trigger structure through spatial dependence (confined area) and temporal dependence (same area for different events), even if no observed data are available for the relevant precipitation events. The latter can be a single trigger in case of event-based inventories or the cumulative effect of several triggers for inventories associated with a large time span. This assumption can be valid when storms have a clear spatial pattern characterised by a transition in precipitation regimes from the “epicentre” to the peripheries of the cloudburst (Lombardo et al., 2018a). However, the same assumption may not hold in our study area because of its size (less than $10 \times 10 \text{ km}^2$) and the absence of significant orographic gradients. An alternative explanation that may relate more closely to our study case could be that the spatial dependence captures the unknown effect of land use or land cover, driven, e.g., by changing economy and agricultural practices, or the dependence induced by large landslides destabilising more than one SU.
## 5.8. Implementation and model validation

To assess the models' performance and their interpretability from an explanatory perspective, we chose to implement an initial modelling stage where we fitted the baseline LGCPs (Mod1 and Mod2) and the three extensions featuring the latent fields (Mod3, Mod4, and Mod5) using 100% of the dataset. Subsequently, we assessed the predictive performance of each model. To do this, we implemented two separate cross-validation (CV) schemes.

The first approach is a spatial 10-fold CV, whereby we randomly leave out 10% of the dataset for testing and keep the complementary 90% for fitting, and we repeat this procedure for 10 (non-overlapping) test sets. The proportion of held-out data may seem relatively small in comparison to commonly used values in the literature (e.g., 30% to 50%, see Reichenbach et al., 2018), but we consider it as a sensible value to allow our complex models to learn about the spatial structure from the training data. We constrained the random extractions to avoid any SU to be sampled more than once, and we used the same SUs across time intervals. In other words, the combination of the ten CV complementary subsets reproduces the original dataset. Notably, additional CV routines could have been implemented. For instance, Goetz et al. (2011) presented a different type of spatial cross-validation technique where the random subset is selected as geographically stratified and disjoint sub-areas (Brenning et al., 2012). Specifically, they used a k-means clustering algorithm (see also, Ruß and Brenning, 2010) to randomly select coherent portions of the geographic space. This procedure is particularly useful in the context of pixel-based susceptibility models where selecting at random few pixels may result in a negligible change in the data. However, in our case, since the geographic space is partitioned into SUs, the spatial dependence—that one should all the more account for in grid-cell models—is less pronounced. Furthermore, because our spatial models consider the adjacency among SUs, removing continuous portions of the landscape would have a large impact in the model estimates. For this reason, we adopted here a purely spatial CV scheme, whereby SUs at subsampled at random.

The second CV approach is a temporal “leave-one-out” procedure whereby six models are fitted, each one calibrated on five time intervals and tested on the remaining one, regardless of the temporal position of the time periods used, i.e., the fitted model does not account for possible landslide heritage effects (Samia et al., 2017a, 2017b).

We assessed the accuracy of the estimated landslide intensities both for the fit and cross-validation procedures by comparing observed and predicted landslide counts for each model, and for each temporal interval. Here, predicted counts for SU i and time interval j are defined as the posterior mean $\hat{\Lambda}_j(s_i)$ of the corresponding intensity model $\Lambda_j(s_i)$. We also maintained a link to the most common landslide susceptibility model, i.e., the probability of observing at least one landslide in a given SU at a given time (Chung and Fabbri, 1999; Guzzetti et al., 1999; Van Westen et al., 1999). Thanks to the assumed Poisson distribution of landslide counts, the fitted intensity $\hat{\Lambda}_j(s_i)$ for the i-th SU and j-th temporal interval can be converted into the landslide susceptibility $S_j(s_i)$ via the following equation:

$$S_j(s_i) = 1 - e^{-\hat{\Lambda}_j(s_i)} \qquad (15)$$

Notice that according to our definition in (15), the susceptibility increases to one as the SU area increases to infinity. Therefore, our approach describes the “absolute” rather than “relative” susceptibility/ intensity, which is sensible because it is closer to the actual total destructiveness of landslides in a given slope unit. This is in line with the notion of “hazard”. For slope planning and management, authorities indeed need to know the total budget to spend to stabilise certain slopes, not the relative budget.

Hence, we computed performance metrics and visual model diagnostics for models Mod1 to Mod5, both in terms of landslide intensity and susceptibility, for each of the two CV schemes. Specifically, we

consider two metrics, the Area Under the Curve (AUC, Fawcett, 2006) for susceptibilities, and $\chi^2$-statistics for intensities. The $\chi^2$-statistic is defined as the sum of squared differences between observed and predicted counts, divided by the corresponding estimated variances; here, the predicted landslide count and its variance for a given SU are both equal to the corresponding estimated intensity. We also visually compare observed and predicted count values.

## 5.9. Classification of intensity and susceptibility estimates

In the statistically-based landslide susceptibility literature, there is no agreement on how to reclassify and show in map form the continuous spectrum of probability values that result from a classification model into few meaningful classes (Reichenbach et al., 2018); a fundamental step for a susceptibility zonation to be used in practical applications (Guzzetti et al., 2000; Galli et al., 2008). The simplest option is to divide the entire probability range [0, 1] into two classes of predicted “stable” and “unstable” conditions; with the stable conditions predicted not to have landslides, and the unstable conditions predicted to have landslides. Even for this simple twofold division, several approaches are found in the literature with authors arbitrarily setting the probability cutoff. For presence-absence balanced datasets, examples exist where the cutoff is set to 0.5 (Dai and Lee, 2002). At times, this seems not to be supported in the text with a specific explanation (e.g., Süzen and Kaya, 2012), or it is justified because it corresponds to the mean value between the two extremes of the theoretical probability range (e.g., Lombardo et al., 2016a). The approach is problematic, because it sets the cutoff where the model is most uncertain (Rossi et al., 2010a; Reichenbach et al., 2018). Frattini et al. (2010) pointed out that the choice of a cutoff value depends on the proportion of the presence-absence cases, leading to (much) smaller probability cutoff in datasets with a larger prevalence. For instance, Van Den Eeckhaut et al. (2006) reported an upper cutoff of 0.0012 for the high susceptibility class in a study case in the Flemish Ardennes where landslides were rare. In an attempt to address the issue, Castro Camilo et al. (2017) proposed to select a cutoff value that maximises the model accuracy, obtained testing thousands of probability values from the estimated susceptibility.

The problem becomes even more complex when the classification scheme involves more than two classes, typically three to seven (Reichenbach et al., 2018). A number of authors have used quantiles to segment the continuous probability estimates into discrete susceptibility classes, e.g., from “low” to “high” susceptibility. As an example, Lombardo and Mai (2018) used 2.5, 25, 50, 75, and 97.5 percentiles (i.e., where the three central values are motivated from the classical “boxplot”). Other authors segmented their probability estimates by counting the number of observed landslides in each probability class. As an example, Petschko et al. (2014) counted the number of landslides in each susceptibility “bin”, and set probability thresholds corresponding to 5%, 25%, and 70% of the total observed landslides. An alternative approach was proposed by Chung and Fabbri (2003), who ranked their probability estimates based on an “effectiveness ratio”, i.e., the ratio between the proportion of total landslide area, $A_{LT}$ in each susceptibility class and the proportion of the susceptibility class in the study area. Guzzetti et al. (2006b) adopted this approach to show the result of their susceptibility zonation for the Collazzone study area. The mentioned approaches have their rationale, but examples exist in the literature for which there is no clear justification for the adopted classification system, especially when the aim is to produce a predictive map rather than measuring the model performance. For instance, Ayalew and Yamagishi (2005) reported that slicing the probability domain into equally-spaced bins was not optimal in their case, and suggested using bins equal to the standard deviation of the probability. Similarly, Pourghasemi et al. (2012) used the natural break tool in ArcGIS® without providing an explanation for this choice. This is the current, controversial, and unclear state of play in the landslide susceptibility
literature. We further note that all these methods fail to recognise that probability values near 0.5 reveal the inability of the model to determine if a mapping unit (e.g., a SU) is stable or unstable, and do not represent “medium” or “intermediate” susceptibility conditions, or levels (Rossi et al., 2010a; Reichenbach et al., 2018).

In this work, we propose an innovative strategy to classify and rank the intensity and the susceptibility estimates, and theirs associated uncertainties, provided by our five models. We adopt the following four-classes ranking scheme:

* Clearly stable SUs: intensity ≤ 0.05, equivalent to susceptibility ≤ 0.05. This class corresponds to an intensity/susceptibility range, where the model predicts a lack of landslide occurrences with high probability. More precisely, for the SUs in this class, the probability of having no landslide is estimated to be more than 95%.

* Uncertain Type 1: 0.05 < intensity ≤ 1, equivalent to 0.05 < susceptibility ≤ 0.63. This class characterises SUs that have a probability of landslide occurrence being greater than 5%, while their expected landslide count is less than one (on average).

* Uncertain Type 2: 1 < intensity ≤ 3, equivalent to 0.63 < susceptibility ≤ 0.95. This class characterises SUs that have a probability of landslide occurrence being less than 95%, while their expected landslide count is more than one (on average).

* Clearly unstable SUs: intensity > 3, equivalent to susceptibility > 0.95. This class corresponds to an intensity/susceptibility range, where the model predicts landslide occurrences with high probability. More precisely, for the SUs in this class, the probability of having at least one landslide is estimated to be more than 95%, and the expected landslide count is more than 3 per SU.

This new classification is applied and discussed below in §6.7.

# 6. Results

In this section, we present the results of our modelling effort. First, in §6.1, we briefly outline the results of our simplest, baseline model Mod1 in (6), showing the model intensity and susceptibility estimates. Model Mod1 is the closest LGCP counterpart to more “traditional” susceptibility models, and it represents a good reference against which

to compare the more complex models, which we do in §6.2. This is followed by the presentation of a rigorous comparative assessment of the models' goodness-of-fit (i.e., “within-sample”) performance (§6.3), of the latent temporal (§6.4) and linear covariates (§6.5) effects, of the models' predictive (i.e., “out-of-sample”) skill assessed by two different cross-validation schemes (§6.6), and of what we consider our best intensity-susceptibility model (Mod5) (§6.7). Lastly, in §6.8, we provide information on the computational requirements to fit our models.

## 6.1. Baseline intensity and susceptibility estimates

For each of the 889 SUs that partition our study area, Fig. 4 shows, in map form, the fitted estimates obtained by model Mod1 exploiting the morphometric and thematic covariates listed in Table 2, and without considering any spatial or temporal latent dependency in the data.

Fig. 4A shows the spatial distribution of the landslide intensity, i.e., of the number of landslides that, on average, can be expected in each SU, based on the used covariates. More precisely, the map shows the temporal average, from T1 to T6 (see Fig. 3), of the modelled intensities for each time interval, i.e., their posterior means. We note that model Mod1 estimates up to 6.4 landslides per SU; a figure that is rather low in the light of 35 observed landslide count values of 10 or more, and a maximum number of landslides given as 25, in a SU in the multi-temporal inventory (see Figs. 1 and 3). By taking into account that a Poisson distribution with mean 6.4 has a variance of 6.4, we can see this as a limitation of model Mod1 in realistically representing the largest landslide counts. Fig. 4B shows the spatial distribution of landslide susceptibility, i.e., the propensity (proneness) of each SU to generate landslides, based on the considered local terrain conditions (Table 2). In the map, the susceptibility estimates were obtained from the intensity estimates of Fig. 4A using Eq. (15). In the intensity (Fig. 4A) and the susceptibility (Fig. 4B) maps, the histograms show the frequency distributions of the modelled intensity and susceptibility estimates. Visual comparison of the two histograms reveals that the distributions are significantly different; with the distribution of the intensity estimates positively skewed, and the distribution for the corresponding susceptibility estimates more uniform. This was expected, as the number of SUs that can generate a large or very large number of landslides is limited in

<figure><img src="image_8.png" /><figcaption>Fig. 4. Maps show (A) the fitted baseline landslide intensity, (B) susceptibility, and (C) unified reclassification estimated based on model Mod1 constructed using the morphometric, geologic, and bedding attitude variables listed in Table 2. In (A) and (B), histograms show the frequency of the modelled intensity and susceptibility values. In (C), pie-chart shows the percentage of SUs falling into one of the four considered classes.</figcaption></figure>
the study area, whereas the number of SUs that can generate landslides, i.e., that are potentially “susceptible” to landslides, is large and geographically distributed.

Lastly, Fig. 4C portrays a joint landslide intensity–susceptibility map prepared adopting the four-class ranking scheme proposed in §5.9, which summarises the other two maps. In the map, out of the 889 SUs, 82 (9.2 %) are classified as “clearly stable” (blue), and 47 (5.3 %) as “clearly unstable” (red). Overall, the terrain estimated to be “clearly stable” covers 1.9 km², 2.4 % of the study area, and the terrain estimated to be “clearly unstable” covers 13.3 km², 16.9 % of the study area. In the “clearly stable” SUs, the estimated landslide intensity is expected to be ≤0.05, and the corresponding susceptibility also ≤0.05, i.e., very low. Conversely, in the “clearly unstable” SUs, the landslide intensity is expected to be > 3, i.e., three or more expected landslides, and the susceptibility is > 0.95, i.e., very high. Further inspection of Fig. 4C reveals that the majority of the SUs (58.8 %), covering 28.5 km², 36.1 % of the study area, is considered of “Uncertain Type 1”, followed by 26.7 % of “Uncertain Type 2”, covering 35.2 km² or 44.7 % of the Collazzone area. In these areas, on average, less (respectively more) than one landslide is expected in each SU.

## 6.2. Intensity and susceptibility estimates of complex LGCP models

To highlight the higher flexibility and the better performance of the more complex models Mod3, Mod4, and Mod5, featuring spatial, temporal, and spatio-temporal latent effects, respectively, we show in Fig. 5 the relative intensity maps, i.e., the ratio of the intensity estimates obtained for each time interval (T1 to T6) obtained by models Mod3, Mod4, and Mod5, to the baseline model Mod1, and in Fig. 6 the corresponding relative susceptibility maps, showing the ratio of the susceptibility estimates obtained by the three complex models to the baseline model Mod1.

Interestingly, the values portrayed in the maps shown in the two Figs. 5 and 6 represent the factors by which the intensity and susceptibility baseline estimates should be multiplied to account for the space, time, and space-time dependencies. More precisely, the maps in the 18 panels of Fig. 5 show $\hat{\Lambda}_j^{\text{Mod}3}(s_i)/\hat{\Lambda}_j^{\text{Mod}1}(s_i)$, $\hat{\Lambda}_j^{\text{Mod}4}(s_i)/\hat{\Lambda}_j^{\text{Mod}1}(s_i)$ and $\hat{\Lambda}_j^{\text{Mod}5}(s_i)/\hat{\Lambda}_j^{\text{Mod}1}(s_i)$, estimated from each respective model for the different time intervals $j = 1, ..., 6$. Similarly, using the intensity–susceptibility conversion Eq. (15), Fig. 6 portrays $[1 - \exp\{-\hat{\Lambda}_j^{\text{Mod}3}(s_i)\}]/[1 - \exp\{-\hat{\Lambda}_j^{\text{Mod}1}(s_i)\}]$, $[1 - \exp\{-\hat{\Lambda}_j^{\text{Mod}4}(s_i)\}]/[1 - \exp\{-\hat{\Lambda}_j^{\text{Mod}1}(s_i)\}]$, and $[1 - \exp\{-\hat{\Lambda}_j^{\text{Mod}5}(s_i)\}]/[1 - \exp\{-\hat{\Lambda}_j^{\text{Mod}1}(s_i)\}]$, $j = 1, ..., 6$,

<figure><img src="image_9.png" /><figcaption>Fig. 5. Estimated intensity ratios (IRs) for the three complex models Mod3, Mod4, and Mod5, compared to the baseline model Mod1, for each temporal interval (T1 to T6) (see Fig. 3). The shown values, $\hat{\Lambda}_j^{\text{Mod}3}(s_i)/\hat{\Lambda}_j^{\text{Mod}1}(s_i)$, $j = 1, ..., 6$, express the factor by which the intensity estimates for model Mod1 (shown in Fig. 4A) have to be multiplied to get the estimated intensity for models Mod3, Mod4, and Mod5, respectively. Color bar is uniform but limited to values greater than or equal to a factor of five for graphical purposes. Small graphs show density of intensity ratios (IRs) for each model and temporal interval. Note that x-axes and y-axes in the individual graphs cover different ranges. See text for explanation.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Fig. 6.** Estimated susceptibility ratios (SRs) for the three complex models Mod3, Mod4, and Mod5, compared to the baseline model Mod1, computed for each temporal interval (T1 to T6). The reported values, $\hat{\Lambda}_j^{\text{Mod3}}(s_i)/\hat{\Lambda}_j^{\text{Mod4}}(s_i)$, $j = 1, ..., 6$, express the factor by which the susceptibility estimates for Mod1 (shown in Fig. 4B) have to be multiplied to get the estimated susceptibility for Mod3, Mod4, and Mod5, respectively. Color bar is uniform but limited to values greater than or equal to a factor of five, for graphical purposes. Small graphs show density of susceptibility ratios (SRs) for each model and temporal interval. Note that x-axes and y-axes in the individual graphs cover different ranges. See text for explanation.</figcaption></figure>

respectively. In this way, the different maps highlight the similarities and the differences in intensity and susceptibility with respect to model Mod1, facilitating the interpretation of the performance of the complex models Mod3, Mod4, and Mod5. To facilitate the visual comparison of the patterns of the estimated latent effects, when preparing the maps we choose different colour bars valid for each map in the two Figs. 5 and 6. This was obtained saturating the colour scales at a factor of 5 for the Intensity ratios (IRs), and at a factor of 3 for the Susceptibility ratios (SRs). To further help the comparison, in each map we show the probability density distributions of the estimated values.

For each of the complex models Mod3, Mod4, and Mod5, the estimated intensity and susceptibility ratios strongly differ from one, in several areas and time intervals, suggesting that the inclusion of latent random effects in these complex models is necessary to capture the large variations in intensity and susceptibility across space and time. In other words, these variations cannot be explained solely by the available covariates included in model Mod1. The higher flexibility of the random effect models may thus improve the goodness-of-fit and prediction skills. This will be investigated in more detail in §6.3 and §6.6. We further note a clear resemblance between the intensity and susceptibility ratios of the spatial and spatio-temporal models, namely Mod3 and Mod5. This hints at a greater effect of the spatial dimension

with respect to the temporal dimension in explaining the known distribution of landslides (Fig. 1).

### 6.3. Models' goodness-of-fit

For each combination of SU and time interval, Fig. 7 compares the observed to the estimated landslide counts and presence-absence information. Visual inspection of Fig. 7 reveals a clear pattern where the baseline models, Mod1 and Mod2, strongly underestimate the observed counts larger than one, whereas they often tend to largely overestimate the zeros (i.e., no landslides). By contrast, the complex models, Mod3, Mod4 and Mod5 that account for space, time, and space-time dependencies, respectively, through a relatively large number of latent variables, clearly improve the goodness-of-fit, with points aligned closer to the diagonal, the latter corresponding to a perfect fit. The same pattern can be numerically followed in terms of $\chi^2$, where the highest value is associated with Mod1 and it decreases as the model complexity increases to Mod5, with Mod3 and Mod5 having best scores, and Mod4 showing an intermediate score between Mod1/Mod2 on the one hand and Mod3/Mod5 on the other hand. The Receiver Operating Characteristic (ROC) curve and the Area Under the Curve (AUC) computed for the five models show a similar situation where Mod1 is the weakest,
<figure><img src="image_11.png" /><figcaption>**Fig. 7.** Within-sample performance summary of each model. Fitted landslide counts plotted against observed counts for the five models. For each observed count, we display boxplots of fitted counts to summarise their distribution across all SUs. Each panel also reports the associated $\chi^2$ values. Bottom right panel shows the results in a susceptibility framework, plotting the ROC curve for each model and summarising the within-sample goodness-of-fit performance with the corresponding AUC values. The higher the AUC value, the better the model fit.</figcaption></figure>

followed by Mod2, whose performance is slightly better due to the
contribution of the multiple temporal intercept. Overall, Mod3, Mod4,
and Mod5 perform much better, improving the two baseline results

from acceptable to outstanding, according to Hosmer and Lemeshow (2000). While it is natural that more complex models (e.g., Mod3, Mod4, Mod5) show a better fit to the data used for estimation, in
Section 6.6 we further corroborate the superior performance of these models through diagnostics based on cross-validated data.

## 6.4. Temporal effects

To investigate the temporal dynamics driving landslide occurrence in our study area, we now focus on model Mod4 in (12), in which we decompose the temporal effect into global multiple intercepts (with one coefficient for each time interval), assumed to be a *priori* independent across time, and latent temporal effects (LTEs) for each SU, assumed to be driven by an autoregressive temporal dependence structure (Blangiardo and Cameletti, 2015; Opitz, 2017). While the multiple intercepts are constant in space and capture abrupt changes in the overall landslide intensity over time (e.g., due to triggers of different magnitudes), the LTE is designed to capture local, SU-specific changes that are smoother in time, and we thus make the assumption that the LTE carries information about “clustering” and “repellency” effects in each SU.

The plot in Fig. 8A shows the posterior distribution of the multiple intercepts.

Panel B shows posterior means of the latent temporal effects (LTE) estimated for all 889 slope units in the study area. Strict temporal clustering effects (with the LTE increasing monotonically) are shown in red (10 SUs), and strict temporal repellency effects (with the LTE decreasing monotonically) are shown in blue (34 SUs). Grey lines are the remaining SUs for which a strict classification cannot be made (845 SUs).

Inspection of the plot reveals a sudden increase of the multiple intercept during T2 (1941–1954). This is the result of a severe regional rainfall event that hit Central Italy in December 1937 (Reichenbach et al., 1998), resulting in numerous landslides in the Collazzone area, which was captured in the multi-temporal inventory interpreting aerial photographs taken in 1941. In this sense, the multiple intercepts carry the strength of the overall effect of the triggers in the six time intervals. Conversely, the LTE captures more localised effects in each SU, estimating the relation between landslide counts in a given time interval and the number of landslides in the following time interval.

In Fig. 8B, we show the temporal evolution of the posterior means of the LTE for the 889 SUs in the study area, in a single plot. Inspection of the plot reveals that most of the SUs (845, i.e., 95.0 %) exhibit erratic temporal trends, with LTE increasing or decreasing “randomly” in time (grey lines in Fig. 8B). Closer inspection of the plot reveals that (i) a small number of SUs (34, i.e., 3.8 %) exhibits a monotonically decreasing trend of the LTE (blue lines in Fig. 8B), and (ii) an even smaller number of SUs (10, i.e., 1.2 %) exhibit a monotonically increasing trend of the LTE (red lines in Fig. 8B). From a geomorphological perspective, the first group encompasses SUs characterised by landslide temporal “repellency”, where the presence of a landslide in a time period has hampered the occurrence of new landslides in the future periods in the same SU, whereas the second group encompasses SUs characterised by temporal “clustering”, where new landslides have continually followed previous landslides in the same SU, for the entire considered period (T1–T6). The later result agrees with the findings of Samia et al. (2018) who have identified a “temporal path dependency” of new landslides on pre-existing landslides in the Collazzone area. Interestingly, the temporal response of landslide path dependency identified by Samia et al. (2018) disappears after about 10–15 year in the study area, i.e., within most of the time periods considered in this study. This explains the reduced number of SUs characterised by distinct temporal “clustering” found in this study.

## 6.5. Effects of covariates

Fig. 9 shows a summary of the posterior distribution of all the estimated regression coefficients that appeared to be significant in at least one of the five models, where significance is measured using 95% credible intervals; if an interval does not contain 0, the corresponding

covariate is considered significant. In the plots, we also show the estimated coefficients for Mod1 and Mod2 to highlight an issue common to all regression models in which residual dependence is not accounted for appropriately. In fact, the 95% credible interval of the regression coefficients estimated for Mod1 and Mod2 (with no latent effects included) is narrower than the credible intervals for the models with the LSE (orange), LTE (yellow) and LSTE (brown). This is a result of the model structure, where the simpler models (Mod1 and Mod2) are overconfident of the information carried by the observations, whereas the models that incorporate spatial (Mod3), temporal (Mod4) and spatio-temporal (Mod5) dependencies better represent SU-wise overdispersion, and thus estimate more realistic credible intervals. In our case, the differences are small, and the pattern shows that the five models assign analogous posterior mean values to each covariate, both

<figure><img src="image_12.png" /><figcaption>Fig. 8. Panel A shows multiple intercepts for each time interval, from T1 to T6, for model Mod4. Dots show posterior intercept means; bars show pointwise 97.5 and 2.5 posterior percentiles bracketing 95 % pointwise credible intervals.</figcaption></figure>
<figure><img src="image_13.png" /></figure>

(caption on next page)
**Fig. 9.** Regression coefficients of all covariates, whose coefficients were significant in at least one of the fitted models, Mod1 (pink), Mod2 (blue), Mod3 (orange), Mod4 (yellow), and Mod5 (brown). Dots show posterior means; bars show pointwise 97.5 and 2.5 posterior percentiles bracketing 95 % pointwise credible intervals. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

in absolute value and sign. We consider this an evidence of the overall goodness-of-fit of the different models, in the sense that the estimated coefficients are relatively robust to model misspecification, e.g., when disregarding spatial and temporal dependencies in Mod1 and Mod2.

Determining whether a predictive model is useful in practice depends, among other factors, on the interpretability of the estimated covariates' effects. Out of all the 29 covariates used to construct the models (see Table 2), 15 were significant in at least one model (see Fig. 9). As the estimated regression coefficients were similar across the five models, we now provide a unique interpretation. The *Mean Slope* and *Mean TWI* variables gave the strongest contribution to the models, with coefficients much larger (in absolute value) than the coefficients of all other covariates. Both variables contributed to increase the landslide intensity (hence the susceptibility) in all the SUs. From a geomorphological perspective, terrain slope controls the balance of the retaining and the destabilising forces acting on a slope (Taylor, 1948; Wu and Sidle, 1995; Donnarumma et al., 2013), and in many areas and for landslides of the slide type (Hungr et al., 2014) like the one prevalent in the Collazzone study area, terrain slope and associated parameterizations (e.g., mean slope, slope range, standard deviation of slope) are known to be positively (albeit not necessarily linearly) correlated to the presence and abundance of landslides, and hence to landslide susceptibility (Carrara et al., 1991, 1995a; Fabbri et al., 2003; Budimir et al., 2015; Lombardo and Mai, 2018). TWI is a proxy for the ability of a given area to retain surface water as a function of the terrain gradient and the upslope contributing area, favouring infiltration and the increase of the pore water pressure at depth, and, hence, slope instability (e.g., Yilmaz, 2009; Cama et al., 2017).

Curvature primarily controls convergence and divergence of overland flows which is often linked to slope stability (e.g., Ohlmacher, 2007). Here, laterally-concave *Planar Curvature* is estimated to contribute to landslide-prone conditions, whereas mean upwardly-concave *Profile Curvature* conditions and their variability within a SU increase the expected number of landslides. This effect is exacerbated by the standard deviation of *Elevation* which is a proxy for terrain roughness. Our five models concurred that an increase in the standard deviation of *Elevation* within a SU contributed to an increase in the estimated number of landslides (i.e., a larger intensity), and hence to a larger susceptibility. The mean *Relative Slope Position (RSP)* was significant only for Mod1 and Mod2, although larger RSP values contributed to increasing the estimated intensity, for all five models. The RSP is a continuous index which essentially assigns 0 to lowland and flat areas, and up to 1 to mountain tops. Thus, a positive regression coefficient suggests that SUs located mostly in the high portion of the local topography are more prone to landslides than SUs located chiefly in the lower part of the local topography.

Terrain aspect, jointly measured by the *Eastness* and *Northness* covariates, both expressed by their mean and standard deviations, played a significant role albeit with a small amplitude (i.e., small absolute coefficient values). According to their sign, SUs facing North or West were related to larger landslide intensities and larger susceptibility estimates. We note here that the decomposition of the terrain aspect into its two main components (*Eastness* and *Northness*) was a numerically convenient way to handle the nonlinear and cyclic exposure signal on slope stability/instability conditions. However, by decoupling the aspect into a linear combination of sine and cosine components, we lost the original interpretation of the overall effect of the aspect expressed in degrees within the [0,360°] range. To compensate for this, in Fig. 10 we show the overall reconstructed effect of the terrain aspect, for each model. Precisely, in the Figure we plot $\beta_{\text{Eastness}} \cos(\theta) + \beta_{\text{Northness}} \sin(\theta)$ as a function of $\theta \in [0,360°]$, where

$\beta_{\text{Eastness}}$ and $\beta_{\text{Northness}}$ denote the *Eastness* and *Northness* coefficients, respectively, estimated from each model. Inspection of the Fig. 10 reveals that the effect of terrain aspect on landslide intensity and susceptibility is significant, and when back-transformed to its original scale, the W-NW components mentioned above reveal a clear positive contribution to the landslide counts (i.e., landslide intensity), which changes to a negative effect when moving towards E-SE components. This was known in the study area, and depends on the geometric and geomorphological interaction between the prevalent attitude of the bedding planes that characterise the study area and the orientation and geometry of the slopes (Marchesini et al., 2015; Santangelo et al., 2015b).

We have opted for a linear influence of all the quantitative covariates in our models, even though the joint inclusion of the SU-based mean and standard deviation of the original covariates allows capturing to some extent a certain type of “nonlinearity”. Our Bayesian framework would easily allow including more general nonlinear response functions for covariates, which can be achieved through the specification of random effects, similar to the one of the time dimension; see Lombardo et al. (2018a) for an example. However, for the sake of clarity and of the focus on aspects related to the spatial and temporal random effects, we did not explore such nonlinear effects in this work.

## 6.6. Predictive performance of models

The goodness-of-fit (i.e., “within-sample” performance) of the baseline models Mod1 and Mod2 was weak in terms of fitted counts, but acceptable in terms of binary metrics. Conversely, the more complex random effect models Mod3, Mod4, and Mod5 showed outstanding within-sample performances both in terms of expected landslide counts, and fitted presence-absence probabilities (see Fig. 7). Thus, it is natural

<figure><img src="image_14.png" /></figure>

**Fig. 10.** Estimated effect (solid curves) of terrain aspect for models Mod1 (pink), Mod2 (blue), Mod3 (orange), Mod4 (yellow), and Mod5 (brown), obtained by combining the effects of *Eastness* and *Northness*. The plot shows $\beta_{\text{Eastness}} \cos(\theta) + \beta_{\text{Northness}} \sin(\theta)$ as a function of $\theta \in [0,360°]$, where $\beta_{\text{Eastness}}$ and $\beta_{\text{Northness}}$ here denote (with some abuse of notation) the *Eastness* and *Northness* coefficients, respectively, estimated from each model. Dashed lines show corresponding 95 % credible bands. Curves for models Mod1 and Mod2 cannot be distinguished. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
<figure><img src="image_15.png" /><figcaption>**Fig. 11.** Out-of-sample performance summary of each fitted LGCP model Mod1 to Mod5 (top to bottom). The plot space is divided into two. The first two columns correspond to count performance whereas the third and fourth columns report binary performance. Then the first and third columns summarise the spatial 10-fold CV results whereas the second and fourth columns summarise the temporal leave-one-out CV scheme. For each observed count distributed across all SUs, we graphically show the agreement between observed and predicted number of landslides. For each cross-validation we also report the $\chi^2$ statistics. In analogy, we graphically show the ROC curves and report the associated AUC values for each cross validation.</figcaption></figure>
to wonder whether or not the more complex LGCP models overfit the data. In case of overfitting, their predictive (“out-of-sample”) performance would be low.

To assess this, and to quantify the predictive performance of each model, we designed two cross-validation (CV) procedures as explained in Section 5.8. Because the data are spatio-temporal in nature, we considered both a spatial CV scheme and a temporal CV scheme. We measured the spatial predictive performance using a 10-fold cross-validation procedure (abbreviated *Space 10-Fold*), and the temporal predictive performance using a leave-one-out cross-validation procedure (abbreviated *Time leave-one-out*). We here briefly recall our approach. The spatial 10-fold CV consists in splitting the original dataset into 10 complementary subsets at random, each comprising 10 % of the original SUs. The model is then fitted using nine subsets (i.e., 90 % of the SUs), and subsequently validated using the left-out subset (i.e., 10 % of the SUs). The procedure is repeated by leaving out each subset. The complementary constraint ensures that every SU in the Collazzone study area is predicted exactly once for all time intervals during the CV routine. As for the temporal leave-one-out CV, we leave out the data from one of the six time intervals, then we fit the model using the five remaining intervals, and finally we validate the model using the data from the time interval that was left out. The procedure is repeated by leaving out the data from each time interval. Essentially, this corresponds to a temporal 6-fold CV scheme, where each test set consists of a single time interval.

For both CV schemes, we examined the models' performances both in terms of intensity (i.e., predicted counts) and susceptibility (i.e., predicted probability of landslide occurrence), using the same summary measures used in Fig. 7. Results are summarised in Fig. 11 where the two main vertical panels represent the intensity and susceptibility results, and the sub-columns summarise the results for the *Space 10-fold* and *Time leave-one-out* procedures. The different rows correspond to Mod1, Mod2, Mod3, Mod4, and Mod5 (from top to bottom). We opted to avoid colour coding the *Time leave-one-out* intensities by time to improve readability of the figure.

Inspection of Fig. 11 reveals that the agreement between the observed and the predicted landslide counts increases significantly from our simple models (Mod1, Mod2) to the complex models that include latent variables (Mod3, Mod4, Mod5), which are better capable of predicting the actual number of landslides in each SU. The cross-

validated $\chi^2$ measures tend to decrease strongly from Mod1/Mod2 to Mod3/Mod4/Mod5, with best values attained for Mod3 and Mod5. In Mod3 and Mod5, the match between observed and estimated counts is reasonably good, up to 14 landslides. Conversely, from this number to the maximum (25 landslides) in the dataset, the predicted counts tend to underestimate the actual observations. We underscore that this is only an indication, but not necessarily a formal proof of bad fit of the models, since predicted counts are the *expected values* of the corresponding Poisson distribution of landslide counts, which do not convey the natural variability of this Poisson distribution. This underestimation for large counts is however a specific consequence of our dataset, which comprises a large number of SUs with no or a few landslides and very few SUs, often spatially and temporally isolated, with many landslides. Therefore, it is very difficult for a model to deterministically predict such large counts from the information in adjacent SUs and time intervals, where counts are usually much smaller. Rather, a model may allow for such high counts through the uncertainty carried by the Poisson distribution of observations. Inspection of Fig. 3 reveals that the landslide dataset has only a few isolated samples larger than 10 counts, which is where the model starts to perform poorly. Also, the *Time leave-one-out* is the CV scheme that deviates most strongly from the performance obtained for the fit. An explanation is that the *Space 10-fold* removes 533 SUs per iteration (889 × 6/10 ≈ 533), whereas a temporal CV removes 889 SUs. Moreover, the temporal CV disrupts the coherence in the data more strongly, since removing one time interval removes either 50 % or 100 % of the direct temporal neighbours for a large number of data points. Therefore, it is the most challenging CV scenario we could devise. Nevertheless, we note that the overall performance shown in the susceptibility case still falls in the “excellent” class of Hosmer and Lemeshow (2000), with outstanding AUC metrics. The improvement from Mod1 to Mod5 is measured quantitatively by the models' AUCs, which justify the inclusion of the latent effects.

## 6.7. Best intensity-susceptibility predictive model—Mod5

We examined the results or our modelling effort in terms of (i) estimated landslide counts, i.e., of predicted landslide intensity (Lombardo et al., 2018a, 2019b,a), and of (ii) estimated binary presence-absence of landslides, i.e., of predicted landslide susceptibility (Brabb, 1985; Guzzetti et al., 1999; Reichenbach et al., 2018). In doing

<table><caption>**Table 3**
Summary of the unified intensity-susceptibility classes, from T1 to T6. See Fig. 3 for coverage of time intervals.</caption><thead><tr><th>Time interval</th><th>Class</th><th>Slope units (SUs) Count (#)</th><th>%</th><th>Total SU area km<sup>2</sup></th><th>%</th></tr></thead><tbody><tr><td rowspan="4">T1<br>before 1941</td><td>Clearly stable</td><td>105</td><td>11.8</td><td>3.0</td><td>3.8</td></tr><tr><td>Uncertain Type 1</td><td>514</td><td>57.8</td><td>40.3</td><td>51.1</td></tr><tr><td>Uncertain Type 2</td><td>227</td><td>25.5</td><td>25.5</td><td>32.3</td></tr><tr><td>Clearly unstable</td><td>43</td><td>4.8</td><td>10.1</td><td>12.8</td></tr><tr><td rowspan="4">T2<br>1941–1954</td><td>Clearly stable</td><td>97</td><td>10.9</td><td>3.5</td><td>4.4</td></tr><tr><td>Uncertain Type 1</td><td>409</td><td>46.0</td><td>25.6</td><td>32.4</td></tr><tr><td>Uncertain Type 2</td><td>243</td><td>27.3</td><td>26.3</td><td>33.3</td></tr><tr><td>Clearly unstable</td><td>140</td><td>15.8</td><td>23.6</td><td>29.9</td></tr><tr><td rowspan="4">T3<br>1954–1977</td><td>Clearly stable</td><td>106</td><td>11.9</td><td>3.6</td><td>4.6</td></tr><tr><td>Uncertain Type 1</td><td>511</td><td>57.5</td><td>31.7</td><td>40.1</td></tr><tr><td>Uncertain Type 2</td><td>202</td><td>22.7</td><td>26.1</td><td>33.0</td></tr><tr><td>Clearly unstable</td><td>70</td><td>7.9</td><td>17.6</td><td>22.3</td></tr><tr><td rowspan="4">T4<br>1977–1996</td><td>Clearly stable</td><td>169</td><td>19.0</td><td>6.8</td><td>8.7</td></tr><tr><td>Uncertain Type 1</td><td>602</td><td>67.7</td><td>49.8</td><td>63.1</td></tr><tr><td>Uncertain Type 2</td><td>100</td><td>11.3</td><td>17.8</td><td>22.5</td></tr><tr><td>Clearly unstable</td><td>18</td><td>2.0</td><td>4.6</td><td>5.8</td></tr><tr><td rowspan="4">T5<br>1997 (snow)</td><td>Clearly stable</td><td>153</td><td>17.2</td><td>5.7</td><td>7.2</td></tr><tr><td>Uncertain Type 1</td><td>586</td><td>65.9</td><td>45.4</td><td>57.6</td></tr><tr><td>Uncertain Type 2</td><td>129</td><td>14.5</td><td>23.5</td><td>29.8</td></tr><tr><td>Clearly unstable</td><td>21</td><td>2.4</td><td>4.3</td><td>5.4</td></tr><tr><td rowspan="4">T6<br>1998–2014</td><td>Clearly stable</td><td>100</td><td>11.3</td><td>2.6</td><td>3.3</td></tr><tr><td>Uncertain Type 1</td><td>522</td><td>58.7</td><td>37.2</td><td>47.2</td></tr><tr><td>Uncertain Type 2</td><td>186</td><td>20.9</td><td>22.2</td><td>28.1</td></tr><tr><td>Clearly unstable</td><td>81</td><td>9.1</td><td>16.9</td><td>21.4</td></tr></tbody></table>
so, we considered the models' spatial and temporal structures (Figs. 5, 6), and their fitting (Fig. 7) and predictive (Fig. 11) performances. In terms of intensity, Mod3 and Mod5 perform equally well with an analogous match between observed and estimated landslide counts, and with similar $\chi^2$ values. Conversely, when converting intensity values into their susceptibility counterparts, Mod5 performed slightly better, especially for the temporal leave-one-out scheme. As a result, we consider Mod5—which jointly accounts for spatial and temporal dependence—as our “best” overall model with increased performance compared to the other four models (Mod1 to Mod4). In addition, the temporal autocorrelation coefficient $\hat{\beta}^{\text{Mod5}}$ is estimated to be 0.68 and is significantly different from 0 (considering its 95% credible interval), which strongly suggests that Mod5 fits the data better than Mod3. Mod5 provides comparable patterns to Mod3 in terms of the predicted landslide counts over time, but also performs well in the binary (i.e., susceptibility) predictions, similarly to Mod4 and in contrast to Mod3. We conclude that Mod5 is our best model, and we select it to generate a classification summarising both the landslide intensity and susceptibility for our study area, adopting the ranking scheme proposed in §5.9. We summarise the results of Mod5, for each of the six time intervals (T1–T6), in Table 3 and Fig. 12. Inspection of results reveals some degree of temporal variability in the combined intensity-susceptibility patterns. However, the general spatio-temporal pattern remains about the same. Similarly to the outcomes of our baseline model Mod1 (§6.1), the number of SUs that Mod5 predicts capable of generating a large or very large number of landslides is limited, whereas the number of SUs that can generate landslides, i.e., that are potentially “susceptible” to slope failures, is large and geographically distributed. This is reasonable from a geomorphological and landscape evolution perspectives.

Overall, and for the entire considered period (Fig. 3), model Mod5 classifies the (relative or absolute) majority of the SUs, from 409 (T2, 46.0 %) to 602 (T4, 67.5%), as of “Uncertain Type 1” (§5.9). In each of these SUs, covering collectively between 25.6 (32.4 %) and 49.8 (63.1 %) km², the estimated landslide intensity, i.e., the expected number of landslides, is in the range (0.05,1], on average. The second class encompasses from 100 (T4, 11.3 %) to 227 (T1, 25.5 %) SUs classified as of “Uncertain Type 2”. In these SUs, covering collectively between 17.8 and 25.5 km² (22.5 to 32.3 %), the number of expected landslides is in the range (1,3], on average (Table 3). With a few exceptions, model Mod5 classifies the smallest number of SUs, from only 18 (T4, 2.0 %) to 140 (T6, 15.8 %), as “Clearly Unstable”, covering between 4.6 (5.8 %) and 23.6 (29.9 %) km², followed by from 97 (T2, 10.9 %) to 169 (T4, 19.0 %) SUs classified as “Clearly Stable”, covering between 3.5 (4.4 %) and 6.8 (8.7 %) km². The predicted classification estimates (Fig. 12, Table 3) reflect the true spatio-temporal variability of the landslide intensity-susceptibility at the spatial scale and in the temporal range considered by our modelling experiment. Such estimates are, however, subject to different sources of uncertainty, including (i) measurement uncertainty (e.g., older data are typically of lower quality); (ii) model uncertainty (our model might have been chosen differently if the data had been different); and (iii) estimation uncertainty (a fitted model will never perfectly match the true model, even if the correct family of models is known). In our case, measurement uncertainty is likely to be relatively negligible because the study area has been constantly surveyed for over 40 years and orthophotos are available for about 80 years (recall Fig. 1 and Table 2), and because we work here with SUs, which smooth out positional errors and other spatial biases. Moreover, model uncertainty is also quite small because our choice of covariates is here driven by geomorphological rather than statistical considerations, and prior distributions for hyperparameters are specified in a way to “let the data speak for themselves” while getting stable model fits. Estimation uncertainty is perhaps the most critical one, but because of our rich spatio-temporal dataset and our careful selection of prior distributions, it can also be considered relatively minor overall.

Another natural source of uncertainty arises when predicting future landslide scenarios from our fitted model: stochastic uncertainty (even if a model is known to describe landslide occurrences perfectly, future values would inevitably vary). This should be kept in mind when interpreting our results, and in particular Fig. 12 and Table 3.

<figure><img src="image_16.png" /><figcaption>Fig. 12. Intensity-susceptibility classification of the Collazzone study area, Umbria, Central Italy, based on model Mod5 constructed using morphometric, geologic, bedding attitude, and space-time dependencies. Each map corresponds to one of the six time intervals shown in Fig. 3. Pie-charts show the percentage of SUs falling into one of the four considered classes. See §5.9 for an explanation of the adopted classification and ranking scheme.</figcaption></figure>
## 6.8. Computational requirements

All models presented in this paper can be fitted on any state-of-the-art computer running one of the standard operating systems for which R-INLA is available. CPU capacities (number of cores, CPU clock etc.) are not crucial, but we recommend that memory of at least 1 Gb is available for INLA alone, in addition to memory used by other processes including the core R software. Computing times for a single model were less than one minute (baseline models Mod1, Mod2), several minutes (Mod4), and around 6 h (Mod3, Mod5). For estimation of the full model with all data and for the cross-validation models with some of the landslide counts held out, computing times were of comparable order of magnitude in all cases. In particular, we consistently observed long computing times of several hours for Mod3 and Mod5 due to longer computations related to Laplace approximations, but the diagnostic output of R-INLA did not indicate any instabilities. Memory requirements were less than 1 Gb in all cases when using 2 cores in parallel for each job. Several models can be run in parallel on the same machine (e.g., using parallel instances of the R software), in which case available capacities (memory, number of cores) should be roughly multiplied by the number of parallel jobs.

The main computational cost of our modelling procedure stems from the Laplace approximations, performed repeatedly during an estimation run with INLA. A general rule of thumb is that computation times increase with the number of observations (5334 with our dataset), the number of latent variables (29 parameters for Mod1, 34 parameters for Mod2, and 5368 parameters in the Mod3, Mod4 and Mod5 comprising random effects), and the complexity of the additive predictor model comprising components for fixed and random effects. A major factor in the overall runtime are the computing times for linear algebra operations on sparse matrices (i.e., matrices containing mostly 0 entries), whose size is determined by the number of latent variables and the number of observations. In 3 dimensions (2 for space, 1 for time) with a model possessing *n* latent variables, we typically encounter a numerical complexity of $O(n^2)$; see Rue et al. (2009) and van Niekerk et al. (2019). In practice, complex spatial and spatio-temporal log-Gaussian Cox processes have been successfully fitted with models having a resolution *n* of several hundred thousands; see Gómez-Rubio et al. (2015), Lombardo et al. (2018a), Opitz et al. (2020), and Pimont et al. (2020) for examples. When there is a tendency towards confounding of effects, i.e., when different additive components can provide similar contributions to the predictor, computations can become less stable such that computation times increase, or estimation may even fail, which did not happen with our models.

The mapping units for which event counts are recorded may have higher resolution or may span larger areas than here, with up to several hundreds of thousands of observations (e.g., Lombardo et al., 2018a). Running INLA on such datasets, with models comprising several thousands of latent variables, then typically takes several hours, or even several days in extreme cases. With memory requirements of R-INLA easily exceeding 16 Gb in such high-dimensional models, estimation is usually carried out on machines dedicated to scientific computing. In general, Bayesian hierarchical modelling demands considerably higher computing resources than the classical frequentist approaches that do not explicitly model the spatial random effects not captured by the available covariates, but it also provides significant benefits as shown in this work. In the context of geomorphological applications, the estimation results usually have to be established only once, without any need to reestimate models for continuously updated data, such that this increased computational burden remains manageable in practice. Moreover, the library R-INLA provides several choices of less accurate approximation schemes to speed up estimation and reduce memory usage. Its recent integration of the powerful PARDISO library for numerical matrix computations further increases its potential for solving very high-dimensional problems (van Niekerk et al., 2019), e.g., where the spatial discretisation of the study area (and of time) may lead to

several hundreds of thousands of observed count variables.

## 7. Discussion

We now discuss the results of our modelling experiment. We first focus on what we consider the main advantages and the limitations of the new modelling framework, and its potential applicability to other areas (§7.1). Next, we make specific and general considerations on the results of our work for landslide hazard assessments (§7.3). This is followed by a critical analysis of the modelling approach for geomorphological and slope stability inference (§7.4). Lastly, we provide a perspective on further developments of landslide predictive modelling and zoning (§7.5).

### 7.1. A new landslide predictive modelling framework

In our work, we experimented with an innovative, Bayesian modelling framework for the spatio-temporal prediction of landslides of the slide type (Hungr et al., 2014) (§5). Results showed that the adopted framework was capable of predicting the temporal, the spatial, and the spatio-temporal distributions of known landslides that occurred in our study area in the period from before 1941 to 2014 (Fig. 1). Results also showed that considering the existing albeit not explicitly identified (i.e., “latent”), landslide dependencies in space and time improved substantially the model predictive performance (§6.6 and Fig. 11), when compared to a simpler (“traditional”) model, exemplified by our baseline model Mod1 (Fig. 4) which does not consider the spatio-temporal dependencies among landslides. The differences between Mod3 (with temporally replicated and independent but spatially structured random effects) and Mod5 (with temporal dependence, in addition to Mod3) are relatively small in comparison, but Mod5 offers additional benefits, such as the interpretation of the estimated temporal autocorrelation coefficient $\hat{\beta}^{\text{Mod5}}$, whose estimated value 0.68 was significantly different from 0. In a situation where the landslide counts were not observed for some of the slope units in some time intervals, Mod5 (and also Mod4) could improve the estimated intensities for these slope units by increasing predictive strength through the landslide counts that have been observed in these slope units during the preceding or following time units.

The proposed Log-Gaussian Cox Process (LGCP) model is a “doubly” stochastic process, with the stochasticity given by (i) a Poisson component, describing the number of landslides in each SU, and by (ii) a Gaussian component, which describes the landslide intensity (i.e., the expected count per SU) on the logarithmic scale. LGCPs allow incorporating several types of random effects in the additive log-intensity function, namely linear but also nonlinear effects of observed covariates, and nonlinear effects at a latent level that may be viewed as effects of unobserved or unavailable predictor variables; see Lombardo et al. (2018a) for examples of nonlinear effects of observed covariates. Therefore, the LGCP approach provides a convenient framework to investigate and interpret morphometric and thematic covariates' influence on slope instability in the study area (Fig. 9). It also brings to light unobserved dependencies that influence the landslide intensity function $\Lambda_i(s)$ (§5.2), and the derived landslide susceptibility estimates. We consider this as important progress with respect to the other existing, statistically-based landslide prediction modelling tools currently available in the literature (Reichenbach et al., 2018). These alternative approaches usually do not cope with latent effects, and typically do not explicitly model temporal or spatio-temporal stochastic landslide dependencies, such that variability stemming from random effects is not systematically incorporated into the model; in many landslide studies, estimation uncertainty was not studied at all. More specifically, state-of-the-art machine learning approaches such as Neural Networks or Random Forests can capture nonlinear effects and complex interactions of covariates, and could even be used with spatial and temporal
coordinates as covariates to capture predictor components not ex-
plained by the other covariates, or space/time-varying influence of
covariates. However, such models need very careful tuning and vali-
dation due to their “black-box”-nature. They work best for binary data
and are more difficult to adapt to the modelling of count data, and their
parameters and model components are less straightforward to interpret
than our LGCP models.

To construct our models, we exploited a multi-temporal landslide inventory comprised of numerous (3,379) landslides, which occurred in a significantly long period over our study area (Figs. 1, 3). The geographical (“cartographic”) and thematic (“geomorphological”) detail and the accuracy of the landslide mapping were key to inform properly our models, and to evaluate their performances (Fig. 7) and prediction skills (Fig. 11). This may be seen as a limitation of the proposed framework, which requires accurate landslide data to provide reliable intensity and susceptibility estimates. However, we maintain that in order to predict landslides in space and time such detailed and accurate information is necessary (mandatory), albeit it may be costly and time consuming to obtain it (Galli et al., 2008; Guzzetti et al., 2012). While accurate multi-temporal landslide information is available, together with relevant thematic data (e.g., Table 2), the added effort to construct and run a complex model (e.g., our model Mod5) compared to a simpler model (e.g., Mod1) is negligible, both in terms of GIS pre-processing and data preparation, and for the statistical modelling. Indeed, with the R-INLA library, we can run very complex models using a simple syntax. In our case, the difference between Mod1 and Mod5 was two additional lines of code.

Inspection of the model fitting (Fig. 7) and predictive (Fig. 11) performances further reveals that the more complex models (Mod3, Mod4, Mod5) where generally better at predicting the spatial (i.e., “where”) rather than the temporal (i.e., “when”) component. We maintain that this is due to the combined effect of (i) the inherent short-term temporal viability—and related unpredictability—of landslide phenomena at the SU scale, at least in our study area (Samia et al., 2017a, 2017b, 2018), and (ii) the number and length of the considered temporal periods and the number of landslides in each period (Fig. 3), which depend on the temporal frequency of landslides in our study area. The latter, is in turn controlled by the frequency of the landslide triggering forcing events (e.g., severe or prolonged rainfall periods, rapid snow-melt events) (Rossi et al., 2010b; Witt et al., 2010). This has hazard and geomorphological consequences, which we address below in §7.3 and §7.4.

The models accounting for the spatial landslide dependencies in-
volve smoothing residuals across adjacent SUs. This may have in-
troduced some local inconsistency or error at the boundary between the
SUs. However, as mentioned in §5.5, the models are flexible, and let the
data prevail without introducing too much noise to the intensity esti-
mates. Should one need to account for the effect of a “rigid” barrier
(e.g., a river, a major divide, a main lithologic or tectonic discontinuity)
on landslide intensity or susceptibility, two solutions are possible.
Bakka et al. (2019) have developed a model for incorporating physical
barriers, which can be fitted using R-INLA, although their method
relies on a different type of spatial effect as the one exploited in this
work. Alternatively, one can remove manually the links between ad-
jacent SUs in the adjacency matrix (Fig. 2). Since the residual
smoothing process is governed by the adjacency matrix, removing ap-
propriate links will prevent the latent effect from “propagating” from
one SU to its direct neighbors.

The computational burden of a space-time LGCP model in R-INLA
depends on, and scales well with, the size of the dataset by exploiting
random effects with sparse precision (i.e., inverse covariance) matrices.
Relatively small areas like the one used for our experiment can be in-
vestigated effectively with a standard, modern personal computer even
for spatio-temporal models (e.g., model Mod5). Larger datasets covering
large and very large areas need proportionally larger computer facil-
ities. We note here that the adoption of the SUs as the mapping unit of

reference, or of other similar terrain mapping units (Guzzetti et al.,
1999; Guzzetti, 2005; Van Westen et al., 2006), as opposed to “grid
cells” or pixels, has reduced greatly the computational burden. In
general, use of SUs facilitates the construction of complex models even
for large and very large areas covering thousands of square kilometres
(Alvioli et al., 2016). We conclude that the main limitation for per-
forming complex, space-time LGCP landslide modelling is mostly due to
the lack of accurate datasets and detailed multi-temporal landslide in-
ventories, and not to the computational requirements. This should
guide geomorphologists interested in landslide prediction modelling in
their time allocation and resource investment (Guzzetti et al., 1999).

We further note that we successfully tested the new framework (Fig. 11) for landslides predominantly of the “slide” type (Hungr et al., 2014), which are common and abundant in our study area, and in similar areas in Central Italy and elsewhere in similar physiographical settings. We acknowledge that further efforts are required to test the framework with different landslide types, since their temporal and spatial dependencies may vary. However, we do not see any geomorphological or statistical reason that should limit or hamper the applicability of the proposed framework to other landslide types.

The predictions made by all our models are valid under the general assumptions that (i) the driving forces that control the landslide processes in the study area are known and captured through the covariates used in the models, (ii) the driving forces will remain nearly the same in the foreseeable future (Fabbri et al., 2003; Guzzetti et al., 2006a), and (iii) the landslide information shown in the multi-temporal inventory is representative of the general landslide processes in the study area (Reichenbach et al., 2018). We maintain that the three assumptions hold in our study area, but the same key assumptions should be considered thoroughly when similar models are constructed in other areas.

Ultimately, we stress that all the considerations made above do not
depend on the choice of mapping unit. In other words, the same
modelling framework could be applied for grid-cells, unique condition
units, hydrological units or any other mapping units used in the lit-
erature (Reichenbach et al., 2018). Notably, if the selected mapping
unit is the grid-cell, we suggest a slight variation to the model presented
here. In fact, one could use a hierarchical structure where the count
data and covariate values are expressed at the pixel level, but the latent
effects are expressed at the slope unit level (see, Lombardo et al.,
2019a). In such a way, the resulting model could feature characteristics
of the two mapping units (Van Den Eeckhaut et al., 2009b), i.e., it
would continuously subdivide the space in a fine squared lattice and it
would contextually feature spatio-temporal effects acting over the
whole slope, while also keeping a reasonable computational burden.

7.2. *Statistical considerations*

A rigorous implementation of a model based on spatial point pattern
theory would have required that each landslide was treated as a pre-
cisely geolocated point. For practical implementation, it is usually sa-
tisfactory to know in which mapping unit a landslide occurred.
However, a few landslides in the multi-temporal inventory (Fig. 1),
mostly present in the T1 and T2 periods (Fig. 3), have a large or very
large area (for T1 and T2: $A_L \ge 2.2 \times 10^2 \text{ m}^2$; for T3, T4 and T5: $A_L \ge
5.8 \times 10^2 \text{ m}^2$), and intersect multiple SUs. Treating such large land-
slides as single “points” would have been a severe forcing from a geo-
morphological perspective. We were then faced with the choice of
conflicting either (i) with the conditional independence assumption of
points in our modelling tool, or (ii) with the empirical, geomorpholo-
gical field evidence.

The conditional independence assumption states that observed
landslide counts are independent if we know the value of the predictor
comprising the covariate information and random effects (the latter
only if they are part of the model). This assumption is common to all
well-established spatial statistical models for discrete data, irrespective
of the choice of a susceptibility model (i.e., using a Bernoulli distribution
<figure><img src="image_17.png" /><figcaption>**Fig. 13.** Panel A shows the relation between raw/observed landslide counts and the corresponding aggregated landslide extent per SU, from T1 to T6 (ordered from left to right in the same color scheme adopted throughout the manuscript). Panel B shows the relation between estimated landslide counts via Mod5 and the corresponding aggregated landslide extent per SU, from T1 to T6 (ordered from left to right in the same color scheme adopted throughout the manuscript).</figcaption></figure>
for presence-absence data) or an intensity model (i.e., using a Poisson distribution for count data), and seems difficult to abandon, especially as it is a critical assumption for using INLA.

By analogy with susceptibility studies in the existing literature (e.g., Guzzetti et al., 2006a), we chose to respect the field evidence, and we counted the presence of a landslide—or of a portion of a landslide—in each SU if the landslide area exceeded 2% of the SU area, a percentage that accounts for possible mapping errors (Carrara et al., 1991, 1995b). We acknowledge that this approach creates some dependence between events in nearby SUs and can lead to local clustering patterns of events that cannot be fully captured by models that obey the conditional independence assumption. This entails two major problems: first, the model lacks realistic small-scale behavior and may underestimate components of landslide risk at relatively small spatial scales; second, statistical inference is flawed by considering dependent observations as independent, which will cause an underestimation of uncertainty, for instance by declaring certain covariates as significant while they are not.

However, by using random effects as in our models we can substantially alleviate these problems by capturing local dependence and clustering structures that cannot be explained by geomorphological covariates alone. In other words, the latent spatial random effect can capture (part of) the dependence induced by the largest landslides affecting several SUs. While classical generalised linear models (GLMs) have only fixed effects and assume complete independence of observations, our models are based on the less restrictive assumption of conditional independence with respect to the combination of fixed and random effects. In this Bayesian framework, we can model the propagation of landslide counts over neighboring SUs, which also includes how far a single landslide extends in space (i.e., “how large” it is). One should note that non-Bayesian modelling frameworks, and estimation approaches not resorting to INLA, do exist for random effect models (i.e., for Generalised Linear Mixed Models—GLMMs, or Generalised Additive Mixed Models—GAMMs), but they often use estimation algorithms prone to providing biased predictions, and may lack flexibility in terms of available dependence forms of space-time random effects (see the general discussion in Rue et al., 2009, and Lindgren et al., 2011, for technical details). By working with intensities instead of susceptibilities, i.e., with count data instead of presence-absence data, we further reduce the loss of information in data and models when opting for larger mapping units, where the phenomenon of single landslides stretching over several units becomes less frequent. As a result, the procedure of using an intensity framework with random effects brings our models closer to the accepted definition of landslide hazard (Varnes, 1984; Guzzetti et al., 1999, 2005).

### 7.3. Hazard considerations

As mentioned in §2, the prediction of landslide hazard proposed by Varnes (1984), and later modified by Guzzetti et al. (1999) and Guzzetti et al. (2005), requires the anticipation—in probabilistic terms—of “where” (spatial component), “when” (temporal component), and “how large” or destructive (magnitude component) landslides are expected to be in an area (Guzzetti, 2005). Our new modelling framework, and specifically our more complex model Mod5, fulfils the definition, to a large extent. Model Mod5 accounts for the spatial and temporal dependencies of landslides, including latent effects not explicitly described by other covariates. We acknowledge that the space and time components of the model Mod5 are conditioned by the quality and completeness of the multi-temporal landslide information, and that the time component is further conditioned by the length of the period covered by the landslide inventory, and by the length of the sub-periods between the (irregularly spaced) “temporal slices” in the multi-temporal inventory (Fig. 1). The magnitude component of the hazard is also present in model Mod5, given by the expected number of landslides in each SU, i.e., by the landslide intensity. We note here that the number of

landslides was used as a measure of landslide event magnitude, e.g., by Keefer (1984) for earthquake-induced landslides and by Malamud et al. (2004) for landslides caused by weather and geophysical triggers. Furthermore, landslide size characteristics, including landslide area (Hovius et al., 1997; Guzzetti et al., 2002; Malamud et al., 2004), volume (Malamud et al., 2004; Brunetti et al., 2009), area-to-volume ratio (Guzzetti et al., 2009; Larsen et al., 2010), and length-to-width ratio (Taylor et al., 2018), which can all be associated to the vulnerability to landslides (Galli and Guzzetti, 2007) and hence to the landslide destructive power, are known to be empirically related to the number of landslides in an area. Furthermore, because we also considered the landslides' extent when counting slope instabilities per slope units, our models Mod3 and Mod5 are even informed by the latent fields on the persistence (a proxy for size) of landslide counts over space.

One could argue with our approach and assume that the landslide counts may not directly be linked to the landslide area, and therefore that the intensity may not address a fundamental component of the hazard definition (Guzzetti et al., 2005). For instance, this could be the case if a single landslide count would have a much larger extent than the extent of many landslide counts combined. Fig. 13 clears any doubt, at least for our study area in Central Italy (Fig. 1). In fact, the observed and predicted counts appear strongly correlated to the landslide area. The relation actually appears quite linear. This implies that the theoretical hazard map would look very close to our intensity map because they would roughly only differ through a multiplicative (positive) constant. We stress that the literature dealing with landslide magnitude unanimously models the area expressed on a logarithmic scale. Conversely, in Fig. 13 the intensity estimated via Mod5 is directly correlated to the area expressed in m², hence on the same scale as the real process, as it occurs in nature.

As a result, we conclude that the landslide intensity framework proposed by Lombardo et al. (2018a, 2019b,a) for spatial predictions, and extended here in time and space-time fulfils the requirements given by the standard definition of landslide hazard, and capable to do so within a single model. This is an noteworthy advancement over previous hazard models that considered the spatial, temporal, and the landslide area (in logarithmic scale) components separately (Guzzetti et al., 2005, 2006a), and had to further assume the independence of the three components to properly estimate landslide hazard in probabilistic terms.

Our experiment revealed that only a few SUs in our study area exhibited a continual landslide clustering or repellency trend over the entire considered period, and that most of the SUs exhibited a (randomly) varying clustering/repellency signal (Fig. 8). We take this as empirical evidence of the fact that the temporal prediction of the landslide distribution over relatively long periods, longer than about 15 years in our case, is problematic and inherently uncertain. Moreover, our model does not strictly inform us on what may happen in 50 years, and it does also not predict individual landslides, but rather their distribution in our spatio-temporal domain (from T1 to T6).

Samia et al. (2017a, 2017b, 2018), working in the same study area, identified a landslide heritage effect—which they called “landslide path dependency”—that conditions the occurrence of new landslides dependent on the location of previous landslides in the same SU, over periods of less than 15 years. Both empirical findings have consequences for hazard assessment. Neglecting the temporal dependence on landslides will underestimate hazard in SUs characterised by a clustering effect, and will overestimate hazard in SUs characterised by a repellent effect. We further note that common approaches to predict future landslide occurrences over large areas, including the definition of empirical landslide thresholds for the possible initiation of landslides from landslide and rainfall records (Aleotti, 2004; Guzzetti et al., 2007, 2008; Saito et al., 2010; Ko and Lo, 2018; Segoni et al., 2018), and the calculation of return periods from time-series of triggering events, chiefly rainfall or precipitation events (Frattini et al., 2009), assume the stationarity of the landslide processes over time. However, evidence
shows that landslide processes are not stationary in our study area, and arguably in other similar areas in Central Italy and in other similar physiographical and climatic settings. The finding poses questions on the reliability of landslide forecast and prediction models based on past landslide and rainfall records (Rossi et al., 2010b; Witt et al., 2010; Segoni et al., 2018; Guzzetti et al., 2020).

Lastly, we note that our work introduced several novel advance-
ments in landslide hazard modelling, both enabled by giving preference
to modelling intensities instead of susceptibilities. Since intensities and
susceptibilities are linked by the relationship in Eq. (15), they can easily
be transformed back and forth using the mapping units of the statistical
model, and then the additivity property of intensities can be used to
transform intensities and susceptibilities to any other mapping unit.
Moreover, we provided a robust and relatively intuitive way to classify
landslide susceptibility, obtained from the landslide intensity using Eq.
(15). The approach avoids the common problem of interpreting inter-
mediate probability estimates as a measure of intermediate or “mean”
or “moderate” susceptibility levels, which is incorrect, conceptually and
operationally, and can lead to serious problems if the susceptibility
models and associated zonings are used for practical applications
(Guzzetti et al., 2000; Galli et al., 2008; Reichenbach et al., 2018).
Finally, we proposed a new way of portraying in a single map the in-
formation provided jointly by the landslide intensity and the landslide
susceptibility estimates. We maintain that the use of a single carto-
graphic representation summarising the intensity-susceptibility in-
formation facilitates the use of the zonation for practical applications,
and the design of landslide protocols for land planning and manage-
ment (Guzzetti et al., 2000; Reichenbach et al., 2018). To the best of our
knowledge, our proposed intensity-susceptibility classification in-
troduced in §5.9 and exemplified in Fig. 12 for our best model Mod5, is
unique in the landslide hazard modelling literature.

7.4. *Geomorphological considerations*

The performance of statistically-based landslide prediction models
depends entirely on the models structure and on the data used to inform
them. If the data (the “covariates”) are accurate and meaningful, an
analysis of the model results can provide valuable insights into the
geomorphic processes that control the spatio-temporal distribution of
landslides in area, provided the modelling framework is geomorpho-
logically sound.

Concerning data, to inform our models we used covariates that are
known to represent geomorphic conditions that favour or hamper the
formation of landslides in our study area (Guzzetti et al., 2006a, 2006b;
Ardizzone et al., 2007; Galli et al., 2008), in similar geologic, physio-
graphic, and climatic settings (Carrara et al., 1991, 2003; Carro et al.,
2003; Guzzetti, 2005; Marchesini et al., 2014), and even in very dif-
ferent landscapes (Budimir et al., 2015; Goetz et al., 2015; Lombardo
et al., 2016a; Reichenbach et al., 2018). With this respect, we maintain
that our morphometric, lithological, and structural covariates (Table 2)
are sound, accurate, and meaningful landslide predictors, and that they
contribute to explain the known spatio-temporal distribution of land-
slides in our (Fig. 1) and in similar study areas.

Concerning the model structure, the LGCP framework assumes that
individual landslides in a complex landscape are the result of a point
process, in space and time. In this framework, a single landslide, i.e., a
single element of a large population of landslides, is represented by a
“point” (*s*<sub>b</sub> *t*<sub>i</sub>) defined by its spatial (*s*<sub>i</sub>) and temporal (*t*<sub>i</sub>) location, i.e.,
“where” and “when” the “point” landslide occurred in the investigated
area (Fig. 1) and period (Fig. 3). The model further assumes that the
spatio-temporal distribution of landslide points is the result of an un-
observed intensity function (λ(s, t)) that varies over space and time. It is
the stochastic variation of this intensity function that determined the
location and temporal occurrence of the landslides. The last assumption
is that the model incorporates effects carried directly by the data, i.e.,
by the covariates, and by unobserved random effects not explained

directly by the covariates. In our case, such random effects include, e.g.,
the fact that geomorphologically similar and adjacent SUs behave si-
milarly in their ability to generate landslides, and the fact that land-
slides tend to repeat in time in the same places where they occurred in
the past (Samia et al., 2018). Overall, these modelling assumptions are
reasonable, from a geomorphological perspective.

Most of the landslides in the Collazzone study area have an area smaller than Mo = 5,648 m², about 0.01 % of the size of the study area. Even the largest landslide, extending for 1.5 × 10⁶ m², covers less than 2 % of the study area. In the study area landslides are caused chiefly by severe weather events, each covering a small or very small fraction of a year, and hence an even much smaller fraction of the multi-decadal period considered by our analyses. We conclude that for the spatial and temporal evolution of the landscape that characterises the Collazzone study area, individual landslides are—or can be safely considered as—"point" events, both in space and time.

It is known that landslides in the study area are not distributed randomly in space (Fig. 1) (Guzzetti et al., 2006a, 2006b; Ardizzone et al., 2007; Galli et al., 2008), and that the size and type of the landslides are controlled by the interaction between the geometry of the slopes (chiefly terrain gradient and aspect) and the attitude of the main lithological layers (i.e., the strike and dip of sand and gravel levels, and clay laminations) (Guzzetti et al., 2006a; Marchesini et al., 2015; Santangelo et al., 2015b). Thus, adjacent “anaclinal” slopes tend to generate similar, large, deep-seated slides, or complex and compound landslides, whereas adjacent “cataclinal” slopes tend to generate similar, small shallow slides and minor rotational landslides. It is also known that landslides in the area do not occur randomly in time. As mentioned before, Samia et al. (2017a, 2017b, 2018), who worked in the same area, identified a landslide heritage effect that conditions the occurrence of new landslides dependent on the location of previous landslides over periods of less than 15 years. Our own results confirm that this heritage effect is limited in time, with only a minority of the SUs exhibiting a continual, long term clustering or repellency trend, with the vast majority of the SUs showing fluctuating dependence signals through time (Fig. 8). Indeed, this is a reasonable and expected behaviour for the medium to long-term evolution of slopes, and more generally of landscapes shaped by mass wasting processes. We conclude that the assumption that there exist “unobserved” latent effects that control and explain the spatio-temporal distribution of landslides is geomorphologically sound, and it matches and explains the existing empirical evidences.

We see two main limitations of geomorphological relevance of our
current LGCP framework. First, we do not explicitly consider the size
(e.g., area, volume) of the predicted landslides in each SU, with con-
sequences on the possibility to exploit the modelling results for erosion,
sediment and landscape evolution modelling. Second, the applicability
of the model over very long periods (centuries or millennia) remains to
be determined. For the former, new modelling frameworks will have to
be devised, and tested. For the latter, we not only lack long-term past
landslide data to train sound models, but we also lack a proper un-
derstanding of how climate may change and influence future slope in-
stabilities, in the same general area (Alvioli et al., 2018), and in other
areas (Gariano and Guzzetti, 2016). The main problem to overcome
both limitations lays in the lack of accurate, spatially distributed, multi-
temporal landslide datasets. However, the rapidly improving methods
and techniques for the automatic or semi-automatic detection and
mapping of landslides over large areas from remotely sensed data
promise to bridge this data acquisition gap (Guzzetti et al., 2012).

7.5. *Perspective*

We see a number of possible future improvements to our work, with
further specific and general modelling, hazard, and geomorphological
implications. For the specific case of the Collazzone study area, we
envision adding new covariates to the model, including covariates
describing (i) land use and land cover types, which are known to influence the size, abundance, and frequency of slope failures in the study area, and (ii) the morphometric and hydrological settings of the individual SUs, which can also influence the presence and evolution of landslides in layered sediments (Carrara et al., 1991). An additional improvement will be to add covariates describing spatio-temporal environmental variations, including, e.g., space-time changes in land use and land cover driven by different agricultural or forest practices. We also envision improving our modelling of the spatial latent effect introduced by SUs with similar or different lithological, hydrological, or structural characteristics. For the purpose, we could experiment with the incorporation of physical barriers (e.g., lithological or structural domain boundaries) using the advanced modelling proposed by Bakka et al. (2019); or we could select/deselect manually the links between adjacent SUs (Fig. 2) to consider local physical—strong (permeable) or weak (impermeable)—barriers. However, the latter solution will be tedious to implement, and may introduce unnecessary subjectivity to the modelling. Lastly, we envision using information on the size of the landslides in each SU, a relevant information not currently used by our models.

More generally, we envision testing our proposed modelling framework in other areas, considering similar and different landslide types, and similar and different spatio-temporal environmental information. This will measure the applicability and flexibility of the modelling framework in different physiographic and climatic settings. As an example, where a multi-temporal landslide inventory is available for a large area, even with a coarser temporal resolution than the multi-temporal inventory available for Collazzone, we envision using covariates describing the spatio-temporal evolution of precipitation (e.g., rainfall totals, rainfall duration, rainfall intensity, number of rainy days) to establish a complex functional link between the medium to long term evolution of the precipitation characteristics, and the occurrence (or lack of occurrence) of landslides. We expect this to improve the currently limited ability to understand landslides in the changing climate, and to provide better climate-driven landslide projections (Gariano and Guzzetti, 2016). Similarly, we foresee the possibility to test the modelling framework in areas where landslides are caused by repeated geophysical (e.g., earthquake) and severe meteorological (e.g., typhoons) triggers. Where event inventory maps can be prepared after each main triggering event, which is now feasible over large areas with the existing remote sensing and image processing technologies (Guzzetti et al., 2012), we expect this to improve our ability to model the evolution of complex landscapes dominated by mass-wasting processes under multiple geophysical and weather forcing (Burbank et al., 2003; Dadson et al., 2003; Lavé and Burbank, 2004; Gabet, 2007; Larsen et al., 2010; Booth et al., 2013).

## 8. Conclusions

We proposed a novel Bayesian modelling framework for the spatio-temporal prediction of landslides. The framework exploits a Log-Gaussian Cox Process (LGCP), which assumes that individual landslides in an area are the result of a stochastic point process driven by an unknown intensity function. We tested the modelling framework in the Collazzone area, Umbria, Central Italy, for which a detailed multi-temporal landslide inventory covering the period from before 1941 to 2014, and lithological and bedding data are available. We used this complex space-time geomorphological and geological information to prepare five statistical models of increasing complexity. Our “baseline” model (Mod1) solely relies on the information carried by morphometric and thematic properties, and does not account for the relative influence of spatial and temporal clustering of the landslide process. The second model (Mod2) is similar, but it allows for time-interval-specific regression constants. The next two models are more complex, and account for spatial (Mod3) and temporal (Mod4) latent effects. Lastly, our model Mod5 jointly accounts for latent temporal effects between

consecutive inventories and latent spatial effects between adjacent SUs. We maintain that our most complex model Mod5 fulfils the definition of landslide hazard given in the literature. Quantification of the spatial and the temporal predictive performances of the five models revealed that our most complex Mod5 performed generally better than the others model, and it also adequately captures both spatial and temporal dependencies in a single model. From these considerations, we concluded that Mod5 is our best model, and we selected it to generate a classification summarising both landslide intensity and susceptibility for our study area, providing more information than traditional susceptibility zonations for land planning and management.

Based on the results of our study using complex models for landslide counts, we draw the following general conclusions.

* The landslide intensity framework introduced by Lombardo et al. (2018a, 2019a, 2019b) for spatial predictions, and extended in this work for time and space-time domains, performs well, and it fulfils the requirements of the standard definition of landslide hazard within a single model. This is a significant advancement over previous landslide hazard modelling frameworks (Guzzetti et al., 2005, 2006a).

* For regional geomorphological evaluations or hazard assessments, individual landslides can be considered as “point” events, both in space and time, and a Log-Gaussian Cox Process (LGCP), or a similar model, is fully adequate for the statistical modelling of the spatial and temporal evolution of landslides in landscapes dominated by mass-wasting processes.

* Our study provides strong evidence that latent or “unobserved” effects (i.e., not captured by covariates) exist and they influence the spatio-temporal distribution of landslides. Considering these latent space-time landslide dependencies significantly improves the model predictive performance, compared to simpler models that neglect the space-time structure of the landslide process.

* The main limitation for complex, space-time landslide modelling resides in the availability of accurate data, and chiefly of detailed multi-temporal landslide inventories, and not in the availability of complex statistical modelling tools, which are available, or in the computational requirements, which can be relatively easily fulfilled in typical applications. This consideration should guide those interested in space-time predictive modelling of landslides in the allocation of their research time and their resource investments (Guzzetti et al., 1999, 2012).

We expect our novel approach to the spatio-temporal prediction of landslides to enhance the ability to evaluate landslide hazard and its temporal and spatial variations, to lead to better projections of future landslides, and to improve our collective understanding of the evolution of landscapes dominated by mass-wasting processes under geophysical and weather drivers. To promote reproducible analysis and replicable experiments in different geomorphological contexts, we share as Supplementary Material the dataset, the adjacency matrix and the R code used in this study.

## Declaration of Competing Interest

None.

## Acknowledgement

We thank Håvard Rue, the main developer of the R-INLA project, and Haakon Bakka for the continuous discussions and technical support through the initial stage of this research. We are grateful to the CNR IRPI geomorphology research group (http://geomorphology.irpi.cnr.it) who have provided and updated the multi-temporal landslide inventory for the Collazzone area. Without their enduring effort, this work would
not have been possible. Ultimately, we would like to genuinely thank the anonymous reviewer and Prof. Alexander Brenning as the second reviewer of the submitted manuscript. Some of the points raised during

the reviewing stage have triggered exceptionally interesting discussions within our team and contributed to improve the quality of our work even further than we originally envisioned.

## Appendix A. Variables, symbols, and acronyms

Here, we list the variables, symbols and acronyms used in the text.

<table><thead><tr><th>Variable</th><th>Units</th><th>Explanation</th></tr></thead><tbody><tr><td>β</td><td></td><td>Regression coefficient</td></tr><tr><td>ε</td><td></td><td>Innovation term in the definition of random effects</td></tr><tr><td>κ</td><td></td><td>Unconditional precision parameter of random effects</td></tr><tr><td>λ</td><td></td><td>Landslide intensity</td></tr><tr><td>τ</td><td></td><td>Conditional precision parameter of random effects</td></tr><tr><td>Λ</td><td>#</td><td>Integrated intensity, i.e., expected landslide count</td></tr><tr><td>Λ̂</td><td>#</td><td>Estimated intensity, i.e., estimated landslide count</td></tr><tr><td>A</td><td>m<sup>2</sup></td><td>Surface area of a SU</td></tr><tr><td>A<sub>L</sub></td><td>m<sup>2</sup></td><td>Surface area of a single landslide</td></tr><tr><td>A<sub>LT</sub></td><td>m<sup>2</sup></td><td>Total landslide surface area</td></tr><tr><td>N</td><td>#</td><td>Number of landslides in each SU</td></tr><tr><td>S</td><td>-</td><td>Susceptibility</td></tr><tr><td>Λ̄</td><td></td><td>Normal distribution</td></tr><tr><td>W</td><td>-</td><td>Spatial/Temporal/Spatio-temporal random effect</td></tr></tbody></table>

<table><thead><tr><th>Symbol</th><th>Explanation</th></tr></thead><tbody><tr><td>μ</td><td>Mean</td></tr><tr><td>Mo</td><td>Mode</td></tr><tr><td>σ</td><td>Standard deviation</td></tr><tr><td>sd</td><td>Standard deviation</td></tr><tr><td>s</td><td>Space</td></tr><tr><td>t</td><td>Time</td></tr><tr><td>Acronym</td><td>Explanation</td></tr><tr><td>AUC</td><td>Area Under the Curve</td></tr><tr><td>CV</td><td>Cross-Validation</td></tr><tr><td>DEM</td><td>Digital Elevation Model</td></tr><tr><td>GSD</td><td>Ground Sampling Distance</td></tr><tr><td>IR</td><td>Intensity Ratio</td></tr><tr><td>LGCP</td><td>Log-Gaussian Cox Process</td></tr><tr><td>LPS</td><td>Leica Photogrammetry Suite</td></tr><tr><td>LSE</td><td>Latent Spatial Effect</td></tr><tr><td>LTE</td><td>Latent Temporal Effect</td></tr><tr><td>LSTE</td><td>Latent Spatial and Temporal Effect</td></tr><tr><td>PC</td><td>Penalised Complexity</td></tr><tr><td>INLA</td><td>Integrated Nested Laplace Approximation</td></tr><tr><td>ROC</td><td>Receiver Operating Characteristic</td></tr><tr><td>RSP</td><td>Relative Slope Position</td></tr><tr><td>SR</td><td>Susceptibility Ratio</td></tr><tr><td>SU</td><td>Slope Unit</td></tr><tr><td>TWI</td><td>Topographic Wetness Index</td></tr><tr><td>VHR</td><td>Very High Resolution</td></tr></tbody></table>

## Appendix B. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.earscirev.2020.103318.

## References

Aleotti, P., 2004. A warning system for rainfall-induced shallow failures. Eng. Geol. 73 (3-4), 247-265.

Alvioli, M., Baum, R.L., 2016. Parallelization of the TRIGRS model for rainfall-induced landslides using the message passing interface. Environ. Model. Softw. 81, 122-135 (C).

Alvioli, M., Marchesini, I., Reichenbach, P., Rossi, M., Ardizzone, F., Fiorucci, F., Guzzetti, F., 2016. Automatic delineation of geomorphological slope units with r.slopeunits v1.0 and their optimization for landslide susceptibility modeling. Geosci. Model Dev. 9 (11), 3975-3991.

Alvioli, M., Melillo, M., Guzzetti, F., Rossi, M., Palazzi, E., von Hardenberg, J., Brunetti, M.T., Peruccacci, S., 2018. Implications of climate change on landslide hazard in Central Italy. Sci. Total Environ. 630, 1528-1543.

Amato, G., Eisank, C., Castro-Camilo, D., Lombardo, L., 2019. Accounting for covariate distributions in slope-unit-based landslide susceptibility models. a case study in the

alpine environment. Eng. Geol. 260 (In print).

Ardizzone, F., Cardinali, M., Galli, M., Guzzetti, F., Reichenbach, P., 2007. Identification and mapping of recent rainfall-induced landslides using elevation data collected by airborne lidar. Nat. Hazards Earth Syst. Sci. 7 (6), 637-650.

Ardizzone, F., Fiorucci, F., Santangelo, M., Cardinali, M., Mondini, A.C., Rossi, M., Reichenbach, P., Guzzetti, F., 2013. Very-high resolution stereoscopic satellite images for landslide mapping. In: Landslide Science and Practice. Springer, pp. 95-101.

Ayalew, L., Yamagishi, H., 2005. The application of GIS-based logistic regression for landslide susceptibility mapping in the Kakuda-Yahiko Mountains, Central Japan. Geomorphology 65 (1), 15-31.

Badoux, A., Andres, N., Techel, F., Hegg, C., 2016. Natural hazard fatalities in Switzerland from 1946 to 2015. Nat. Hazards Earth Syst. Sci. 16 (12), 2747-2768.

Bakka, H., Rue, H., Fuglstad, G.-A., Riebler, A., Bolin, D., Illian, J., Krainski, E., Simpson, D., Lindgren, F., 2018. Spatial modeling with R-INLA: a review. Wiley Interdisc. Rev. 10 (6), e1443.

Bakka, H., Vanhatalo, J., Illian, J.B., Simpson, D., Rue, H., 2019. Non-stationary Gaussian models with physical barriers. Spatial Stat. 29, 268-288.
Banerjee, S., Carlin, B.P., Gelfand, A.E., 2014. Hierarchical Modeling and Analysis for Spatial Data. Chapman and Hall/CRC.

Basu, S., Dassios, A., 2002. A Cox process with log-normal intensity. Insurance 31 (2), 297-302.

Baum, R.L., Savage, W.Z., Godt, J.W., 2008. Trigs— A Fortran Program for Transient Rainfall Infiltration and Grid-based Spatial Slope-stability Analysis, Version 2.0. Open-File Report 2008-1159. U.S. Geological Survey.

Besag, J., 1974. Spatial interaction and the statistical analysis of lattice systems. J. R. Stat. Soc. Ser. B Methodol. 36 (2), 192-225.

Beven, K., Kirkby, M.J., 1979. A physically based, variable contributing area model of basin hydrology/Un modèle à base physique de zone d'appel variable de l'hydrologie du bassin versant. Hydrol. Sci. J. 24 (1), 43-69.

Blangiardo, M., Cameletti, M., 2015. Spatial and Spatio-temporal Bayesian Models with R-INLA. John Wiley & Sons.

Böhner, J., Selige, T., 2006. Spatial prediction of soil attributes using terrain analysis and climate regionalisation. Gottinger Geograph. Abhandlungen 115, 13-28.

Booth, A.M., Roering, J.J., Rempel, A.W., 2013. Topographic signatures and a general transport law for deep-seated landslides in a landscape evolution model. J. Geophys. Res. 118 (2), 603-624.

Bout, B., Lombardo, L., van Westen, C., Jetten, V., 2018. Integration of two-phase solid fluid equations in a catchment model for flashflows, debris flows and shallow slope failures. Environ. Model. Softw. 105, 1-16.

Brabb, E.E., 1985. Innovative approaches to landslide hazard and risk mapping. In: International Landslide Symposium Proceedings, Toronto, Canada. vol. 1. pp. 17-22.

Brabb, E.E., 1989. Landslides: Extent and economic significance in the United States. In: 28th International Geological Congress Symposium on Landslides, pp. 25-50.

Brabb, E.E., 1991. The world landslide problem. Episodes 14 (1), 52-61.

Brenning, A., 2005. Spatial prediction models for landslide hazards: review, comparison and evaluation. Nat. Hazards Earth Syst. Sci. 5 (6), 853-862.

Brenning, A., Long, S., Fieguth, P., 2012. Detecting rock glacier flow structures using Gabor filters and IKONOS imagery. Remote Sens. Environ. 125, 227-237.

Brunetti, M., Guzzetti, F., Rossi, M., 2009. Probability distributions of landslide volumes. Nonlinear Process. Earth Syst. Sci. 16 (2), 179-188.

Budimir, M.E.A., Atkinson, P.M., Lewis, H.G., 2015. A systematic review of landslide probability mapping using logistic regression. Landslides 1-18.

Burbank, D.W., Blythe, A.E., Putkonen, J.K., Pratt-Sitaula, B.A., Gabet, E.J., Oskin, M., Barros, A.P., Ojha, T.P., 2003. Decoupling of erosion and precipitation in the Himalayas. Nature 426 (6967), 652-655.

Cama, M., Lombardo, L., Conoscenti, C., Agnesi, V., Rotigliano, E., 2015. Predicting storm-triggered debris flow events: application to the 2009 Ionian Peloritan disaster (Sicily, Italy). Nat. Hazards Earth Syst. Sci. 15 (8), 1785-1806.

Cama, M., Lombardo, L., Conoscenti, C., Rotigliano, E., 2017. Improving transferability strategies for debris flow susceptibility assessment: Application to the Saponara and Itala catchments (Messina, Italy). Geomorphology 288, 52-65.

Campbell, R.H., 1973. Isopleth map of landslide deposits, Point Dume Quadrangle, Los Angeles County, California; an experiment in generalizing and quantifying areal distribution of landslides. In: Technical report. U.S. Geological Survey.

Cardinali, M., Ardizzone, F., Galli, M., Guzzetti, F., Reichenbach, P., 2000. Landslides triggered by rapid snow melting: the December 1996-January 1997 event in Central Italy. In: Proceedings 1st Plinius Conference on Mediterranean Storms, pp. 439-448.

Carrara, A., 1983. Multivariate models for landslide hazard evaluation. J. Int. Assoc. Math. Geol. 15 (3), 403-426.

Carrara, A., Guzzetti, F., 2013. Geographical information systems in assessing natural hazards. Volume 5. Springer Science & Business Media.

Carrara, A., Cardinali, M., Detti, R., Guzzetti, F., Pasqui, V., Reichenbach, P., 1991. GIS techniques and statistical models in evaluating landslide hazard. Earth Surf. Process. Landf. 16 (5), 427-445.

Carrara, A., Cardinali, M., Guzzetti, F., Reichenbach, P., 1995a. GIS technology in mapping landslide hazard. In: Geographical Information Systems in Assessing Natural Hazards, Advances in Natural and Technological Hazards Research. Kluwer, Springer, Dordrecht, pp. 135-175 (ISBN 978-90-481-4561-4 978-98-015-8404-3).

Carrara, A., Cardinali, M., Guzzetti, F., Reichenbach, P., 1995b. GIS technology in mapping landslide hazard. In: Geographical information systems in assessing natural hazards. Springer, pp. 135-175.

Carrara, A., Crosta, G., Frattini, P., 2003. Geomorphological and historical data in assessing landslide hazard. Earth Surf. Process. Landforms 28 (10), 1125-1142.

Carro, M., De Amicis, M., Luzi, L., Marzorati, S., 2003. The application of predictive modeling techniques to landslides induced by earthquakes: the case study of the 26 September 1997 Umbria-Marche earthquake (Italy). Eng. Geol. 69 (1-2), 139-159.

Castro Camilo, D., Lombardo, L., Mai, P., Dou, J., Huser, R., 2017. Handling high predictor dimensionality in slope-unit-based landslide susceptibility models through LASSO-penalized Generalized Linear Model. Environ. Model. Softw. 97, 145-156.

Chakraborty, A., Goswami, D., 2016. State of the art: Three dimensional (3D) slope-stability analysis. Int. J. Geotech. Eng. 10 (5), 493-498.

Chiu, S.N., Stoyan, D., Kendall, W.S., Mecke, J., 2013. Stochastic geometry and its applications. John Wiley & Sons.

Chung, C.-J.F., Fabbri, A.G., 1999. Probabilistic Prediction Models for Landslide Hazard Mapping. J. Photogram. Remote Sens. 65 (12), 1389-1399.

Chung, C.-J.F., Fabbri, A.G., 2003. Validation of spatial prediction models for landslide hazard mapping. Nat. Hazards 30 (3), 451-472.

Cox, D.R., 1965. On the estimation of the intensity function of a stationary point process. J. R. Stat. Ser. B Methodol. 27 (2), 332-337.

Cox, D.R., Isham, V., 1980. Point processes. Volume 12. CRC Press.

Crosta, G.B., Frattini, P., 2000. Rainfall thresholds for triggering soil slips and debris flows. In: Mediterranean Storms 2000 - Proceedings Second EGS Plinius Conference, Siena, Italy, October. vol. 2000. pp. 463-487.

Crovelli, R., Coe, J., 2009. Probabilistic estimation of numbers and costs of future landslides in the San Francisco Bay region. Georisk 3 (4), 206-223.

Cruden, D., Fell, R., 1997. Quantitative risk assessment for slopes and landslides—the state of the art. Landslide Risk Assessment pp. 3-12.

Cruden, D. M. and Varnes, D. J. (1996) Landslide types and processes. In Landslides and Engineering Practice, Number 247 in Transportation Research Board, National Academy of Sciences, pp. 36-75. Washington, D.C.: Transportation Research Board, U.S. National Research Council, (1996 edition).

Dadson, S., Hovius, N., Chen, H., Dade, W.B., Hsieh, M.-L., Willett, S.D., Hu, J.-C., Horng, M.-J., Chen, M.-C., Stark, C.P., 2003. Links between erosion, runoff variability and seismicity in the Taiwan orogen. Nature 426 (6967), 648-651.

Dai, F., Lee, C., 2001. Frequency-volume relation and prediction of rainfall-induced landslides. Eng. Geol. 59 (3-4), 253-266.

Dai, F., Lee, C., 2002. Landslide characteristics and slope instability modeling using GIS, Lantau Island, Hong Kong. Geomorphology 42 (3-4), 213-228.

Das, I., Stein, A., Kerle, N., Dadhalw, V.K., 2012. Landslide susceptibility mapping along road corridors in the Indian Himalayas using Bayesian logistic regression models. Geomorphology 179, 116-125.

DeGraff, J., Canuti, P., 1988. Using isopleth mapping to evaluate landslide activity in relation to agricultural practices. Bull. Int. Assoc. Eng. Geol. 38 (1), 61-71.

Dietrich, W.E., Bellugi, D., De Asua, R.R., 2001. Validation of the shallow landslide model, SHALSTAB, for forest management. Water Sci. Appl. 2, 195-227.

Diggle, P.J., Ribeiro, P.J., 2007. Model-based Geostatistics. Springer, New York.

Diggle, P.J., Moraga, P., Rowlingson, B., Taylor, B.M., et al., 2013. Spatial and spatio-temporal log-Gaussian Cox processes: extending the geostatistical paradigm. Stat. Sci. 28 (4), 542-563.

Donnarumma, A., Revellino, P., Grelle, G., Guadagno, F.M., 2013. Slope Angle as Indicator Parameter of Landslide Susceptibility in a Geologically Complex Area. Springer, Berlin, Heidelberg, pp. 425-433.

Dowling, C.A., Santi, P.M., 2014. Debris flows and their toll on human life: a global analysis of debris-flow fatalities from 1950 to 2011. Nat. Hazards 71 (1), 203-227.

Duman, T.Y., Can, T., Emre, Ö., Keçer, M., Doğan, A., Ateş, Ş., Durmaz, S., 2005. Landslide inventory of northwestern Anatolia, Turkey. Eng. Geol. 77 (1-2), 99-114.

Erener, A., Düzgün, H.S.B., 2012. Landslide susceptibility assessment: what are the effects of mapping unit and mapping method? Environ. Earth Sci. 66 (3), 859-877.

Fabbri, A.G., Chung, C.-J.F., Cendrero, A., Remondo, J., 2003. Is prediction of future landslides possible with a GIS? Nat. Hazards 30 (3), 487-503.

Fawcett, T., 2006. An introduction to ROC analysis. Pattern Recogn. Lett. 27 (8), 861-874.

Fell, R., Corominas, J., Bonnard, C., Cascini, L., Leroi, E., Savage, W.Z., et al., 2008. Guidelines for landslide susceptibility, hazard and risk zoning for land-use planning. Eng. Geol. 102 (3-4), 99-111.

Formetta, G., Rago, V., Capparelli, G., Rigon, R., Muto, F., Versace, P., 2014. Integrated physically based system for modeling landslide susceptibility. Proc. Earth Planet. Sci. 9, 74-82.

Formetta, G., Simoni, S., Godt, J.W., Lu, N., Rigon, R., 2016. Geomorphological control on variably saturated hillside hydrology and slope instability. Water Resour. Res. 52 (6), 4590-4607.

Frattini, P., Crosta, G., Sosio, R., 2009. Approaches for defining thresholds and return periods for rainfall-triggered shallow landslides. Hydrol. Process. 23 (10), 1444-1460.

Frattini, P., Crosta, G., Carrara, A., 2010. Techniques for evaluating the performance of landslide susceptibility models. Eng. Geol. 111 (1), 62-72.

Froude, M.J., Petley, D., 2018. Global fatal landslide occurrence from 2004 to 2016. Nat. Hazards Earth Syst. Sci. 18, 2161-2181.

Gabet, E.J., 2007. A theoretical model coupling chemical weathering and physical erosion in landslide-dominated landscapes. Earth Planet. Sci. Lett. 264 (1-2), 259-265.

Galli, M., Guzzetti, F., 2007. Landslide vulnerability criteria: A case study from umbria, central italy. Environ. Manag. 40 (4), 649-665.

Galli, M., Ardizzone, F., Cardinali, M., Guzzetti, F., Reichenbach, P., 2008. Comparing landslide inventory maps. Geomorphology 94 (3-4), 268-289.

Gariano, S.L., Guzzetti, F., 2016. Landslides in a changing climate. Earth Sci. Rev. 162, 227-252.

Gebregziabher, M., DeSantis, S.M., 2010. Latent class based multiple imputation approach for missing categorical data. J. Stat. Plan. Infer. 140 (11), 3252-3262.

Glade, T., Crozier, M., Smith, P., 2000. Applying probability determination to refine landslide-triggering rainfall thresholds using an empirical “Antecedent Daily Rainfall Model”. Pure Appl. Geophys. 157 (6-8), 1059-1079.

Goetz, J.N., Guthrie, R.H., Brenning, A., 2011. Integrating physical and empirical landslides susceptibility models using generalized additive models. Geomorphology 129 (3-4), 376-386.

Goetz, J., Brenning, A., Petschko, H., Leopold, P., 2015. Evaluating machine learning and statistical prediction techniques for landslide susceptibility modeling. Comput. Geosci. 81, 1-11.

Gómez-Rubio, V., Cameletti, M., Finazzi, F., 2015. Analysis of massive marked point patterns with stochastic partial differential equations. Spatial Stat. 14, 179-196.

Gorsevski, P.V., Gessler, P.E., Boll, J., Elliot, W.J., Foltz, R.B., 2006. Spatially and temporally distributed modeling of landslide susceptibility. Geomorphology 80 (3-4), 178-198.

Grahn, T., Jaldell, H., 2017. Assessment of data availability for the development of landslide fatality curves. Landslides 14 (3), 1113-1126.

Guzzetti, F., 2005. Landslide hazard and risk assessment. Ph.D. Thesis) In: Mathematich-Naturwissenschaftlichen Fakultät der Rheinischen Friedrich-Wilhelms-Universität. University of Bonn, Bonn, Germany, pp. 33-65.

Guzzetti, F., Carrara, A., Cardinali, M., Reichenbach, P., 1999. Landslide hazard evaluation: a review of current techniques and their application in a multi-scale study,
Central Italy. Geomorphology 31 (1), 181-216.

Guzzetti, F., Cardinali, M., Reichenbach, P., Carrara, A., 2000. Comparing landslide maps: a case study in the upper Tiber River Basin, central Italy. Environ. Manag. 25 (3), 247-263.

Guzzetti, F., Malamud, B.D., Turcotte, D.L., Reichenbach, P., 2002. Power-law correlations of landslide areas in central Italy. Earth Planet. Sci. Lett. 195 (3-4), 169-183.

Guzzetti, F., Reichenbach, P., Cardinali, M., Galli, M., Ardizzone, F., 2005. Probabilistic landslide hazard assessment at the basin scale. Geomorphology 72 (1-4), 272-299.

Guzzetti, F., Galli, M., Reichenbach, P., Ardizzone, F., Cardinali, M., 2006a. Landslide hazard assessment in the Collazzone area, Umbria, Central Italy. Nat. Hazards Earth Syst. Sci. 6 (1), 115-131.

Guzzetti, F., Reichenbach, P., Ardizzone, F., Cardinali, M., Galli, M., 2006b. Estimating the quality of landslide susceptibility models. Geomorphology 81 (1-2), 166-184.

Guzzetti, F., Peruccacci, S., Rossi, M., Stark, C.P., 2007. Rainfall thresholds for the initiation of landslides in central and southern Europe. Meteorog. Atmos. Phys. 98 (3-4), 239-267.

Guzzetti, F., Peruccacci, S., Rossi, M., Stark, C.P., 2008. The rainfall intensity-duration control of shallow landslides and debris flows: an update. Landslides 5 (1), 3-17.

Guzzetti, F., Ardizzone, F., Cardinali, M., Rossi, M., Valigi, D., 2009. Landslide volumes and landslide mobilization rates in Umbria, central Italy. Earth Planet. Sci. Lett. 279 (3-4), 222-229.

Guzzetti, F., Mondini, A.C., Cardinali, M., Fiorucci, F., Santangelo, M., Chang, K.-T., 2012. Landslide inventory maps: New tools for an old problem. Earth Sci. Rev. 112 (1-2), 42-66.

Guzzetti, F., Gariano, S.L., Peruccacci, S., Brunetti, M.T., Marchesini, I., Rossi, M., Melillo, M., 2020. Geographical landslide early warning systems. Earth Sci. Rev. 200, 102973. https://doi.org/10.1016/j.earscirev.2019.102973. (In press).

Hansen, A., Franks, C., Kirk, P., Brimcombe, A., Tung, F., 1995. Application of GIS to hazard assessment, with particular reference to landslides in Hong Kong. In: Geographical Information Systems in Assessing Natural Hazards. Springer, pp. 273-298.

Heerdegen, R.G., Beran, M.A., 1982. Quantifying source areas through land surface curvature and shape. J. Hydrol. 57 (3-4), 359-373.

Hosmer, D.W., Lemeshow, S., 2000. Applied Logistic Regression, Second edition. Wiley, New York.

Hovius, N., Stark, C.P., Allen, P.A., 1997. Sediment flux from a mountain belt derived by landslide mapping. Geology 25 (3), 231-234.

Hungr, O., Evans, S., Hazzard, J., 1999. Magnitude and frequency of rock falls and rock slides along the main transportation corridors of southwestern British Columbia. Can. Geotech. J. 36 (2), 224-238.

Hungr, O., Leroueil, S., Picarelli, L., 2014. The Varnes classification of landslide types, an update. Landslides 11 (2), 167-194.

Illian, J.B., Sørbye, S.H., Rue, H., et al., 2012. A toolbox for fitting complex spatial point process models using integrated nested Laplace approximation (INLA). Ann. Appl. Stat. 6 (4), 1499-1530.

Keefer, D.K., 1984. Landslides caused by earthquakes. Geol. Soc. Am. Bull. 95 (4), 406-421.

Kirschbaum, D., Adler, R., Hong, Y., Lerner-Lam, A., 2009. Evaluation of a preliminary satellite-based landslide hazard algorithm using global landslide inventories. Nat. Hazards Earth Syst. Sci. 9 (3), 673-686.

Ko, F.W., Lo, F.L., 2018. From landslide susceptibility to landslide frequency: a territorial-wide study in Hong Kong. Eng. Geol. 242, 12-22.

Krainski, E.T., Gómez-Rubio, V., Bakka, H., Lenzi, A., Castro-Camilo, D., Simpson, D., Lindgren, F., Rue, H., 2018. Advanced Spatial Modeling with Stochastic Partial Differential Equations Using R and INLA. Chapman and Hall/CRC, New York.

Lanni, C., McDonnell, J., Hopp, L., Rigon, R., 2013. Simulated effect of soil depth and bedrock topography on near-surface hydrologic response and slope stability. Earth Surf. Process. Landf. 38 (2), 146-159.

Lari, S., Frattini, P., Crosta, G., 2014. A probabilistic approach for landslide hazard analysis. Eng. Geol. 182, 13-24.

Larsen, I.J., Montgomery, D.R., Korup, O., 2010. Landslide erosion controlled by hillslope material. Nat. Geosci. 3 (4), 247.

Lavé, J., Burbank, D.W., 2004. Denudation processes and rates in the Transverse Ranges, southern California: erosional response of a transitional landscape to external and anthropogenic forcing. J. Geophys. Res. 109 (F1), F01006.

Leoni, G., Barchiesi, F., Catallo, F., Dramis, F., Fubelli, G., Lucifora, S., Mattei, M., Pezzo, G., Puglisi, C., 2009. GIS methodology to assess landslide susceptibility: application to a river catchment of Central Italy. J. Maps 5 (1), 87-93.

Lindgren, F.K., Rue, H., Lindström, J., 2011. An explicit link between gaussian fields and gaussian markov random fields: the stochastic partial differential equation approach (with discussion). J. R. Stat. Soc. 73 (4), 423-498.

Lindgren, F., Rue, H., et al., 2015. Bayesian spatial modelling with r-inla. J. Stat. Softw. 63 (19), 1-25.

Lombardo, L., Mai, P.M., 2018. Presenting logistic regression-based landslide susceptibility results. Eng. Geol. 244, 14-24.

Lombardo, L., Cama, M., Märker, M., Rotigliano, E., 2014. A test of transferability for landslides susceptibility models under extreme climatic events: application to the Messina 2009 disaster. Nat. Hazards 74 (3), 1951-1989.

Lombardo, L., Cama, M., Conoscenti, C., Märker, M., Rotigliano, E., 2015. Binary logistic regression versus stochastic gradient boosted decision trees in assessing landslide susceptibility for multiple-occurring landslide events: application to the 2009 storm event in Messina (Sicily, southern Italy). Nat. Hazards 79 (3), 1621-1648.

Lombardo, L., Bachofer, F., Cama, M., Märker, M., Rotigliano, E., 2016a. Exploiting Maximum Entropy method and ASTER data for assessing debris flow and debris slide susceptibility for the Giampilieri catchment (north-eastern Sicily, Italy). Earth Surf. Process. Landf. 41 (12), 1776-1789.

Lombardo, L., Fubelli, G., Amato, G., Bonasera, M., 2016b. Presence-only approach to assess landslide triggering-thickness susceptibility: a test for the Mili catchment (north-eastern Sicily, Italy). Nat. Hazards 84 (1), 565-588.

Lombardo, L., Opitz, T., Huser, R., 2018a. Point process-based modeling of multiple debris flow landslides using INLA: an application to the 2009 Messina disaster. Stoch. Env. Res. Risk A. 32 (7), 2179-2198.

Lombardo, L., Saia, S., Schillaci, C., Mai, P.M., Huser, R., 2018b. Modeling soil organic carbon with Quantile Regression: dissecting predictors' effects on carbon stocks. Geoderma 318, 148-159.

Lombardo, L., Bakka, H., Tanyas, H., van Westen, C., Mai, P.M., Huser, R., 2019a. Geostatistical modeling to capture seismic-shaking patterns from earthquake-induced landslides. J. Geophys. Res. Earth Surf. 124 (7), 1958-1980.

Lombardo, L., Opitz, T., Huser, R., 2019b. 3 - Numerical recipes for landslide spatial prediction using R-INLA: a step-by-step tutorial. In: Pourghasemi, H.R., Gokceoglu, C. (Eds.), Spatial Modeling in GIS and R for Earth and Environmental Sciences. Elsevier, pp. 55-83 (ISBN 978-0-12-815226-3).

Lombardo, L., Tanyas, H., Nicu, I.C., 2020. Spatial modeling of multi-hazard threat to cultural heritage sites. Eng. Geol., 105776.

Malamud, B.D., Turcotte, D.L., Guzzetti, F., Reichenbach, P., 2004. Landslide inventories and their statistical properties. Earth Surf. Process. Landf. 29 (6), 687-711.

Marchesini, I., Ardizzone, F., Alvioli, M., Rossi, M., Guzzetti, F., 2014. Non-susceptible landslide areas in Italy and in the Mediterranean region. Nat. Hazards Earth Syst. Sci. 14 (8), 2215-2231.

Marchesini, I., Santangelo, M., Guzzetti, F., Cardinali, M., Bucci, F., 2015. Assessing the influence of morpho-structural setting on landslide abundance. Georisk 9 (4), 261-271.

McCalpin, J., 1984. Preliminary age classification of landslides for inventory mapping. In: Proceedings 21st annual Enginering Geology and Soils Engineering Symposium, pp. 5-6.

Mergili, M., Marchesini, I., Rossi, M., Guzzetti, F., Fellin, W., 2014. Spatially distributed three-dimensional slope stability modelling in a raster GIS. Geomorphology 206, 178-195.

Meusburger, K., Alewell, C., 2009. On the influence of temporal change on the validity of landslide susceptibility maps. Nat. Hazards Earth Syst. Sci. 9, 1495-1507.

Møller, J., Syversveen, A.R., Waagpetersen, R.P., 1998. Log Gaussian Cox processes. Scand. J. Stat. 25 (3), 451-482.

Montgomery, D.R., Dietrich, W.E., 1994. A physically based model for the topographic control on shallow landslides. Water Resour. Res. 30 (4), 1153-1171.

Moraga, P., 2019. Geospatial Health Data: Modeling and Visualization with R-INLA and Shiny. Chapman & Hall/CRC Biostatistics Series, Boca Raton, FL.

Moreiras, S.M., 2004. Landslide incidence zonation in the Rio Mendoza valley, Mendoza province, Argentina. Earth Surf. Process. Landf. 29 (2), 255-266.

Murillo-Garcia, F.G., Alcantara-Ayala, I., Ardizzone, F., Cardinali, M., Fiourucci, F., Guzzetti, F., 2015. Satellite stereoscopic pair images of very high resolution: a step forward for the development of landslide inventories. Landslides 12 (2), 277-291.

Nadim, F., Kjekstad, O., Peduzzi, P., Herold, C., Jaedicke, C., 2006. Global landslide and avalanche hotspots. Landslides 3 (2), 159-173.

van Niekerk, J., Bakka, H., Rue, H., Schenk, O., 2019. New Frontiers in Bayesian Modeling Using the INLA Package in R. (arXiv preprint arXiv:1907.10426).

Nilsen, T., Brabb, E.E., 1977. Slope stability studies in the San Francisco Bay region, California. Geol. Soc. Am. Rev. Eng. Geol. 3, 235-243.

Ohlmacher, G.C., 2007. Plan curvature and landslide probability in regions dominated by earth flows and earth slides. Eng. Geol. 91 (2), 117-134.

Opitz, T., 2017. Latent Gaussian modeling and INLA: A review with focus on space-time applications. J. French Stat. Soc. 158 (3), 62-85.

Opitz, T., Huser, R., Bakka, H., Rue, H., 2018. INLA goes extreme: Bayesian tail regression for the estimation of high spatio-temporal quantiles. Extremes 21, 441-462.

Opitz, T., Bonneu, F., Gabriel, E., 2020. Point-process based Bayesian modeling of space-time structures of forest fire occurrences in Mediterranean France. Spatial Stat. 100429.

Pereira, S., Zêzere, J.L., Quaresma, I., 2017. Landslide societal risk in Portugal in the period 1865-2015. In: Workshop on World Landscape Forum, pp. 491-499.

Petschko, H., Brenning, A., Bell, R., Goetz, J., Glade, T., 2014. Assessing the quality of landslide susceptibility maps—case study lower austria. Nat. Hazards Earth Syst. Sci. 14 (1), 95-118.

Pimont, F., Fargeon, H., Opitz, T., Ruffault, J., Barbero, R., Martin-StPaul, N., Rigolot, E., Riviere, M. and Dupuy, J.-L. (2020) Prediction of regional wildfire activity with a probabilistic Bayesian framework. bioRxiv https://doi.org/10.1101/2020.05.20. 105767.

Posner, A.J., Georgakakos, K.P., 2015. Normalized landslide index method for susceptibility map development in El Salvador. Nat. Hazards 79 (3), 1825-1845.

Pourghasemi, H.R., Pradhan, B., Gokceoglu, C., 2012. Application of fuzzy logic and analytical hierarchy process (AHP) to landslide susceptibility mapping at Haraz watershed, Iran. Nat. Hazards 63 (2), 965-996.

R Core Team, 2014. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria.

Raia, S., Alvioli, M., Rossi, M., Baum, R., Godt, J., Guzzetti, F., 2014. Improving predictive power of physically based rainfall-induced shallow landslide models: a probabilistic approach. Geosci. Model Dev. 7 (2), 495-514.

Reichenbach, P., Cardinali, M., De Vita, P., Guzzetti, F., 1998. Regional hydrological thresholds for landslides and floods in the Tiber River Basin (central Italy). Environ. Geol. 35 (2-3), 146-159.

Reichenbach, P., Galli, M., Cardinali, M., Guzzetti, F., Ardizzone, F., 2005. Geomorphologic mapping to assess landslide risk: concepts, methods and applications in the umbria region of central italy. In: Landslide Risk Assessment. John Wiley, Chichester, pp. 429-468.