Introduction

Rainfall-induced shallow landslides are slope instabilities of a mass of soil and/or debris, which involve the most superficial layers until around 2.0 m from ground level. Although they involve small volumes (101–105 m3) of soil, they can be densely distributed across small catchments, contributing a lot of sediments to the river network, developing into devastating debris flows, provoking significant damages to cultivations and infrastructures, and, sometimes, causing the loss of human lives (Lacasse et al. 2010). Hence, there is a pressing need for developing and implementing actions of risk mitigation and early-warning system strategies to reduce the negative effects of the occurrence of these slope instabilities at local and regional scales (Segoni et al. 2018a).

A preliminary and fundamental step is represented by the reconstruction of a reliable model of assessment of the spatial and temporal probability of occurrence of these slope instabilities (Van Westen et al. 2006; Guzzetti et al. 2020). The spatial component determines the most prone areas according to a set of predisposing factors. The temporal component defines the moment or how frequently the slope instability occurs according to particular triggering factors (Corominas et al. 2014).

The spatio-temporal prediction of rainfall-induced shallow landslides is usually performed by the reconstruction of rainfall thresholds that represent the rainfall conditions that caused the triggering of shallow landslides in a particular geological, geomorphological, and environmental setting (Guzzetti et al. 2008; Segoni et al. 2018a; Piciullo et al. 2020). This method can be easily implemented over large areas, where significant amounts of rainfall and landslide data are available, and allows to reconstruct a robust statistical model able to define the rainfall scenarios that are representative of the real triggering conditions of shallow landslides (Guzzetti et al. 2008; Bogaard and Greco 2018; Melillo et al. 2018). Soil features and geomorphological and hydrological predisposing factors are generally not considered (Chang et al. 2008; Corominas et al. 2014; Wicki et al. 2020). Furthermore, these models are highly dependent on the availability and quality of rainfall and landslide time-series of data across the analyzed study area (Nikolopoulos et al. 2014).

A second approach concerns the use of a physically based method that combines a hydrological and a slope stability model to estimate the response of the hillslope soils towards an input rainfall event (Montgomery and Dietrich 1994; Iverson 2000; Lu and Godt 2013; Montrasio et al. 2014; Canli et al. 2018; Kang et al. 2020). These models are dynamic, allowing to analyze the change in unstable zones across a test-site during a particular rainfall event (Schilirò et al. 2015; Zhuang et al. 2017). They require several soil geotechnical and hydrological input data that are not easy to be measured over large areas (Gorsevski et al. 2006). Furthermore, the hydrological boundary conditions in soil and bedrock materials are sometimes assumed as steady in wide areas (Corominas et al. 2014).

Few attempts have been performed to exploit the possibility of using data-driven techniques which integrate static predisposing and dynamic triggering factors (Wang and Sassa 2006; Lee et al. 2020; Lombardo et al. 2014, 2020). These approaches aimed to fill the intrinsic gaps of rainfall thresholds and physically based models, especially to: (i) take into account both predisposing geomorphological, geological, and hydrological features and rainfall triggering factors; (ii) be effective at local and regional scales; and (iii) represent the variation of unstable areas in time during a particular rainfall event.

Most of these approaches insert parameters of a rainfall event, collected through rain gauges or radar instruments, within the set of predictors of a data-driven algorithm to model the probability of occurrence of the triggered slope instabilities (Dai and Lee 2003; Ayalew and Yamagishi 2005; Wang and Sassa 2006; Chang et al. 2008; Chang and Chiang 2009; Capecchi et al. 2015; Lee et al. 2020). The defined algorithms were not tested further with triggering or not-triggering events different from that one used to build up the model, limiting the assessment of the predictive capability of the method and its possible application to any type of rainfall scenario. A second type of approach consists in the definition of a matrix to get information on the spatial and temporal probability of occurrence towards shallow landslides at local or regional scale (Hong and Adler 2008; Kirschbaum et al. 2012; Segoni et al. 2018b; Wei et al. 2018; Monsieurs et al. 2019; Pradhan et al. 2019). In these methods, spatial probability of shallow landslide occurrence is quantified by means of data-driven susceptibility models, while the temporal occurrence is assessed through rainfall thresholds. Instead, the definition of the rainfall thresholds implemented in these models was made with typical empirical-statistical approaches (Brunetti et al. 2010; Piciullo et al. 2020) that did not investigate in details that rainfall attributes were the most influencing for the triggering of shallow slope failures in a particular context. Furthermore, they do not analyze the effects of physical and hydrological soil characteristics on the rainfall amount required to trigger shallow landslides.

Vasu et al. (2016) and Park et al. (2019) developed a data-driven model able to assess the spatio-temporal probability of occurrence of rainfall-induced shallow failures in Gangwon Province (Southern Korea). This method combined statistically a susceptibility model and another data-driven method of the temporal probability of occurrence of these phenomena based on cumulative rainfall amount of the event, 20-days antecedent rainfalls, and physical-hydrological soil features. This method demonstrated a good effectiveness, but it neglected the soil moisture conditions for modeling the temporal occurrence of shallow landslides. Observed soil moisture conditions at the beginning of a rainfall event do not correspond often to antecedent rainfall indexes (Brocca et al. 2008), due to the complex subsurface hydrological processes which lead to landslide triggering (Lu and Godt 2013; Bordoni et al. 2015; Thomas et al. 2018; Picarelli et al. 2020). This also may cause a significant amount of false positives and/or false negatives limiting the reliability of these models (Mirus et al. 2018).

To overcome this gap, soil moisture data acquired over large areas could be inserted in a similar framework. Recent research efforts (Ray and Jacobs 2007; Ray et al. 2010; Brocca et al. 2012, 2016; Krogli et al. 2018; Dahigamuwa et al. 2018; Marino et al. 2020) demonstrated how the use of remotely sensed soil moisture could allow to extract useful information regarding the relationship between soil moisture in steep terrains and the triggering moments. Moreover, it is also fundamental testing the feasibility of implementing remote sensing observations of rainfall data in this type of models, in order to exploit the advantages and the limits of these products (Ray and Jacobs 2007; Kirschbaum et al. 2012; Brunetti et al. 2018).

For these reasons, this work aims to develop a method able to integrate the spatial and the temporal probability of occurrence of shallow landslides at a small scale. This method, built-up on a data-driven approach, has to be able to incorporate also remotely sensed data of soil moisture and rainfall, in order to assess the feasibility of these products for the estimation of spatio-temporal probability of occurrence.

Material and methods

Test-sites

