Multiseasonal probabilistic slope stability analysis of a large area of unsaturated pyroclastic soils

The analysis of slope stability over large areas is a demanding task for several reasons, such as the need for extensive datasets, the uncertainty of collected data, the difficulty of accounting for site-specific factors, and the considerable computation time required due to the size of investigated areas, which can pose major barriers, particularly in civil protection contexts where rapid analysis and forecasts are essential. However, as the identification of zones of higher failure probability is very useful for stakeholders and decision-makers, the scientific community has attempted to improve capabilities to provide physically based assessments. This study combined a transient seepage analysis of an unsaturated-saturated condition with an infinite slope stability model and probabilistic analysis through the use of a high-computing capacity parallelized platform. Both short- and long-term analyses were performed for a study area, and roles of evapotranspiration, vegetation interception, and the root increment of soil strength were considered. A model was first calibrated based on hourly rainfall data recorded over a 4-day event (December 14–17, 1999) causing destructive landslides to compare the results of model simulations to actual landslide events. Then, the calibrated model was applied for a long-term simulation where daily rainfall data recorded over a 4-year period (January 1, 2005–December 31, 2008) were considered to study the behavior of the area in response to a long period of rainfall. The calibration shows that the model can correctly identify higher failure probability within the time range of the observed landslides as well as the extents and locations of zones computed as the most prone ones. The long-term analysis allowed for the identification of a number of days (9) when the slope factor of safety was lower than 1.2 over a significant number of cells. In all of these cases, zones approaching slope instability were concentrated in specific sectors and catchments of the study area. In addition, some subbasins were found to be the most recurrently prone to possible slope instability. Interestingly, the application of the adopted methodology provided clear indications of both weekly and seasonal fluctuations of overall slope stability conditions. Limitations of the present study and future developments are finally discussed.


Introduction
The stability conditions of steep slopes covered by shallow soils are influenced by a number of factors, including soil mechanical properties, vegetation, rainfall, and evapotranspiration. Some local predisposing or instability factors may play a role such as cut slopes, mountain roads, water runoff management, local erosion, hydrogeological factors such as springs and differences in soil permeability.
A fundamental issue concerns the time scale of governing processes, which, in the case of rainfall-induced landslides, are infiltration and runoff. Both processes are controlled by soil conductivity, which also depends on soil suction (matric suction) when (and where) soil is unsaturated, i.e., the higher the suction, the lower the soil conductivity. Matric suction is defined as the difference between air pressure (u a ) and pore water pressure (u w ) and affects soil shear strength, i.e., the higher the suction, the higher the shear strength. However, the amount of water infiltrating the ground surface also depends on the hydraulic gradient, which changes within soil during a rainfall event. From this, a first important difference arises between (i) short heavy rainfall events (in the range of few hours) and (ii) less intense, prolonged rainfall (tens of hours up to several days). The former events tend to promote runoff while the latter generally result in infiltration. In addition, soil suction at the start of any rainfall event is also a key factor and could even be identified as a seasonal factor. In fact, in some steep shallow deposits of loose unsaturated soils, different seasonal average values of suction have been detected as a fundamental predisposing factor of superficial erosion (in the high suction season) and shallow landsliding (in the low suction season) as reported by Cascini et al. (2014). This is the case because (i) high suction increases soil strength (and favors slope stability) and reduces soil conductivity (and thus water infiltration) and thus promotes water runoff and superficial erosion while (ii) low suction corresponds to low initial soil strength and promotes water infiltration, which decreases soil suction and may result in slop weakening until failure.
Rainfall is also characterized by seasonal intensity-duration trends and by different sizes of affected areas over a year. Cascini et al. (2005Cascini et al. ( , 2014 showed that in southern Italy (specifically in the Campania region), (a) short heavy rainstorms occur in specific small areas at the end of the dry season when soil suction is high and mainly cause soil erosion and (b) prolonged rainfall at the end of the rainy season with low suction in the soil is that responsible for shallow landslides, which may become widespread over large areas.
Other than rainy periods, sunny days are also important for slope safety since the slopes recover in terms of soil suction and shear strength. The governing process is evapotranspiration, which depends on seasonal solar irradiations, weather conditions, and the specific features of plants and grass (including their spatial distributions, types, leaf coverage, root frameworks, and age), air and ground temperature, wind, slope exposure, agricultural practices and human activities, and slope management such as the periodic removal of grass, the cutting of plants, and irrigation.
In all of these chained processes, time scale and seasonal effects are paramount factors. Over a yearly scale, low and high suction periods alternate as well as heavy and prolonged rainfall periods, and seasonal temperatures play a role.
Taking into account all of these factors, features and processes in a rigorous (i.e., physically based) manner is still challenging given (i) the large amount of data required and (ii) the considerable computation time involved. This is particularly true when large areas and long periods are considered. However, recently developed methodological approaches have been designed to achieve an optimized characterization of slopes and soils over large areas. Furthermore, recent efforts in parallel calculation show promising enhancements of such applications Salvatici et al. 2018).
Based on these recent literature contributions, this paper proposes a methodological approach to assess the stability of a large hilly area (Cervinara territory, Campania, Italy) over different time periods and time steps: (i) 4 days (14-17 December 1999) considering hourly rain data and (ii) 4 years (January 1, 2005-December 31, 2008) considering daily rain data. The calibration and validation of the model are done for past inventoried landslides, and then, slope stability is assessed over a multiyear period. Contributions to the slope stability of vegetation and evapotranspiration were considered in the analysis.

Materials and methods
HIRESSS model for slope stability analysis Physically based distributed slope stability simulator HIRESSS (HIgh REsolution Slope Stability Simulator; Rossi et al. 2013) is a model developed to analyze shallow landslide triggering conditions over large scales at high spatial and temporal resolutions using parallel calculation Salvatici et al. 2018). The model is composed by two different models: hydrological and geotechnical ). The hydrological model is based on the hydraulic diffusivity concept and uses an analytical solution of an approximated form of the Richards equation for wet conditions (Richards 1931) with rainfall as input data while geotechnical stability analysis is based on an infinite slope stability model, accounting for the effect of strength and cohesion increases due to matric suction.
The HIRESSS model requires the following spatially distributed input data: slope gradient, effective cohesion (c′), root cohesion (c r ), friction angle (φ′), dry unit weight (γ d ), soil thickness, hydraulic conductivity (k s ), initial soil saturation (S), pore size index (l), bubbling pressure (h s , the air-entry value of the soil, i.e., matric suction where air starts to enter the largest soil pores), porosity (n), and residual water content (θ r ) and rainfall intensity data.
The HIRESSS model computes the factor of safety at each selected time step (and not only at the end of a rainfall event) and at different depths of the soil layer, and it can operate at any spatial resolution. Monte Carlo simulations are implemented in the model to manage the uncertainty of input data, and consequently, the output is the probability of failure.
Detailed information about the HIRESSS model and its applications to case studies can be found in Rossi et al. (2013), Mercogliano et al. (2013), Tofani et al. (2017), Salvatici et al. (2018), and Bicocchi et al. (2019). Salvatici et al. (2018) use a modified version of the geotechnical model that is novel in its inclusion of root cohesion (c r ) as an increase in soil effective cohesion (c′).
The new FS equation for unsaturated conditions is written as follows (Eq. 1): where φ is the friction angle, α is the slope angle, c tot the sum of c r and c′, γ d is the dry soil unit weight, y is depth, γ w is the water unit weight, h is the pressure head, h b is the bubbling pressure, and λ is the pore size index distribution. For saturated conditions, the FS equation ) becomes the following (Eq. 2): where γ sat is the saturated soil unit weight.