The models were developed and tested in two representative sectors of the Oltrepò Pavese hilly area, corresponding to the northern termination of the Italian Apennines (Fig. 1). Oltrepò Pavese is composed by a complex set of tectonic units, structured as follows (Meisina et al. 2006; Bosino et al. 2019): (i) Internal Ligurian Units in the southern sector, characterized by cretaceous interstratified and massive limestones; (ii) External Ligurian Units in the central part of the area, composed by cretaceous and eocenic flyshes and claystones; (iii) Epiligurian Units, in stratigraphic and unconformable contact with the Ligurian Units in central and northern part of Oltrepò Pavese, formed by melanges with a block-in-matrix texture, marls, and interstratified rocks and sandstones; and (iv) Messinian and post-Messinian deposits in the northern part of the study area, composed by a succession of gypsum, marls, and poorly cemented sandstones and conglomerates. In this area, bedrock formations have typically a medium and high susceptibility to be affected by widespread slope instabilities, involving colluvial covers above bedrock or bedrock materials themselves (Bertolini et al. 2002).

Fig. 1
figure 1

Location of the test-sites (a); map of bedrock geology (b); map of the bedrock lithology (c); main geomorphological features and multi-temporal shallow landslides distribution (d)

Full size image

The representative test-sites of Oltrepò Pavese corresponded to to the catchments of Scuropasso and Versa (83 km2 wide), in north-eastern Oltrepò Pavese and to the Ardivestra catchment (47 km2 wide) in central Oltrepò Pavese.

In Scuropasso and Versa catchments, a bedrock succession composed by poorly cemented sandstones and conglomerates overlying marls and evaporitic deposits is present in the northern part. Superficial soils, derived from bedrock weathering, are mostly clayey or clayey–sandy silts, with thickness between few tens of centimeters and 2 m as measured in trenches, micro-boreholes, and landslide source areas. Hillslopes are steep, with an average slope angle between 15 and 35°. The central and southern parts are characterized by marly and calcareous flyshes, alternated with sandstones and marls. In this area, due to the different lithology of the bedrock and for the presence of slow-moving landslide deposits, colluvial soils have thicknesses ranging between 1 m and more than 4 m, and hillslopes have steepness typically ranging between 8 and 20°. About 40% of the area is covered by vineyards, with a significant presence (> 20%) of woodlands and shrub lands developed in hillslopes where agricultural activities were abandoned.

In Ardivestra catchment, alternated sandstones and marls and claystones represent the bedrock lithology and determine the development of thick soil covers, as in central and southern part of Versa-Scuropasso. Hillslopes have medium steepness, especially between 8 and 15°, and are covered mostly by sowed fields (> 50%), pastures, and shrub lands.

In both the study areas, slope altitude ranges between 60 and 600 m a.s.l.

According to Koppen’s classification, the climatic regime of these catchments is temperate/mesothermal, with a mean yearly temperature of 11–12 °C and an average yearly rainfall amount between 669 (Ardivestra catchment) and 684 mm (Versa and Scuropasso catchments), according to rainfall data collected from rain-gauge stations of ARPA Lombardy monitoring network in the period 2004-2019.

The area is significantly prone to shallow landslides (Bordoni et al. 2015, 2019). In the last 20 years, more than 2500 shallow landslides (Fig. 1) occurred as a consequence of 143 rainfall triggering events (Bordoni et al. 2019). According to Cruden and Varnes’ (1996) classification, most of the shallow landslides are classified as complex phenomena (Fig. 2), starting as roto-translational slides and evolving into flows, with ratio between length and width of 1.0–7.1 (length between 10 and 500 m, width between 10 and 70 m). These slope instabilities have sliding surfaces generally located at 1 m in depth (Bordoni et al. 2015).

Fig. 2
figure 2

Examples of typical shallow landslides in the test-sites

Full size image

Assessment of the spatio-temporal probability of occurrence of shallow landslides

Data-driven method for the assessment of the spatial probability of occurrence

Figure 3 shows the flowchart of the methodological approach adopted for the assessment of the spatio-temporal probability of occurrence of rainfall-induced shallow landslides exploiting also satellite-derived soil moisture and rainfall.

Fig. 3
figure 3

Methodological scheme for the assessment of the spatio-temporal probability of occurrence of shallow landslides

Full size image

The susceptibility (SUSC) map of a study area was calculated through a data-driven model that considered 11 predictors, according to their capacity to outline destabilizing factors of shallow landslides (Table 1). Nine morphological and hydrological predictors were extracted by a digital elevation model (DEM). These variables corresponded to: slope angle (SL), slope aspect (ASP), planform curvature (PLA), profile curvature (PRO), catchment area (CA), catchment slope (CS), topographic wetness index (TWI), topographic position index (TPI), and terrain ruggedness index (TRI). Also, land use (LU) and bedrock geology (GEO) were considered in SUSC model as categorical predictors.

Table 1 Explanatory variables used in the data-driven model for the spatial probability of occurrence of shallow landslides and their influences on shallow landslide occurrence
Full size table

Multicollinear analysis based on variance inflation factor was performed on the numerical predictors. As stated in different data-driven models for landslide susceptibility assessment (Bui et al. 2016), if this factor is equal or higher than 10, collinear variables were identified and, then, excluded from the data-driven model. Instead, all the considered variables were not collinear in the study areas and could be included in the model for the assessment of susceptibility.

A multi-temporal inventory of the shallow landslides occurred in each of the test-sites was considered a response variable. They consisted of polygons which bound the whole landslide perimeter. A semi-automatic method (Galve et al. 2015) was implemented to extract each landslide source area, which corresponded to 25% of the pixels with the highest elevation in each landslide.

A data matrix grouping predictors and response variables were prepared for each study area for the application of the susceptibility data-driven model that was based on Multivariate Adaptive Regression Splines (MARS) technique (Friedman 1991).

MARS fits complex and linear/non-linear relationships between the response variable and predictors. Despite other multivariate statistical methods, MARS divides the range of each predictor into regions and fitting a regression equation to each of them. Each region defines into the predictor variable axis a hinge function delimited by knots. MARS can also create more complex functions, defined by two or more predictors, in order to improve the comprehension of the real situation and to solve complex problems more easily than other multivariate techniques. The general structure of MARS algorithm can be written as (Eq. 1):

$$ y=\alpha +{\sum}_{n=1}^N{\beta}_n{h}_n(x) $$
(Eq. 1)

where y is the dependent variable, α is a constant, N is the number of terms, each formed by a coefficient βn and hn(x) is an n-th single basis function or a product of two or more basic functions of the independent variable x. Each basis function has the form max (0, x–k) or max (0, k–x), where x is a predictor and k is a knot.

MARS was applied in few data-driven models for landslide susceptibility assessment (Vorpahl et al. 2012; Felicísimo et al. 2013; Conoscenti et al. 2016; Pourghasemi and Rahmati 2018; Rotigliano et al. 2019; Vargas-Cuervo et al. 2019; Pourghasemi et al. 2020).