Soil thickness data collection and processing
In this study, we used an integrated procedure (Matano et al. 2016) to assess and map at a detailed scale (1:5000) the thickness of pyroclastic cover beds resting on mainly calcareous bedrock. The method relies on input data (thickness measurements) collected via fieldwork based on probings, dynamic penetration tests, and hand-dug pits with an average spatial distribution of approximately 160 measures per km 2 . The interpretation is supported by information gathered from geological and geomorphological field surveys, geophysical surveys, orthophoto interpretation, and digital terrain model (DTM) analysis. The appropriate use of geostatistical techniques for the interpolation of thickness data using a kriging technique is required combined with manual expert adjustments made in a geographical information system (GIS) environment. Probings were carried out with an iron rod (1.8 cm in diameter and up to 306 cm in length) inserted into pyroclastic cover beds by hand or with the aid of a hammer (weight of 0.03 kN). This allowed for depth measurements from the surface of contact between the pyroclastic cover beds and underlying calcareous bedrock or debris.
Dynamic penetration tests provided a measure of rock and soil in situ resistance to penetration up to 14 m in depth and, if rejected, an indirect measure of cover bed thickness. The test was performed by driving an iron rod (tip section of 10 cm 2 ) into the ground by repeatedly striking it with a 0.3-kN weight dropped from a distance of 20 cm.
Hand-dug pits (usually not exceeding 2 m in depth and 20 × 20 cm in width) allowed for the collection of further data on cover bed stratigraphy and basal contact with calcareous rock with an average spatial distribution of approximately 14 pits per km 2 .
Geophysical seismic surveys provided information about the geometry of layers forming the cover deposits and on the shape of contact between them and calcareous substratum. The seismic surveys recorded reflection data to up to 10 m in depth from 3 bursts (direct, reverse and intermediate) through the use of 20 geophones spaced 5 m apart and placed on the ground surface in a straight line with an average spatial distribution of approximately 3 surveys per km 2 .
On the steepest portion of the slopes, where in situ investigations were not available, the boundaries of outcropping calcareous rock (i.e., absence of cover bed deposits) were accurately mapped using slope attributes derived by DTM processing combined with air-photo interpretation and geological field observations.
Both measured data and derived information from in situ investigations and geomorphological surveys were collected into a database and managed in a GIS environment. The interpolation of measured pyroclastic cover bed thickness was carried out using the ordinary kriging technique implemented in Esri ArcGis 10.2®. Then, a heuristic procedure was applied to the interpolation results, which involved integrating and amending the kriging results based on qualitative data such as those gathered from geological field surveys along with air photos and DTM-derived geomorphological interpretations (Matano et al. 2016).
Case study of Cervinara (Southern Italy)

Geological context
The study area is located along the southern border of the Valle Caudina valley within the southern Apennines mountain chain (Campania, southern Italy) approximately 20 km NE from the town of Avellino (Fig. 1). Along the base of the slope massif is the village of Cervinara, which was affected by landslides of the December 1999 hydrogeological event. The landslides were triggered along the northern slope of the Monti di Avella massif, which is characterized by a difference in altitude of approximately 1200 m (from 290 m a.s.l. at the base to 1591 m a.s.l. at the summit) over less than 4 km of horizontal distance.
Mesozoic carbonate rocks belonging to the Alburno tectonic Unit constitute the skeleton of the mountain massif (Bonardi et al. 2009).
The relief consists of a N-NE-dipping monocline composed of steeply stratified (35-70°) Mesozoic limestone and dolomitic limestone beds. A Late Miocene pelitic-arenaceous unit (Castelvetere Fm.) locally overlies the bedrock. A widespread mantle of heterogeneous pyroclastic deposits, colluvial soils, and slope debris of several meters thick covers the carbonate bedrock. The pyroclastic deposits were produced by explosive eruptions from Vesuvius (i.e., 3.7 kyr Avellino eruption; Mastrolorenzo et al. 2006) and Phlegrean (i.e., 4.1 kyr Agnano-Monte Spina eruption; Orsi et al. 2009) volcanic vents located approximately 40 km southwest along the coastal sector of the Campania region. The pyroclastic deposits consist of pumice and ash alternated with buried soils. They vary in thickness (up to few meters) and their geotechnical properties relate to characteristics of the parent eruption and to weathering and geomorphic evolution (Fiorillo et al. 2001). At the base of the massif, coalescent alluvial-debris fans and thick colluvial-debris talus mark the transition from the slope base to the Valle Caudina plain and are filled with ancient lacustrine and alluvial deposits interbedded with Ignimbrite Campana tuff layers (ISPRA-Servizio Geologico d'Italia 2010). Slope characterization and landslide inventory For our Cervinara case study, in situ investigations were carried out from September to December 2011 by the Autorità di Bacino dei fiumi Liri-Garigliano e Volturno (AdB-LGV) project team (AdB-LGV 2013). The investigations involved 2250 probings, 213 dynamic penetration tests, 236 man-made pits, and 51 seismic surveys distributed across an area extending approximately 16.8 km 2 (Fig. 2). The steepest portions of the study area where slope angles exceed 40°were not investigated because they are characterized by widespread rock outcrops.
The mountain slopes of the Cervinara study area show a complex spatial distribution of pyroclastic cover bed thickness classes. A map (Fig. 3) with 10 depth classes was created to best represent the spatial variability of bedrock depths across the study area. The first class relates to outcropping calcareous rock; the following 6 classes of pyroclastic cover thicknesses have ranges of 0.5 m, while the last 3 classes have ranges of 1 m. A map analysis reveals that pyroclastic cover thicknesses change rapidly over short distances. Spatial variability is mainly controlled by the morphology of the calcareous bedrock and by the intensity, activity, and types of geomorphological processes occurring on slopes and the hydrographic network.
The site was divided into subzones suitable for the aggregation of the results. Subzones are defined by regrouping cells belonging to the upslope contributing area of a common flow node. Due to the average dimension of landslides occurring in the area (among 300 m 2 ) and to obtain a reasonable number of subzones, we used nodes that at least 400 pixels flow into. The subzones cover an average area of 15,000 m 2 . The subzones were assigned to three macrozones (zones 1, 2, and 3) based on the latitude coordinates of their centroid for a useful regrouping for graphs and descriptions in this text.
On December 15-16, 1999, heavy rainfall affected the slopes of the Monti di Avella massif and triggered approximately 70 shallow landslides within the weathered pyroclastic cover mantle along the calcareous mountain slopes. This landslide database was validated in the 1999 triggering sectors by interpreting air photographs and through geologic-geomorphologic field surveying (Fig. 4). The 1999 landslide source areas affected a total of approximately 22,120 m 2 , which corresponds to approximately 885 cells (each 5 × 5 m), which are used below as a magnitude descriptor for the whole 1999 landslide event. The landslides originated along the slopes of three small catchments located in the Cervinara urban area to the southeast (Fig. 4). The landslides resulted from the evolution of shallow translational slides into debris avalanches and debris flows with velocities ranging from very rapid to extremely rapid as defined in Hungr et al. (2014). Specifically, Cascini et al. (2011) reported that first (from 22:15 of December 15 to 00:40 of December 16) three debris floods occurred in the San Gennaro, Ioffredo, and Castello catchments (Fig. 5) and traveled distances of up to 2.5 km in the piedmont zone, covering an overall area of approximately 4.7 km 2 . Approximately 3 h later (at approximately 01:20 on December 16), a major debris avalanche was triggered in the Ioffredo catchment and then it split in two debris flows where the larger one propagated 1.5 km, affecting approximately 0.2 km 2 of the piedmont zone. Advanced numerical propagation modeling based on the SPH Smooth Particle Hydrodynamics technique ) confirms the in situ observation that relevant erosion processes occurred along the slope in the debriscolluvial deposits, increasing the mass of the propagating flows. At the base of the slopes, some debris flows impacted urban areas of Cervinara, harming residents and resulting in severe building destruction (Cascini et al. 2011).