The susceptibility model of each study area was built through the following steps. First, a database, consisting of all shallow landslide pixels and the same number of randomly selected non-shallow landslide pixels, was subdivided into 2 subsets. The training set, representing 2/3 of this dataset, was used to build the MARS model, while the test set, including the remaining 1/3 of the dataset, was used to estimate the model accuracy.

In the case of a source area composed of at least two pixels, these were extracted all together and assigned either to the training set or to the test set. As regards the non-landslide pixels, the random choice of these ones was consistent with the use of a multi-temporal inventory of shallow landslides for each study area, where most of the prone sectors to landsliding have been already affected by real shallow failures.

The process of training and test selection was repeated in a 100-fold bootstrap procedure. The 100 fitted bootstrap models were used to extend the prediction to the whole area, to obtain a distribution of shallow landslide spatial probability for each pixel. The mean values of each bootstrap distribution of 100 probability values were used to obtain the final shallow landslide susceptibility map of a particular area. The confidence interval of the obtained average probability was also calculated for each pixel, to consider the uncertainty. The accuracy of the final susceptibility map was evaluated through the receiver operating characteristic (ROC) curve (Hosmer and Lemeshow 2000). In particular, the mean value of the 100 area under the ROC curve (AUROC) samples obtained from the 100-fold bootstrap procedure was calculated.

The probability values were then subdivided into four intervals through Jenks’s (1967) natural break classification: (i) low (0 < susceptibility ≤ 0.25); (ii) medium-low (0.25 < susceptibility ≤ 0.50); (iii) medium-high (0.50 < susceptibility ≤ 0.75); and (iv) high (0.75 ≤ susceptibility ≤ 1).

Moreover, a 2 × 2 a posteriori contingency table of the susceptibility distribution was reconstructed considering 0.5 a threshold to discriminate areas with a low probability of occurrence of shallow landslides and high stability (susceptibility ≤ 0.5) from the ones with a significant probability of occurrence of shallow landslides and intrinsic high instability (susceptibility > 0.5). In this way, the model could be considered a binary classifier. Its performance could be, then, assessed through true positives (TP), false positives (FP), true negatives (TN), and false negatives (FN) (Galanti et al. 2018).

Data-driven method for the assessment of the temporal probability of occurrence

Shallow landslides triggering occur in consequence of particular rainfall patterns coupled with particular soil hydrological conditions, related in particular to high degree of saturation in superficial soil levels (Lu and Godt 2013; Vasu et al. 2016; Bordoni et al. 2019). For these reasons, a data-driven model for evaluating the temporal probability of occurrence of shallow landslides in the test-sites (TEMP) was created coupling rainfall factors to soil saturation degree at the beginning of a rainfall and to the soil types, differentiated in terms of materials derived by different bedrocks.

According to their influence to hydrological conditions leading to shallow landslides (Table 2), the model considered: (i) the rainfall of one day (1-day rainfall); (ii) the cumulative rainfall amounts in different time periods (3, 5, 7, 10, 15, 30, 45, 60 days); (iii) the soil saturation degree measured at 0:00 of a particular day (SAT); and (iv) the bedrock geology (GEO). As regards the cumulative rainfall parameters, the most adopted parameters considered in previous works (Aleotti 2004; Giannecchini et al. 2012; Tien Bui et al. 2013; Vasu et al. 2016; Wei et al. 2018; Kim et al. 2020) to build models of temporal assessment of the occurrence of shallow failures were considered.

Table 2 Explanatory variables used in the data-driven model for the temporal probability of occurrence of shallow landslides and their influences on shallow landslide occurrence
Full size table

The response variable of this model was a database of the period January 2007–December 2018 indicating the days when at least one shallow landslide occurred or no phenomena were triggered, in a circular buffer of 10 km of radius centered from each rain gauge. It was chosen according to the morphology of the study area and to the density of rain gauges around the three study areas (Bordoni et al. 2019). The indications on landslides triggering were, then, transformed in a binary response variable, assigning 0 or 1 in case of landslide absence or presence, respectively. This choice was consistent with the methodologies adopted to reconstruct rainfall thresholds for shallow landslide occurrence at the catchment or regional scale, which usually consider that at least one shallow landslide occurs in a triggering day (Segoni et al. 2018a).

The predictors and the response variable were merged in a matrix for the application of the data-driven methodology. Models based on MARS algorithm were built, each considering a couple of the rainfall predictors together with the saturation degree and the bedrock geology. Each model was reconstructed through the same statistical approach implemented for estimating the spatial probability of occurrence of shallow landslides. Thus, the models were reconstructed since a database, consisting of all triggering days and the same number of days without triggering, subdivided into 2 subsets, corresponding to the training (2/3 of this dataset) and the test sets (the remaining 1/3 of the dataset).

The scenario with the highest AUROC was selected as the best model that allowed to estimate the temporal probability of occurrence within the study areas. Furthermore, a 2 × 2 a posteriori contingency table was reconstructed, considering 0.5 as a threshold to discriminate days with limited possibility of shallow landslide occurrence (temporal probability of occurrence ≤ 0.5), from the ones with a significant probability of shallow landslides triggering (temporal probability of occurrence> 0.5). In this way, the model could be considered a binary classifier, assessing its performance through the same indexes used for the assessment of the predictive capability of the spatial probability of occurrence (TP, FP, TN, FN).

Combination between spatial and temporal probability of occurrence of shallow landslides: Dynamic Landslide Probability Index

Spatio-temporal probability of occurrence is usually obtained through a joint probability of the two components ensured by independence assumption (Guzzetti et al. 2005; Vasu et al. 2016). According to this, Dynamic Landslide Probability Index (DLPI) was calculated from a joint probability of the estimated spatial and temporal probability of occurrence of shallow landslides, carried on multiplying these two probabilities (Cox and Donnelly 2011). The independence assumption was warranted since the data-driven model for the spatial occurrence of landslides did not consider rainfall factors and soil hydrological conditions, while the model for the temporal occurrence did not take into account for the geomorphological and hydrological attributes used to define the susceptibility.

Evaluation of the predictive capability of DLPI to real case studies

The reliability of DLPI was tested for real events that triggered shallow landslides in the analyzed catchments: (1) the days between 25 and 28 April 2009 in Versa and Scuropasso catchments; (2) the days between 27 February and 1 March 2016 in Ardivestra catchment. In Versa and Scuropasso catchments, shallow landslides triggered on 27 and 28 April 2009, while they occurred on 29 February 2016 in Ardivestra area. Furthermore, the evaluation of the feasibility of incorporating satellite-measured rainfall data for the correct prediction of the temporal occurrence of the shallow landslides could be performed for these events, comparing the calculation of DLPI using rainfall data acquired through traditional field rain gauges (blue circles in Fig. 1b) or through remote-sensing techniques.