Determination of soil properties
The deposits are volcanic soils derived from explosive activity from Vesuvius. The area was formerly investigated and details of the origins and later pedogenetic processes of pyroclastic soils are provided by Guadagno and Revellino (2005) while advanced hydromechanical characterization was carried out by Bilotta et al. (2005), Damiano et al. (2012), Sorbino et al. 2011, Cuomo andForesta (2015), and Cuomo and Iervolino (2016).
The soils recovered at the Cervinara test site can be distinguished by grain size distribution: (i) a well-graded ash layer with a high sandy component and a significant amount of nonplastic silt, (ii) altered ash with a large silty component and some clay, and (iii) pumices of decreasing grain size distribution in the domain of sandy gravel or gravelly sand. Correspondingly, all major physical and mechanical properties differ as well as reported in detail by Bilotta et al. (2005) and as confirmed by more recent experiments (Damiano et al. 2012;Cuomo and Foresta 2015).
Based on numerous literature contributions, we can conclude that the soil specific gravity (G s ) ranges from 2.5 and 2.6, the overall porosity of the ash deposits is very high (0.50-0.7), saturation levels (S r ) vary from 0.48 to 0.86, and dry unit weights (γ d ) range from 8.0 to 13.0 kN/m 3 .
Saturated conductivity (k s ) was measured through constant head tests performed using triaxial and oedometric apparatuses on undisturbed specimens derived from volcanic ash, altered ash, and fine pumices under different effective stresses (20-630 kPa). Tests of coarse pumices were performed on reconstructed specimens. It was observed that hydraulic conductivity is strongly influenced by the confining stress and that under the lowest confining stresses, the saturated hydraulic conductivity of ash soils and altered ash were 5.0 × 10 −5 m/s and approximately 10 −7 m/s, respectively (Damiano et al. 2012). The saturated conductivity adopted for simulations is 5.0 × 10 −5 m/s (Table 1), which is consistent with the value adopted by Cuomo and Iervolino (2016).
For unsaturated soil conditions, water retention properties were characterized based on the experimental laboratory results given in Damiano et al. (2012), which are consistent with those obtained by Sorbino et al. 2011 for the Pizzo d'Alvano site located 25 km from the Cervinara site. Experimental data obtained from drying and wetting tests show modest hysteretic effects, and thus, differences in water content between drying and wetting paths are neglected in the model. Thus, the data are well interpolated by Gardner's (1958) model, which is given by Eq. 3 where ψ is the pore water pressure divided by the unit weight of water (the pressure head), θ s is volumetric water, θ r is residual water content, and α is a model parameter for the shape of the function. The set of parameters defining the soil-water characteristic curve (SWCC) is reported in Table 1.
Evidence of mechanical behavior under both saturated and unsaturated conditions is available from triaxial tests performed under drained and undrained conditions on specimens derived from well-graded and altered ash soils. Damiano et al. (2012) and Olivares and Picarelli (2003) reported that the soils exhibited different behaviors and shear strength parameters. The wellgraded ash soil under saturated conditions had no effective cohesion and a friction angle of 38°while altered ash had a low effective cohesion value (11 kPa) and a lower friction angle (31°). Moreover, the well-graded ash exhibits brittle failure associated with unstable volumetric behavior upon shearing and is susceptible to static liquefaction while the altered ash does not. A comparison between the two mechanical behaviors shows that altered ash soils experienced a cementation process. The shear strength parameters adopted are related to well-graded ash, and cohesion was imposed at equal to 5 kPa. We assume that static liquefaction is one of the most frequent triggering mechanisms of rainfall-induced flowslides in loose granular soils (Cascini et al. 2010) and that slight cementation could have occurred also in well-graded ash along the slopes. The slope stability of multilayered shallow deposits has been previously investigated in the literature, which shows that accounting for slope layering involves using time- consuming iterative numerical methods. Such an approach can be used for single rainfall events, but it has not yet been applied for long-term analyses such as those described below. Conversely, Sorbino et al. (2007) and Cuomo and Iervolino (2016) showed that simplified approaches based on the schematization of single homogeneous deposits can be used. Therefore, our analyses were performed while assuming that the study area was covered by a single homogeneous deposit whose hydraulic and mechanical properties were accurately calibrated for this case study.

Analysis of meteorological data
Rainfall intensity data used were those acquired by the rain gauge located at S. Martino Valle Caudina (288 m a.s.l., 2 km far from Cervinara, ID 18897, UTM 469109, 4540511) with an acquisition time of 10 min. Unfortunately, while the rain gauge is positioned at an elevation of 288 m a.s.l., the targeted mountainous massif covers altitudes of 290 m a.s.l. to 1591 m a.s.l. Resulting differences in rainfall intensity were inferred from a previous comparison made by Cascini et al. (2005) for the Pizzo d'Alvano massif site where after major landslides, a rain-gauge was installed at approximately 1000 m a.s.l. in addition to an existing one located at approximately 300 m a.s.l. A comparison of rainfall intensity data collected at different altitudes shows some differences (by a factor of 1.5 to 2) in the short-term data (rainstorms of a few hours) but few discrepancies over longer periods (of a month or longer).
The dataset used for modeling covers a short-term period of 4 days corresponding to the 1999 landslide events and a long-term interval of 4 years (January 1, 2005-December 31, 2008). The resolution of rainfall data was downscaled to an hourly resolution for the 1999 events and to a daily resolution for the long-term event. This approach allowed us to decrease the computation time required while not compromising the accuracy of output results.
For the analysis of long-term daily evapotranspiration (ET), interception (I) and runoff were considered. For each day, water loss from the ground surface was evaluated. ET was defined as combination of two separate yet simultaneous processes: (i) evaporation from soil and plant surfaces and (ii) plant transpiration. ET was evaluated using the FAO-Penman-Monteith method (Allen  et al. 1998) based on an assessment of reference evapotranspiration (ET SZ ). ET SZ is defined as a ground surface covered by a hypothetical 12-cm grass reference crop under well-watered conditions with a fixed surface resistance of 70 s/m (Jensen et al. 1990) and an albedo of 0.23 (Eq. 4) where R n is net radiation calculated at the crop surface, G is the soil heat flux density at the soil surface, T is mean daily or hourly air temperatures at 1.5 to 2.5 m in height, u 2 is mean daily wind speed at 2.0 m in height, e s saturation vapor pressure from 1.5 to 2.5 m in height, e s mean actual vapor pressure from 1.5 to 2.5 m in height, Δ is the slope of the saturation vapor pressure-temperature curve, γ is the psychrometric constant, C n is the numerator constant, which changes by reference type and calculation time step and here is assumed to be 900, and C d is the denominator constant, which changes by reference type and calculation time step and here is set at 0.34 as recommended by Allen et al. (2000).
ET SZ was evaluated using meteorological data from the nearest meteorological station, i.e., at Airola (41°04′ 25.73″ N, 14°35′ 26.02″ E), located 7 km from the Cervinara site. For this station, daily information on temperature, relative humidity, wind speed, and solar radiation are available from the regional territorial office (Campania Region website). Specifically, the latent heat of vaporization (λ) varies slightly within the range of air temperatures that typically occurs in agricultural or hydrologic systems, and a constant value of λ = 2.45 MJ/kg was adopted in this study. Mean atmospheric pressure (P, kPa) is computed from site elevation (z) using a simplified formulation of the Universal Gas Law (Eq. 5) where z is the weather site elevation above mean sea level (m) equal to 284.9 m (mean elevation of Cervinara). The psychrometric constant (γ) is a function of P and has a unit of kilopascal per degree Celsius as reported in Eq. 6. The slope of the saturation vapor pressure-temperature curve (Δ) is shown in Eq. 7 where Δ has a unit of kilopascal per degree Celsius. Saturation water pressure (e s ) is a function of temperature (Eq. 8). Actual vapor pressure (e a ) is used to represent the water content (humidity) of air at the specific site (Eq. 9) where RH mean and T mean are the mean daily relative humidity and mean daily air temperature available from the Campania region website. Net radiation (R n ) is the net amount of radiant energy available at a vegetation or soil surface for evaporating water, heating the air, or heating the surface. R n for each day is provided by the Campania Region website. Soil heat flux density (G) is assumed to equal be to 0 for daily ET evaluation. Wind speed varies with height above the ground surface. For the calculation of ET SZ , wind speed 2 m above the surface is required, and therefore, wind measured at other heights must be adjusted from Eq. 10 where z w is the height of wind measured above the ground surface and u z is the wind speed measured z w m above the ground surface. Here, u z is drawn from the Campania Region website, and z w is assumed to be 2.5 m. ET is obtained from Eq. 11 where k c is the crop coefficient, which is defined as a function of (i) crop type, (ii) the crop development stage (season), and (iii) crop development length (Allen et al. 1998). Here, k c is assumed to be constant and equal to 0.6, i.e., the area is assumed to be covered by broadleaved forest.
ET was evaluated on a daily time step, and then, each daily ET value was subdivided on a 10-min long time step without rainfall, i.e., ET is assumed to be null during rainy time steps.
The capacity for the vegetated surface to intercept and store water was also evaluated. Actual canopy interception (I) at a given moment was calculated from the maximum canopy storage S max (mm) and cumulative rainfall from the beginning of the event R cum (Aston 1979), and for Eq. 12, where 푆 max is the maximum canopy storage (mm), k is a parameter related to canopy openness, and R cum is total precipitation (mm) (see Eq. 13). Here, R cum is defined as the total precipitation for a day. Both k and 푆 max are a function of the leaf area index (LAI), which is fixed at equal to 3 (Eq. 14), where co is canopy openness, which is constant at 0.9 (for Broadleaved Forest). Contrary to what was done for ET, here interception (I) was evaluated daily and then subdivided into 10min time steps when rainfall did occur.  the average I cum /R cum for 1 year were measured as 17% and 6%, respectively. Cumulated values R cum , ET cum , and I cum were also evaluated over 4 days. Over the multiyear period, ET cum /R cum is always lower than 23% (maximum value observed in the summer of 2005) with an average value of approximately Fig. 6 a Rainfall data used for short-term analysis and b rainfall, evapotranspiration, and interception used as input data for long-term analysis Fig. 7 Hourly counts of cells with a factor of safety of FS < 1.2 over four probability classes (20 < P < 40, 40 < P < 60, 60 < P < 80, and 80 < P < 100 with values expressed in %) for December [14][15][16][17]1999 1.66% while I cum /R cum is lower than 1% with an average value of 0.23% (Fig. 6b). Then, ratios ET cum /R cum and I cum /R cum evaluated over 4 days were clustered into a monthly set of observations and their trends were estimated for each month of 2005-2008. In December, the ET cum /R cum ratio trend is fairly constant over the 4 years and is always lower than 2% while I cum /R cum is approximately 1%. Thus, it is reasonable to assume that for December 14-17, 1999 (short-term analysis), evapotranspiration (ET) and interception (I) can be neglected. Thus, the short-term analysis was carried out considering only rainfall data.