The reliability of the reconstructed DLPI maps was quantified calculating the AUROC and the TP, FP, TN, and FN indexes for the days when shallow landslides triggered (27–28 April 2009 in Versa and Scuropasso, 29 February 2016 in Ardivestra). The last four indexes were calculated since a 2 × 2 a posteriori contingency table, considering 0.5 as probability threshold to discriminate stable from unstable areas, as in the reconstructed data-driven models. Moreover, the quantitative assessment of the degree of agreement among the maps of a particular day obtained using different input rainfall data was performed calculating the Cohen’s Kappa index k (Cohen 1960) between each pair of obtained DLPI maps.

Datasets

Maps based on pixel mapping unit of square cells were obtained from the models. The resolution of the pixels had to be consistent with the minimum landslide area detected in the study areas (6–7 m2). According to this, a 5-m resolution was chosen.

The geomorphological and hydrological predictors (SL, ASP, PLA, PRO, CA, CS, TWI, TPI, TRI) were retrieved from a 5-m resolution digital elevation model (DEM) realized through LiDAR data acquired in 2008 and 2010 by the Italian Ministry for Environment, Land, and Sea. LU map came from DUSAF of 2007 of Lombardy Region, at the scale 1:10000. GEO map derived from the bedrock lithological map at scale 1:10000 (Meisina et al. 2006; Bosino et al. 2019). It is worth noting that the shallow landslides of the test-sites did not affect significantly the geomorphological attributes of the affected hillslopes, due to the limited area (tens or few hundreds km2) and to the limited depth of the sliding surfaces (about 1 m from the ground level).