Multitemporal scale modeling
Back analysis of December 1999 landslides An HIRESSS simulation of the study area was carried out for the 1999 event for December [14][15][16][17]1999. The simulation was based on the geotechnical parameters described in "Determination of soil properties" and reported in Table 1. The rainfall data have an hourly resolution (Fig. 6a) and the HIRESSS simulation used a Digital Terrain Model of 5 m resolution, meaning that for each pixel of 5 × 5 m, we have hourly values for the probability of failure. The model's high spatiotemporal accuracy allows us to Table 2 Daily count of cells with a probability (P) of slope instability (FS < 1.2) of the five classes (S1-S5) for December [14][15][16][17]1999 Class P (%) December 14, 1999December 15, 1999December 16, 1999December 17, 1999  monitor the gradual worsening (and later recovery) of slope stability over time during and after rainy days.
The potential of such an approach is here presented with reference to hourly computed probability (P) to produce a factor of safety (FS) for each cell of < 1.2. This value was considered to account for some uncertainty in the model (such as geotechnical parameters and DTM resolution) and to conduct a more robust analysis over several year period. Then, all computed P values were aggregated into five classes (P < 20, 20 < P < 40, 40 < P < 60, 60 < P < 80, and 80 < P < 100 with values expressed in %), and the corresponding cells were counted for December 14-17, 1999 on an hourly (Fig. 7) or daily basis (Table 2).
In Fig. 7, the period with the highest values of the hourly probability of failure runs from December 15 02:00 h to December 16 08:00 h (Fig. 7), which is quite consistent with field evidence and landslide occurrence periods. This result especially pertains to classes 60 < P < 80 and 80 < P < 100. As such a detailed temporal output is not as straightforward for longer time periods (as shown below in our multiyear analysis), we aggregate the results daily. This was done for the back analysis presented in Table 2, which shows that: (a) the two highest classes of P (60 < P < 80 and 80 < P < 100) have the highest number of cells for December 15, and (b) the number of cells that better reproduces field evidence for the 1999 landslides is that simulated for class S5 in Table 2.  The spatially distributed results of the simulation are reported in Fig. 8. The results have been daily aggregated, meaning that for each 5 × 5 pixel and map, the maximum value of failure probability reached over each 24 h is displayed. The simulation results show stable slope conditions for December 14 and December 17 while on December 15 and 16, a worsening of slope stability occurs in several parts of the study area with an increase in pixels with failure probabilities of over 60% from 2720 for December 14, 1999 to 7458 for December 15, 1999 and 6055 for December 16, 1999 ( Table 2).
The HIRESSS results were spatially validated based on real landslides occurring during the event. The 1999 rain event triggered 70 landslides (Fig. 5). Validation was carried out for the basins defined in "Slope characterization and landslide inventory." Landslides occurring during the 1999 event affected 50 of the 1743 basins in the study area. Of course, basins "affected" by landslide(s) can be defined differently. In particular, one should fix the minimum value of failure probability (P fail ) and the minimum number (N fail ) of unstable cells with probability values of higher than P fail . In Fig. 9, different P fail -N fail pairs are adopted, including 80-1, 80-3, 50-15, and 45-15. The first pair (80-1) is very conservative once applied to landslide forecasting, as it requires just one cell to be unstable to declare a subbasin as affected. However, this local approach can be very sensitive to dataset quality and related uncertainties. The second pair (80-3) requires the same threshold for P (> 80%) but at least 3 DTM cells to be affected (an area equal to or greater than 75 m 2 ). This latter value is quite similar to the average size of landslide source areas observed in the field after the 1999 events. In the other two cases, cells of larger than 10 or 15 are considered with probabilities of 50% or 45%, respectively. In these cases, we are looking for subbasins where more than a single landslide occurred. Under this rationale, the so-defined affected basins are computed for December 15 and are shown in Fig. 9.
A visual comparison of subbasins computed as "landslideaffected" to the landslide source area inventory shows a general agreement of results with field evidence. However, the validation results are reported in Table 3 in the form of a contingency matrix where for different combinations of FP and numbers of pixels, true positive (TP), false negative (FP), true negative (TN), and false positive (FP) basins are defined. In particular, "FP/pixels number" represents the combination of the thresholds of failure probability and the number of pixels to consider a basin as unstable (i.e., at least x pixels with a failure probability of higher than y); "forecast landslides" is the number of landslides whose perimeter is at least partially overlaid with basins defined as unstable in the model simulation; "TP" is the true positive percentage (correct alarms), i.e., unstable areas correctly localized by the simulation (the number of unstable basins identified by the model with at least one landslide relative to actual total unstable basins); "FN" is the false negative percentage (missing alarms), i.e., unstable areas not localized by the model (the number of stable basins identified by the model with at least one actual landslide relative to the total number of actual unstable basins); "TN" is the true negative percentage (correct nonalarms), i.e., areas correctly defined as stable by the model (basins identified as stable by the model and without landslides relative to the actual number of stable basins); and "FP" is the false positive percentage (incorrect alarms), i.e., areas incorrectly defined as stable by the model (unstable basins as identified by the model without landslides relative to the total actual number of stable basins).
Our validation from the contingency matrix shows satisfactory results overall. The percentage of forecasted landslides (i.e., the number of landslides whose perimeter is at least partially overlaid with basins found unstable by the model) for all FP/pixel value combinations is higher than 84% (up to 92.8% for combination 80/ 1), the TP percentage is higher than 60%, and the FP never exceeds 10%. However, of the combinations of P fail -N fail pairs, 80/3 seems to be the best solution in terms of validation results and dimensions of detectable landslides since 3 DTM cells (75 m 2 ) cover an area quite similar to the average size of landslide source areas observed in the field after the 1999 events.

Long-term analysis of slope stability conditions
The model was then applied to the long time series of 2005-2008. Particularly, cells with a probability (P > 80%) of slope instability (FS < 1.2) were counted on a daily basis as shown in Fig. 10.
The geotechnical input, spatial resolution of the terrain model (5 × 5 m), and time step of the rainfall data (1 h) used for the longterm analysis are the same those of the back analysis simulation ("Back analysis of December 1999 landslides"). The model computes and provides the failure probability of each pixel for every hourly time step of the 4-year simulation, and therefore for every pixel of the study area, 35,064 changes in failure probability are available.
For the long-term analysis, the daily max failure probability of pixels was considered. For each year, the number of pixels with a max failure probability of higher than 80% recorded daily in the area was calculated (Fig. 10).
To analyze the model's spatial behavior for the long-term multiseasonal simulation, an in-depth analysis was carried out on selected days of 2006, as this year registers a relevant number of days with numerous pixels with a failure probability exceeding 80%. To select significant days, a threshold of 3000 pixels of higher than 80% was set. This value was selected because it is more or less equivalent to the number of unstable pixels (80 < P < 100) found for the most unstable day (December 15, 1999) of the back analysis of the December 1999 event (Table 2). Following this approach, 9 days of 2006 were selected (Fig. 10). For each day of the spatial distribution, depending on the subbasin considered, the number of unstable pixels (failure probability > 80%) was analyzed (Fig. 11). Subbasins are grouped into three zones (zones 1, 2, and 3) from the left-to-right sections of the study area. In general, subbasins in zones 2 and 3 are more prone to the formation of shallow landslides with the highest number of unstable pixels. Furthermore, over the whole time period, the same subbasins of zones 2 and 3 show the highest number of unstable pixels. These areas can be defined as the most susceptible given local geological and geomorphological conditions. By contrast, for the same days, we find differences between subbasins with and without unstable pixels and between the same subbasins depending on the number of unstable pixels for different days. In this case, variations are mainly related to total daily precipitation and its daytime distribution.
This general trend is also shown in Fig. 12, where following the approach used for the back analysis of the 1999 event (Fig. 9), the subbasins are classified as unstable when at least three pixels with a failure probability of more than 80% are present in a basin. The results of this classification for three selected days (June 02, 2006, June 06, 2006, and August 13, 2006 of those analyzed in Fig. 11 are displayed in Fig. 12. Here, the subbasins shown in green are stable (i.e., subbasins with less than three pixels with a failure probability of higher than 80%) while those shown in yellow, orange, and red have increasing numbers of pixels with 80% probability.

Concluding remarks
The physical assessment of landslide susceptibility over large areas is critical to efficiently manage risks of shallow landslides evolving into catastrophic flows. However, field data uncertainties, the large areas to be investigated and analyzed, seasonal trends, and site-specific factors are difficult to account for in a GIS-assisted numerical model. This paper contributes to this area in terms of methodologies proposed through a relevant case study focused on southern Italy.
The proposed methodology involves acquiring detailed yet well spatially distributed data on soil deposit thickness, which heavily affects rainfall infiltration and unsaturated seepage flow. Soil thickness data processing was conducted by kriging interpolation for large areas considering topographical and geomorphological features. A probabilistic Monte Carlo simulation was incorporated into the developed model to manage input data uncertainties, and consequently, the output is the failure probability of each cell, which is then postprocessed for subcatchments of a fixed maximum size. As additional model feature, a highly parallelized algorithm was used to manage hundreds or thousands of iterations over some millions of cells of the digital terrain model (DTM). Five slope failure probability (FP) classes were defined and probability values were computed over short-(several days) and longterm periods (several years) in a multiseasonal analysis. A target factor of safety (FS) of less than 1.2 was applied to provide a robust susceptibility descriptor for long time period (here of 4 years), which is a case when uncertainty is amplified by both spatial and temporal uncertainty. In such a context, the evaluation of slope-atmosphere interactions is an important model feature for slopes subject to both inward (e.g., rainfall) and outward flux (e.g., evapotranspiration) with some rainfall intercepted by trees and canopies. The methodology also involves a postprocessing of numerical simulations. For model calibration, contingency matrixes (based on TP, FP, TN, and FN indicators) were constructed over time (4 days) and with reference to the five probability classes listed above.
Of the three acceptable solutions of the model, all satisfactorily reproducing landslide in situ evidence, one was individuated through a multicriteria procedure and selected as the one to emulate (i) a number of landslide susceptible cells similar to (or slightly higher than) those observed in the field; (ii) the maximum number of simulated susceptible cells within the observed landslide time range; (iii) the maximum percentage of "forecasted landslides," meaning a total or partial Fig. 12 Distribution of unstable subbasins (with at least 3 pixels with a failure probability FP of > 80%) for 3 days (02.06, 06.06, and 13.08) of 2006 overlapping of the simulated susceptible areas with those observed; and (iv) a small-moderate value of overall model error (i.e., a false positive (FP)). Over a long-term analysis, the model results were again classified based on the probability classes and the results were examined in relation to the number of cells with a daily FS of < 1.2.
The methodology was applied to the studied case history to produce a detailed 5-m cell DTM, accurate soil deposit thickness map, and advanced geotechnical dataset. The achieved results were used to back analyze events in 1999 in time and by spatial location. The use of the calibrated model to analyze a period of 4 years revealed a few days when slope stability conditions worsened due to heavy rainfall, which did not lead to slope failure in the examined area but did so in other zones of the region.
The proposed approach is mainly limited in that the stratified conditions of soil deposits (which play a role as indicated in the literature) are not taken explicitly into account but are addressed by means of reasonable approximations (Cuomo and Iervolino 2016) as done in the present case study. However, to overcome this limitation but not only, a probabilistic approach was introduced and the target FS was fixed to 1.2. The model does not state exact moments and locations for slope instabilities, which could be an unrealistic task for large areas and long time series for thin stratified deposits. The model provides failure probabilities for each pixel and rainfall change that can be used by stakeholders and decision-makers to identify areas with high landslide probability with a useful spatial and temporal accuracy.

Funding
Open access funding provided by Università degli Studi di Firenze within the CRUI-CARE Agreement.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/ 4.0/.