Daily rainfall measurements, collected in the period from January 2007 to December 2018, by a network of 21 rain gauges (blue circles in Fig. 1d), were considered predictors of the model for temporal probability of occurrence. SAT data, in correspondence of each considered rain gauge, were obtained through the Advanced SCATterometer (ASCAT) sensor onboard the Metop satellites. ASCAT is a C-band (5.255 GHz) instrument that provides soil moisture (SM) products characterized by a spatial sampling of 12.5 km and daily temporal resolution (Wagner et al. 2013). The retrieval algorithm is based on a change detection technique and allows to estimate the relative saturation in the first centimeters of soil. The quality of the SM estimates is impacted by dense vegetation, high topographic complexity, and frozen soil. In this study, the product provided within the EUMETSAT project H SAF (http://hsaf.meteoam.it/) denoted as H115 has been used.

The response variables of both the models were represented by multi-temporal shallow landslide inventories occurred in each of the test-sites from January 2007 to December 2018. The use of a multi-temporal database of the triggering areas occurred during several years allows generally to reconstruct a more reliable data-driven model and to identify better the geological and geomorphological features which could induce a bigger proneness to shallow landsliding (Bordoni et al. 2020). Three hundred sixty-five failures were detected in Ardivestra catchment (mean density of 7.8 shallow landslides per km2), while 1101 failures were found in Scuropasso-Versa catchments (mean density of 13.3 shallow landslides per km2). In the test-sites, 3% of the slopes were affected twice by shallow landslides triggering during the analyzed time span. These slope failures were detected by means of high resolution aerial photographs, satellite images, and field surveys (Bordoni et al. 2019). The triggering days of these events were obtained from reports of public administrations, newspapers, and online chronicles (Bordoni et al. 2019). There were 212 days where shallow landslides triggered. In this database, the first triggering events occurred in 6–8 February 2009.

Different DLPI maps of the two considered events were reconstructed considering different rainfall data. First, rainfall data acquired by the rain gauges of the considered network in the area of the analyzed catchments were exploited. Rainfall data were also obtained through two different data sources: (1) the Global Precipitation Measurement (GPM) Mission and (2) soil moisture (SM)-derived rainfall through the application of the SM2RAIN algorithm.

The first one is based on the application of the Integrated Multi-Satellite Retrievals for GPM (IMERG, Huffman et al. 2018) algorithm, that provides rainfall data at 0.1° × 0.1° spatial and half-hourly temporal resolutions in three modes, based on latency and accuracy: “early” (IMERG-ER, latency of 4–6 h after observation), “late” (IMERG-LR, 12–18 h), and “final” (IMERG-FR about 3 months). The early and the final runs differentiate in the calibration scheme and in the fact that IMERG-ER has a climatological rain gauge adjustment, whereas the IMERG-FR uses a month-to-month adjustment based on GPCC data. In this study, the IMERG_ER product has been used, hereinafter GPM for the sake of simplicity. The use of a product characterized by a very short latency allows to test the feasibility of such an approach for near-real time applications. This aspect is of paramount aspect and deserves to be investigated, even though the “Late” and the “Final” runs could be more reliable and accurate.

The second one is based on the application of the SM2RAIN algorithm (Brocca et al. 2014) that allows to estimate rainfall starting from SM observation. The algorithm is based on the inversion of the soil water balance equation and it is found to provide reliable rainfall data at different spatial scales (Brocca et al. 2014, 2019; Ciabatta et al. 2017, 2018, 2020; Massari et al. 2020). In this study, the rainfall data obtained through the application of SM2RAIN to ASCAT SM product (https://doi.org/10.5281/zenodo.3405563, SM2RAIN-ASCAT) have been used. More in detail, the dataset has been created by applying the SM2RAIN algorithm to the H115 SM dataset described above. The SM2RAIN parameters have been calibrated against ERA5 (Hersbach et al. 2020). In order to reduce the systematic errors in the retrieval of SM that can impact the quality of derived rainfall, a bias correction technique has been applied to SM2RAIN-derived estimates. For further details, the readers are referred to Brocca et al. (2019).

The spatial resolutions of the considered rainfall data were different, but it was coarse in all the cases. As regards the measures of the rain gauge network, rainfall measures provided in the rain gauges of each area and of the surrounding sectors were interpolated through ordinary kriging. Instead, all the pixels of GPM or SM2RAIN-ASCAT covering the study areas were considered, to create continuous maps for these products for both the study areas.

Rainfall and soil moisture maps had coarse resolution (~ 12 km) and were resampled at 5 m resolution, to make them homogeneous with the other predictors used in the models. The combination between the susceptibility map and the temporal probability of occurrence allowed to obtain dynamic maps of the evolution of DLPI during the analyzed events at 5 m resolution.

Correlation between field and satellite measures of daily rainfall was estimated in the two test-sites for 2013–2018 period. Considered rain gauges were the ones with time-series similar to those ones of GPM and SM2RAIN-ASCAT data (rain gauges 14 and 16 in Fig. 1d for Ardivestra and Versa-Scuropasso, respectively). The correlation between field and satellite measures was good, as testified by values of the correlation coefficient R higher than 0.60 (Fig. 4). SM2RAIN-ASCAT measures were more in agreement with field measures for both the study areas than GPM measures, as testified by values of R for SM2RAIN-ASCAT 0.04–0.05 bigger than the ones for GPM. According to this, rainfall satellite data were able to represent the rainy trends and most of the features of the main events in the analyzed test-sites.

Fig. 4
figure 4

Comparison between daily rainfall measured through rain gauges and satellite products in the two test-sites for 2013-2018 period: a comparison between rain gauge and GPM data in Ardivestra area; b comparison between rain gauge and SM2RAIN-ASCAT data in Ardivestra area; c comparison between rain gauge and GPM data in Versa-Scuropasso area; d comparison between rain gauge and SM2RAIN-ASCAT data in Versa-Scuropasso area

Full size image

Results

Spatial probability of occurrence of shallow landslides

According to Hosmer and Lemeshow’ (2000) classification, modeled susceptibility maps (Fig. 5) had an excellent (AUROC of 0.89 ± 0.01) and outstanding performance (AUROC of 0.90 ± 0.01) for Versa-Scuropasso and Ardivestra, respectively (Table 3). For Versa-Scuropasso area, MARS model improved significantly the estimation of the susceptibility distribution than the data-driven models reconstructed through Genetic Algorithm Method (Persichillo et al. 2017), whose AUROCs were of 0.74-0.77.

Fig. 5
figure 5

Maps of the spatial probability of occurrence (susceptibility) of shallow landslides in Versa-Scuropasso catchments (a) and Ardivestra catchment (b)

Full size image
Table 3 Distribution of the classes and predictive performance of the models of spatial probability of occurrence reconstructed for the different study areas
Full size table

The high predictive capabilities of the reconstructed models were confirmed by high values of TP (88.7–90.8%) and TN (80.6–82.0%) indexes (Table 3). The models failed the correct estimation of only a percentage of 9.2–11.3% of past shallow landslides occurred in the test-sites (Table 3). Overestimation of unstable areas (FP indexes) was higher than 18.0% in both the study areas (Table 3). Consequently, the percentage of area correctly predicted by the model in each test-site was very high, equal to 84.7 and 86.4% in Ardivestra and Versa-Scuropasso, respectively.

Versa-Scuropasso catchments had higher amount of areas with low or medium-low susceptibility class (81.1%) than Ardivestra (78.5%). High susceptibility class occupied more than 9.6% of the entire area in Ardivestra and Versa and 7.8% in Versa-Scuropasso area.

In the test-sites, most of the hillslopes affected in past by shallow landslides fell correctly in areas classified with medium-high (41.4% in Ardivestra, 55.0% in Versa-Scuropasso) or high (47.3% in Ardivestra, 35.8% in Versa-Scuropasso) susceptibility (Table 4).

Table 4 Distribution of the shallow landslides of the multi-temporal inventories in the different classes of spatial probability of occurrence
Full size table

A further validation of the susceptibility models was carried on considering another triggering event not considered in the previous multi-temporal inventory. In Ardivestra catchment, 22 source areas of shallow landslides were detected by field surveys after rainfall events occurred in November 2019, in correspondence of sectors with high susceptibility (spatial probability of occurrence higher than 0.77) not affected in past by shallow failures (Fig. 6).

Fig. 6
figure 6

Validation of the maps of spatial probability of occurrence (susceptibility) through the comparison with shallow landslides not present in the multi-temporal inventory used to build the models: shallow landslides occurred in November 2019 in Ardivestra catchment

Full size image

Temporal probability of occurrence of shallow landslides

The best model for the temporal probability of occurrence was the one that considered the cumulated rainfall amount of 3 (3-day cumulated rainfall) and 30 days (30-day cumulated rainfall), together with soil saturation degree measured by ASCAT satellite and bedrock geology features. This model was, in fact, characterized by an average AUROC of 0.98, representing an outstanding performance that was not reached by models built considering other rainfall attributes (Table 5).

Table 5 Mean and standard deviation of AUROCs, calculated for the entire multi-temporal dataset, of the MARS models reconstructed for the estimation of the temporal probability of occurrence of shallow landslides, each considering a pairs of rainfall attributestogether with soil saturation degree at the beginning of the day and bedrock geology
Full size table

Once defined the best model, TP, TN, FP, and FN indexes were calculated to verify its predictive capability. The high predictive capability of the temporal model on discriminating daily occurrence of shallow landslides triggering from days characterized by stable conditions was confirmed by high values of TP (203 of 212 days, 96.2%) and TN (3598 of 3803 days, 94.6%). The model failed the correct estimation of only 3.8% of days (9 of 212 days) when shallow landslides triggered. Overestimation of the model was 5.4% (205 days of 3803 days). A correlation could be found with saturation degree measured by ASCAT at the beginning of the day. Sixty-six percent of the days without triggering of shallow failures, but with modeled probability of occurrence higher than 0.5, was characterized by saturation degree higher than 80%. Twenty-eight percent of these days had a saturation degree between 50 and 80%, while the remaining 6% had saturation degree lower than 50%.

In the two test-sites, the rainfall attributes and the initial soil saturation degree of the triggering events were similar (Table 6), confirming how a unique data-driven model for the estimation of temporal probability of occurrence could be built for both the catchments.

Table 6 Range of the rainfall features (3-day and 30-day cumulated rainfall) and of the soil saturation degree at the beginning of an event for the shallow landslides triggering events occurred in Versa-Scuropasso and Ardivestra catchments in 2007–2018 time span
Full size table

A further validation of the model was performed by calculating the temporal probability of occurrence of shallow landslides for the period between 1 and 30 November 2019 for Ardivestra catchment. Shallow landslides were triggered in this time span in correspondence of significant rainy days, as testified by the source areas present in Fig. 6. Three-day and 30-day cumulated rainfalls were calculated from the representative rain gauge of this catchment, corresponding to the one labeled as 14 (Fig. 1d). The trend of saturation degree was measured in field through a dielectric sensor installed at 0.2 m from a ground level close to this rain gauge (Bordoni et al. 2013), at about 5 km from the hillslopes affected by the landslides. Saturation degree was estimated for this soil layer as the percentage of the ratio between the measured soil water content θ and the saturated water content θs, corresponding to 0.49 m3/m3 for this soil.

The model assessed the real moments of shallow landslides triggering for all of the 8 days when slope instabilities occurred (Fig. 7c). Other 5 days, when no shallow failures occurred, were characterized by probability of occurrence between 0.6 and 1.0, representing false positives (Fig. 7c). MARS model estimated unstable conditions (temporal probability of occurrence higher than 0.5) when soil saturation degree at the beginning of the day was higher than 80% and 3-day and 30-day cumulated amounts were over 21.8 and 158.4 mm, respectively (Figs. 7a and 7b).

Fig. 7
figure 7

Modeled temporal probability of occurrence of shallow landslides in the period between 1 and 30 November 2019: a 3-day and 30-day cumulated rainfalls; b saturation degree measured in field; c temporal probability of occurrence of shallow landslides (TEMP)

Full size image

Trend of spatio-temporal probability of occurrence of shallow landslides for real events

Only rainfall data of rain gauges and SM2RAIN-ASCAT were available for the 25–28 April 2009 event in Versa and Scuropasso area. The temporal probability of occurrence of shallow landslides increased in 27 April 2009, in consequence of an extreme daily rainfall detected by both field and satellite-derived data (Figs. 8, 9, and 10).

Fig. 8
figure 8

Temporal probability of occurrence between 25 and 28 April 2009 using rain gauge and SM2RAIN-ASCAT mesaures: a trends of 3-day cumulated rainfall; b trends of 30-day cumulated rainfall; c trends of saturation degree measured by ASCAT; d trends of temporal probability of occurrence

Full size image
Fig. 9
figure 9

Daily DLPI maps for the period between 25 and 28 April 2009 in Versa and Scuropasso catchments

Full size image
Fig. 10
figure 10

Distribution of DLPI classes for the period between 25 and 28 April 2009 in Versa and Scuropasso catchments

Full size image

In 27 April 2009, 17.7% of the study area was classified with medium-high or high probability of occurrence of shallow landslides considering rain gauge data, while 16.5% of the entire area was classified in the same classes considering SM2RAIN-ASCAT estimated rainfalls (Figs. 9, 10). 28 April 2009 was another important rainy day that caused a further increase in cumulated rainfalls, quantified similarly by both rain gauge and SM2RAIN-ASCAT in about 28–33 mm (Fig. 8a, b). Saturation degree, measured by ASCAT, identified completed saturated conditions at the beginning of the day (100%), in consequence of the rain fallen in the previous day. These conditions kept very high the temporal probability of occurrence of shallow landslides, which was equal to 1 considering both rain gauge or SM2RAIN-ASCAT data (Fig. 8d). Thus, a similar amount of the test-sites was classified by both the models in medium-high and high DLPI classes (19.9%).

During this event, shallow landslides triggered mostly in 27 April 2009, with further enlargement of the phenomena during 28 April 2009. The performance of DLPI models quantified using different rainfall inputs was the same and in the range of outstanding performance. Mean AUROCs were of 0.89, with very high TP and TN indexes (91.8 and 84.1%, respectively) and low FP and FN indexes (15.9 and 8.2%, respectively).

DLPI maps, reconstructed for the days of this event using different models with different input rainfalls, were in agreement for each modeled day, as average k index ranged between 0.97 and 1 defining an almost perfect correspondence between the maps created using different input rainfall data (Landis and Koch 1977).

For the modeled event in Ardivestra area (27 February–1 March 2016), the temporal probability of occurrence of shallow landslides increased significantly in 28 February, as confirmed by high daily rainfalls measured by rain gauge and SM2RAIN-ASCAT (Fig. 11). GPM did not detect this change, since it measured only 2.1 and 5.0 mm of 3-day and 30-day cumulated rainfalls in 28 February, respectively. DLPI maps of rain gauge and SM2RAIN-ASCAT showed a corresponding increase in medium-high and high probability during this day (Figs. 12, 13).

Fig. 11
figure 11

Temporal probability of occurrence between 27 February and 1 March 2016 using rain gauge, SM2RAIN-ASCAT, and GPM mesaures: a trends of 3-day cumulated rainfall; b trends of 30-day cumulated rainfall; c trends of saturation degree measured by ASCAT; d trends of temporal probability of occurrence

Full size image
Fig. 12
figure 12

Daily DLPI maps for the period between 27 February and 1 March 2016 in Ardivestra catchment

Full size image
Fig. 13
figure 13

Distribution of DLPI classes for the period between 27 February and 1 March 2016 in Ardivestra catchment

Full size image

In the following days, temporal probability of occurrence estimated since field rain gauge and SM2RAIN-ASCAT kept close to 1, while the one calculated since GPM increased until 0.98 since 29 February due to a higher amount of measured rainfalls (Fig. 11).

This had effects on the DLPI maps produced for these days using different rainfall inputs (Figs. 12, 13). The maps of 28 February obtained by models that considered rain gauge and SM2RAIN-ASCAT data were similar, with 17.8–21.5% of the territory classified in medium-high and high class, while the map obtained using GPM rainfalls as input was similar to that one of the previous day. Instead, all the three maps of 29 February classified the territory in the same way, with the same amount of areas in medium-high and high classes (21.5%).

In 1 March, rainfall fell again with lower intensity than the past two days, but the temporal probability of occurrence kept close to 1 considering rainfall measures obtained with different sensors. All the DLPI maps were similar to those ones of the previous day (Figs. 12, 13).

During this event, shallow landslides triggered in 28 and 29 February 2016. All the models identified correctly only 29 February, while 28 February was not identified as unstable day by model built using GPM data. For the inventory of the phenomena occurred in these days, the performance of DLPI models quantified using different rainfall inputs was the same and in the range of outstanding performance. Mean AUROCs were of 0.92, with very high TP and TN indexes (98.8 and 78.5%, respectively) and low FP and FN indexes (21.5 and 1.2%, respectively).

Moreover, DLPI maps were in agreement for 27 and 29 February and for 1 March, as average k index was of 1 defining an almost perfect correspondence between the maps created using different input rainfall data (Landis and Koch 1977).

Instead, k index was lower for the pairs GPM/rain gauge and GPM/SM2RAIN-ASCAT, in the range 0.43–0.47 (moderate agreement, Landis and Koch 1977). The agreement was only in the sectors where both the maps classified the hillslopes with low probability of occurrence. k index of the pair rain gauge/SM2RAIN-ASCAT was 0.97, indicating an almost perfect agreement.

Discussions

The developed methodology integrated different approaches for coupling the spatial and the temporal prediction of rainfall-induced shallow landslide at daily resolution, exploiting also the capability of satellite soil moisture and rainfall products.

The first model estimated the spatial probability of occurrence (susceptibility) of shallow failures. The data-driven method was effective in predicting unstable and substantially stable slopes in different settings. The performance of the data-driven model was excellent-outstanding, with more than 84% of area classified correctly by this model in both the test-sites. The resulting susceptibility maps represented well the sectors where shallow landslides were more probable during intense rainstorms. As demonstrated by the validation of the maps with shallow landslides triggering zones not included in the input multi-temporal inventory, the susceptibility maps classified, in medium-high or high probability of occurrence classes, areas where other triggering events could determine unstable conditions not reached in past. The developed data-driven technique was also very flexible in identifying the predictors which had more influence in predisposing the proneness of the hillslopes. Table 7 shows the number of models in 100-fold bootstrap step that included each predictor. The higher was the number of the models which considered a particular variable; the higher was its significance in influencing the outcomes (Lay et al. 2019). Slope angle and land use were selected by all the models in each study area, resulting the parameters that had the highest influence in controlling the proneness of the territory. Each test-site presented other significant predictors, represented in particular by aspect and topographic position index in Versa-Scuropasso area, and by bedrock geology and catchment area in Ardivestra catchment.

Table 7 Frequency of selection of the predictors in 100-fold bootstrap models used to reconstruct the suceptibility in each study area: SL, slope angle; ASP, slope aspect; PLA, planform curvature; PRO, profile curvature; CA, catchement area; CS, catchment slope; TWI, topographic wetness index; TPI, topographic position index; TRI, terrain ruggedness index; LU, land use; GEO, bedrock geology
Full size table

The importance of slope angle and land use in predicting correctly the susceptibility level of a particular area was confirmed also by the significant reduction of the predictive capability of models which did not consider these parameters in the set of the predictors. A model without slope angle was characterized by a reduction of 0.15–0.22 of AUROC, of 13.8–23.7% of TP and TN, and of 8.9–15.3% of FP and FN. A model without land use was characterized by a reduction of 0.10–0.15 of AUROC, of 7.9–12.5% of TP and TN, and of 5.2–10.1% of FP and FN. For Versa and Scuropasso, a model without slope aspect or without topographic position index resulted in a decrease of 0.05–0.10 of AUROC, of 5.2–9.8% of TP and TN, and of 3.6–9.2% of FP and FN. For Ardivestra, a model without catchment area or without bedrock geology was characterized by an AUROC lower of 0.05–0.15, TP and TN lower than 2.5–15.6%, and FP and FN lower than 3.6–9.8%.

For Versa-Scuropasso, susceptibility was medium-high and high where slope angle was higher than 18°, the slopes faced to east or south-west/west direction, and the land cover corresponded to vineyards and shrub lands. Also, steep hillslopes (slope angle ranging between 27 and 35°) covered by woods were characterized by susceptibility values close to 1. In Ardivestra, the most prone hillslopes were those ones characterized by slope angle higher than 10°, with upslope contributing area lower than 1000 m2 and covered by sowed fields or shrublands. Slopes were more unstable where bedrock geology was formed by marly and clayey deposits.

The effectiveness of the model could be explained by the ability of MARS to detect complex and, sometimes, non-linear relationships between terrain attributes and proneness to shallow landslides (Conoscenti et al. 2016).

The same advantages of MARS algorithm influence the reliable performance of the data-driven model reconstructed for the estimation of the temporal probability of occurrence of shallow landslides. The developed methodology identified the rainfall parameters (3-day and 30-day cumulated rainfalls) which allowed to predict the days with high probability of occurrence of shallow landslides triggering, together with the soil saturation degree at the beginning of the day measured through ASCAT satellite and the bedrock geology features. The predictive capability of this model was in the range of the outstanding performance, with a further validation of its goodness with recent (November 2019) data not included in the multi-temporal inventory used to reconstruct the model.

Instead, false positives were still present, as typically observed in other models of temporal probability of occurrence of slope instabilities as rainfall thresholds. False positives could indicate rainy conditions which effectively caused the triggering of shallow failures in areas far from urbanized areas which could not be detected or reported (Carrara et al. 2003; Giannecchini et al. 2012).

Soil moisture data (ASCAT) used in this model referred to the first centimeters of the soil profile, while shallow landslides occur generally at 1 m from ground level. At that depth, conditions of complete saturation, which lead to shallow landslides triggering in the study area, are reached less frequently than in superficial layers (Bordoni et al. 2015; Mirus et al. 2018). Considering saturation degree at 1 m from ground level could limit the number of false positives identified by the model.

As for the model reconstructed for the estimation of the spatial probability of occurrence, the higher was the number of the models which considered a particular predictor in the 100-fold bootstrap step; the higher was its significance in influencing the outcomes. The predictors of this model were selected for more than 95% of the simulation, with a 100% of selection for 3-day cumulated rainfall and saturation degree. Thus, all these predictors represented well the hydro-meteorological conditions which could lead to shallow landslides triggering in the study area. Three-day cumulated rainfall considered the short-term temporal effect of intense rainfalls and of their infiltration in soil. Thirty-day cumulated rainfall took into account for the effects of antecedent rainfall in increasing the soil water content and in decreasing the soil shear strength of the soil layers. Soil saturation degree and bedrock geology allowed to represent the hydrological response of the soil towards intense rainfalls according to different initial soil water status. These triggering conditions were in agreement with the results in Bordoni et al. (2019), who highlighted the combined effect of intense rainfall events with high soil water contents in promoting the triggering of shallow failures in hilly Oltrepò Pavese.

The two data-driven models were combined to obtain a joint probability index which allowed to assess dynamically the evolution of spatio-temporal probability of occurrence during an event. The reconstruction of trends of this dynamic index during past triggering events demonstrated the reliability of the procedure in identifying both the days of triggering and the areas affected by slope instabilities. In the triggering days, excellent-outstanding performance was obtained, with an estimation of real unstable areas over 90%.

Since the goodness of the results obtained for the two analyzed rainfall events, DLPI was calculated for each pixel of the study areas, along the period between January 2007 and December 2018. Rainfall measures acquired by field rain-gauges network were considered rainfall parameter for the estimation of the temporal probability of occurrence, since they had a continuous time-series of rainfall data for the entire period. DLPI of each pixel for each day of the analyzed period was compared with the triggering or not-triggering of shallow failures in correspondence of that pixel. A value of 0.5 was considered threshold to discriminate days with limited possibility of shallow landslide occurrence in each pixel (DLPI ≤ 0.5), from the ones with significant probability of shallow landslide occurrence triggering (DLPI > 0.5). In this way, the model could be considered a binary classifier and its performance could be assessed also through the same indexes used for the assessment of the predictive capability of the spatial and of the temporal probability of occurrence (TP, FP, TN, FN).

Figure 14 shows the maps of distribution of TP, FP, TN, and FN indexes in the test-sites for the 11-year time span. The effectiveness of DLPI was confirmed by high values of TP and TN along the analyzed period in both the study areas. Only 7.3% of the study areas presented TP and TN lower than 95%. For the same reasons, the other two indexes were very low, since FP and FN were higher than 5% only in 15.3 and 1.6% of the study areas, respectively. These results were in agreement with DLPI values calculated for the two analyzed events. The length of the analyzed time-span confirmed the reliability of the developed methodology in assessing the possible moments when shallow failures could occur, with an effective estimation of the most unstable areas.

Fig. 14
figure 14

Distribution of TP/FN indexes (a) and TN/FP indexes (b) for DLPI modeled in the test-sites during the period between January 2007 and December 2018

Full size image

The reliability of the predictions did not change significantly considering rainfall measured by rain gauges or through remote-sensing techniques. Even if satellite daily rainfalls tended to be underestimated respect to field measures (Brunetti et al. 2018), rainfall parameters obtained through satellite products were effective to model correctly the trend of the temporal probability of occurrence during an event, in terms of both identification of jumps in temporal probability and in real triggering moments. Despite the huge resampling of the rainfall data used to model the spatio-temporal distribution of the probability of occurrence, the implementation in the modeling scheme of different rainfall measures, collected with different spatial resolution, guarantees a similar and good predictive capability. However, the use of rainfall data with higher spatial resolution (e.g., radar measures) could overcome the limitations of other measures, allowing to reduce further the number of false positives or false negatives.

The developed methodology helped in filling several gaps of the other techniques, namely rainfall thresholds and physically based models.

The method guaranteed an estimation of both the temporal occurrence and of the most prone hillslopes where shallow failures could develop, according to meteorological pattern and soil hydrological conditions (Guzzetti et al. 2020). The proposed methodology overcome, also, the intrinsic limitations of physically based methodology, such as their feasible applications over large areas of tens of km2 due to computational burden, availability of a reliable dataset of soil geotechnical and hydrological features and correct reconstruction of geological-hydrological boundary conditions in different settings (Corominas et al. 2015). Moreover, this method represented an improvement respect to previous integrated data-driven approaches developed for similar aims, guaranteeing an index of temporal probability related not only to rainfall patterns but also to the soil hydrological conditions during the rainfall event (Mirus et al. 2018; Thomas et al. 2018).

The developed methodology allowed to obtain feasible results using rainfall and soil moisture data acquired through different devices. In particular, it allowed to identify triggering days using satellite-based rainfall products. For these reasons, the method is strongly flexible and could allow to obtain accurate landslide prediction also in areas with scarce field rainfall measures (e.g., developing countries).

However, for a comprehensive evaluation of the developed methodology and for its possible application in other contexts, it is necessary to highlight some limits that have to be taken into account.

It is necessary to compile a detailed and reliable multi-temporal inventory of past shallow landslide events, which indicates the triggering zones and, at least, the days of occurrence (Guzzetti et al. 2012; Corominas et al. 2015).

Furthermore, even if the methodology can be implemented using saturation degree and rainfall data provided by different tools, maps at a higher resolution of these parameters should be used for the estimation of the temporal probability of occurrence. This is required especially in those areas characterized by orographic effects on rainfall amount due to big variations in altitude and by strong variations in soil saturation due to the presence of very different land use covers (Dahigamuwa et al. 2018). The use of rainfall radar products and of remotely sensed soil moisture and rainfall data with higher resolution (e.g., Bauer-Marschallinger et al. 2018) or downscaling procedures (Wang et al. 2016; Dahigamuwa et al. 2018) may estimate better the rainfall and soil hydrological conditions over large and heterogeneous areas.

Soil moisture measured by satellites corresponds to the water content of the most superficial (less than 5 cm from the ground level) soil horizons (McCabe et al. 2005). Soil moisture information of this layer may not represent truly the entire soil moisture profile until bedrock (Lu and Godt 2013). The developed data-driven model for the estimation of temporal probability can describe the relations between superficial soil moisture, rainfall amounts, and antecedent events, trying to fill this gap. Instead, the model should be improved by adding, as predictor variable, a value of soil moisture/saturation degree representative of deeper soil levels, where the hydro-mechanical processes leading to shallow failures occur. This parameter could be obtained through physically based methods able to link the soil moisture patterns on the surface with the sub-surface water content (Brocca et al. 2016).

Generally, as other spatially distributed modeling approaches, false positives are identified for all the simulations carried on at local or regional scales. The implementation of this method for hazard assessment has to consider the possible overestimation of areas with high hazard level and could be integrated with other modeling approaches at site-specific scale, such as physically based analyses in areas identified as unstable by the models but where shallow landslides have never occurred (Wang and Sassa 2006).

Conclusions

A method combining a susceptibility model and a temporal component was developed and tested to provide spatial and temporal probability of occurrence of shallow landslides at a large scale. The variations of this probability dynamically were evaluated during particular rainfall events triggering shallow landslides, exploiting also the capability of satellite products of soil moisture and rainfall. The methodology overcomes several limitations of the most traditional approaches used for shallow landslides forecast, regarding in particular the use of rainfall thresholds without indications on spatial distribution of where phenomena could occur and the limited-in-space application of physically based methodology for time computation and input data availability.

The models identified the geological, geomorphological, and hydrological parameters influencing the proneness of a particular area, while it reconstructed the statistical relations between rainfall attributes and soil hydrological conditions leading to determine different temporal probability of occurrence of shallow failures. The method was tested on different past events, demonstrating a strong predictive capability of both stable/unstable areas and triggering/not triggering days. The developed methodology allowed to obtain feasible results using traditional rainfall data acquired by field rain gauges or using satellite-based rainfall products.

Considering its reliability, main advantages, and weak points, this method could allow to reconstruct hazard maps, according to rainfall scenarios and soil hydrological conditions at different return time. In this way, it could assist land use planners in making planning decisions for the community development of the study area and for implementing tools and strategies of risk reduction.

The developed methodology can be applied using rainfall data acquired since different sources and tools (e.g., rain gauges, satellite sensors, modeling approaches). Also, predicted rainfall amounts and maps could be implemented in this model, allowing to its possible implementation in a territorial early-warning system. Spatio-temporal probability of occurrence higher than 0.5, estimated by DLPI, could provide daily forecasts of shallow landslide occurrence, helping main stakeholders (e.g., civil protection agencies, local administrations) to act properly in order to reduce the negative consequences for people, infrastructures, and economic activities. For these aims, this model could be coupled with other spatially distributed approaches for the assessment of shallow landslides triggering (e.g., physically based methods at site-specific scale), to improve the prediction in uncertain areas and increase the predictive capability of the forecasting chain. In the frame of a potential use of the method as tool for early warning system strategies, it is required also to discriminate the hillslopes which were affected in past by multiple triggering events, since they could represent areas where higher warning levels and mitigation strategies should be implemented more necessarily (Piciullo et al. 2018). For the same aims, warning levels should be differentiated in terms of landslide intensity, expressed as number of predicted triggered phenomena in an area during a particular day. Moreover, the definition of different probability levels of landslide occurrence and of warning levels could take the effect of landslide path spatial and temporal dependency into account, since the possible occurrence of a shallow failure could be raised in and around already affected hillslopes for few years (Samia et al. 2020).

Its flexibility and feasibility make the method potentially applicable to different geological and geomorphological settings. In particular, thanks to the optimal integration with satellite rainfall and soil moisture products, it could be implemented in areas with scarce field data, operating as a fundamental tool for increasing the management of hazard and the resilience of those territories.