1. Introduction
Large scale occurrence of slope failures can be verified either when the equilibrium between forces acting on the slopes and their own resistances are imminently unbalanced through the action of short term forcings (dynamic triggers) or progressively through the cumulative and repetitive action of inner or external forces [
1,
2].
Periodic fluctuations of the hydraulic pore pressures can drive slope deformations toward failures. These processes are responsible for medium- to long-term variation in strength and can be considered preparatory factors for slope instability. Weathering, damaging or progressive erosion are all considered slow and effective preparatory actions [
3] for landsliding on natural slopes. At the same time, changes in the water content and soil saturation degree in the shallow and unsaturated portion of the slopes can impose suction stresses and control the position of capillary fringes, which provide transitory stabilization [
4].
Changes in hydraulic and physical conditions can instead statically predispose the action of short-term forcings since they control the mechanical strength and the bulk density of the geological layering. In this regard, the role of covers, which derive from vertical weathering profiles of the soil, in the predisposition to shallow landslides when triggering events, such as intense rainfall and earthquakes, are concurrent is not fully constrained. Such uncertainty arises from the objective difficulties in spatializing the thickness that the alteration covers have on the slopes as well as in the reconstruction of the geometries through which they are connected to the alluvial or debris at the valley floor.
The presence over large areas of pre-existing landslides can, in turn, play an effective role in facilitating new instability since residual strengths act within landslide mass. In this context, when a dynamic trigger acts, quantifying the concurrent presence/absence of specific predisposition is important to adequately assess landslide scenarios if associated with the spatial distribution of the intensity of the seismic shaking.
While conditioning of different sets of parameters (e.g., topographic) on seismic slope stability has been parametrically analyzed [
5,
6,
7] and the sensitivity of specific landslides often approached, the predisposing contributions of rainfall prior to or directly after an earthquake were rarely analyzed in literature [
8,
9,
10]; few are the complete and detailed landslide inventories referred to landslides induced by an earthquake under wet-to-saturated soil conditions. Rare are also deterministic studies dealing with the predisposing role of rainfall prior to an earthquake in controlling the spatial trend and distribution of coherent failures [
11,
12]. Among these studies, Jibson [
13] and Keefer [
14] reported cases of coherent slides occurring on a clayey slope, saturated by above-ordinary pre-seismic precipitations.
The evaluation of the importance of causative factors among rainfall-induced landslides (RfIL) and earthquake trigger landslides (EqTLs) should account for the hydraulic response of slope material, whose behavior relies on the local combination of lithology, nature of weathering, unsaturated permeability, rainfall intensity, local topography, and vegetation cover [
15,
16]. Because of the large requirements of input data and the sharp variability of pore pressure distribution, such analyses usually refer to specific landslides and are approached by numerical analyses and physically-based models [
17].
The reliable zonation of thickness and hydraulic properties of soil cover that can be mobilized under seismic shaking is also tricky. These uncertainties require the application of empirical techniques validated by ground evidence, aimed at reconstructing the distribution of the alteration covers [
18,
19], as well as reconstruction of physical initial soil conditions and mechanical parameters represents a challenging and fundamental goal for the analysis of landslide scenarios in a multi-hazard perspective. For this mutual combination of factors, deep insights into this topic and assessment of causes controlling landslide scenarios after earthquakes are not obvious and require probabilistic or simplified parametric approaches to address problems of data availability and uncertainties [
15].
Therefore, for contributing to the better understanding of how combined and transient destabilizing actions could induce slope failure over a wide landslide-predisposed area, we back-analyzed the landslide scenario that occurred after the 2018 Mw 5.1 Montecilfone earthquake (Molise, Italy), where the role of antecedent intense rainfall and the intrinsic tendency of the clayey slopes to evolve into landslides were highlighted by Martino et al. [
20]. To this aim, the main activities carried out to weigh the mutual influence of the pore-pressure and the seismic shaking in the stability conditions of the surficial soil covers will be presented in this paper:
Recognition of the soil cover thickness based on both field observation and literature review;
Field and laboratory permeability tests to obtain reliable permeability coefficient values for the different geotechnical units in which different soil covers were grouped;
1D infiltration numerical models to estimate the water content of the soil covers and water infiltration depth, spatializing the data over time in the entire study area, and reproducing the intensity of the rainfall event that occurred before the earthquake;
Slope stability analysis in the soil covers for conditions referable to the experienced rainy conditions or in the absence of rainfall.
Finally, a comparison between the location of first-time failure and reactivation of pre-existing landslides which reached critical conditions computed by seismic displacement regression model was performed with respect to the observed landslide validating the scenario directly surveyed immediately after the event.
2. The Montecilfone Earthquake Case Study
The Molise region is characterized by a very high frequency of landslides, which impact more than 28,000 sites over an area of about 4400 km
2, including a wide range of processes and landslide size [
21,
22,
23,
24], with a different state of activity (IFFI database;
https://idrogeo.isprambiente.it/app/iffi; last access 10 March 2022) [
25]. The high landslide concentration in the study area and the abundance of shallow landslides (i.e., soil slip), earth slides and solifluction are strictly related to the lithological features, where slope deformations mainly involve quaternary cover deposits and weathered portions of the clayey and marly bedrock. These units are ascribable to a thick Cenozoic marine succession belonging to the Lagonegrese-Molisan Basin [
26,
27], whose lower sequence is mainly composed of scaly clays [
28] with subordinate calciturbiditic intercalations and followed by Neogenic arenaceous-marly and calcareous-marly flysch deposits [
26]. In the North-East sector of the study area Plio-Pleistocene blue marly-clays belonging to the lower horizon of the Mutignano Fm. [
29] crop out, often involved in landslides.
In this landslide-active context, on 16 August 2018, an Mw 5.1 earthquake occurred at 18:19:04 UTC in Montecilfone (CB), Molise region (Central Italy), anticipated by a seismic event of the 14 August 2018 with a Mw 4.6 (SISMIKO Group INGV), striking an area involved in the 2002 San Giuliano di Puglia seismic sequence [
30], located about 20 km south of Montecilfone (
Figure 1a).
The mainshock of the 2018 seismic sequence caused a peak ground acceleration (PGA) of up to 0.12 g (
http://shakemap.rm.ingv.it; last access 10 March 2022), which corresponds to a return time of 745 years according to the Italian seismic hazard map (
http://esse1-gis.mi.ingv.it/; last access 10 March 2022). The earthquake reflected a reduced macroseismic intensity (V–VI; EMS98 [
33]) but induced an unexpected frequency and distribution of ground effects [
20]. Ground crackrelated to landslides were observed within an area of tens of square kilometers, resulting in an induced scenario that suggested a possible combined action of multiple destabilizing actions controlling the landslide scenario. The seismic event was anticipated by intense rainfall that affected the NE part of Molise region, causing a cumulative rainfall of up to 160 mm [
20], with an intensity of 145 mm cumulated in three days in the epicentral area, recorded at the Palata weather station (
Figure 1b).
Field surveys carried out in the days immediately after the earthquake (
Figure 1a) allowed us to recognize and classify 88 earthquake-induced ground effects, grouped into 75 coherent slides (earth slide and earth flow mechanisms), 9 disrupted landslides [
34] (
Figure 1b), and 4 ground cracks mostly related to slides.
Such landslides occurred in the summer, after more than 30 days of absence of intense rainfall, in a period characterized by the landslide’s quiescence [
35]. A multitemporal analysis performed by DIn-SAR interferometry pointed out that the recurrent landslide activity in the area mostly falls in the autumn–winter period up to early spring, where intense rainfall causes intermittent reactivations [
35]. In fact, many of the surveyed earthquake-induced ground effects (up to 37.5%) consisted of reactivated landslides, some of which were already catalogued in the official Italian inventory (IFFI). Instead, the largest portion corresponds to first-time slope failures that involved the clayey slopes (60%), with minor landslides in marls and arenaceous flysch.
While the totality of disrupted landslides occurred at a distance lower than the Keefer [
34] and Martino et al. [
36] distance for disrupted landslides, 43 out of 75 coherent landslides occurred beyond the maximum expected threshold given the Mw 5.1 of the earthquake. Such surveyed scenarios pointed out a specific proportion of first-time failure and reactivations, whose location and relationship with the intensity of the two triggers suggest a mutual contribution of ground shaking and rainfall in the observed evidence (
Figure 2; [
20]). In this context, this research aims at back analyzing the observed landslide scenario, considering the joint contribution of soil thickness distribution, physical state conditions of the unsaturated covers, and the local unbalancing forces provided by the seismic trigger.
3. Materials and Methods
To better frame the predisposing or triggering role of antecedent rainfall in aggravating the landslide scenario and to evaluate potential interference among the two factors, an aerial GIS-based stability analysis was performed adopting a pseudo-dynamic Newmark displacement model by means of Hsieh and Lee [
37] regression for soil materials. For the definition of critical acceleration, soil saturation conditions and suction stresses in the soil covers were accounted for in the static safety factor (SF) [
38].
To achieve these goals, ancillary data, such as geological and landslide inventory maps or available geotechnical parameters of the outcropping lithologies, were primarily collected. Four pre-existing databases (DBs) were considered, i.e., IFFI landslide catalogue [
25], Molise region microzonation studies (MS), geo-lithological DB by Forte et al. [
39], Geoservizi S.r.l. DB, and a 1:100,000 geological map of Molise [
26]. The reference target for validation in the back-analysis of the observed scenario consists of the ground effects induced by the Montecilfone earthquake, where attributes and locations can be downloaded or loaded as a Web Map Service (WMS) in open-access from the new release of the Italian Catalogue of Earthquake-Induced Ground Effect (CEDIT) (
https://gdb.ceri.uniroma1.it/; last access 10 March2022 [
32]).
Afterwards, direct field surveys devoted to soil sampling and altered cover characterization were conducted to fully characterize physical parameters and state conditions required for analytical stages. The workflow adopted to manage the multihazard scenarios consequently includes (i) laboratory tests for geotechnical soil cover characterization; (ii) spatialization of soil thickness measurements and pluviometric data generalization by definition of Thiessen polygons; (iii) 1D time-dependent modelling of rainfall infiltration in the covers through the use of Hydrus-1D open-source software; and (iv) computing of the SF in static dynamic conditions with a hydraulically controlled scheme. Technical aspects and analytical solutions adopted are reported in the following.
3.1. Soil Cover Characterization
Most of the observed landslides in the study area involve eluvial-colluvial soil covers. To characterize the presence and distribution of the weathered covers developed on the different lithologies, more than 50 direct measurements of total cover thickness were collected, while the maximum visible thickness in natural and man-made cuts was surveyed at 264 points to constrain minimum soil cover thickness (
Figure 3). Available stratigraphic logs performed in the frame of microzonation studies were also analyzed to enrich the DB and distribution of soil cover thickness.
To define the spatial distribution of the soil thickness in the study area and constrain the potential volume mobilized by shallow landsliding, two empirical models of the soil cover thickness proposed by Saulnier et al. [
18] were adopted. These simplified models, calculated in a GIS environment, compute the local thickness according to the local slope angle (S-model) and elevation (Z-model) through the linear correlation hereafter reported:
where
is the thickness of the soil cover in the i-cell (m);
is the maximum thickness of the soil cover in the area (m);
is the minimum thickness of the soil cover in the area (m);
is the elevation of the i-cell (m a.s.l.);
is the maximum elevation of the study area (m a.s.l.);
is the minimum elevation of the study area (m a.s.l.);
is the slope angle in the i-cell (°);
is the maximum slope angle in the study area (°), and
is the minimum slope angle in the study area (°).
These two models were merged into one, called the SGZ model (
Figure 3), which was created with the aim of merging the weights of the two factors in the development of the alteration covers and incorporating the contribution given by the different lithologies, adopting a simplified relationship between them:
where the multiplicative factor
incorporates the contribution of the lithology to the thickness of the soil cover. This parameter is defined as the ratio between the average of the thicknesses referring to each lithology, and the maximum thickness of the lithology outcroppings in the area:
where
is the lithological factor of the l-lithology in the i-cell (-),
is the thickness of the l-lithology (m) and
is the maximum thickness among the lithologies outcropping in the area (m).
By using the SGZ model, a soil cover depth map was obtained (
Figure 3) which was useful for conducting the slope stability analyses presented below. The reactivated landslides are also reported on this map, whose sliding surface was deduced according to the semiempirical relation with the volume proposed by Guzzetti et al. [
40].
Soil covers sampling, in situ permeability, and density tests were conducted directly on the field following standard procedures [
41]. Laboratory tests were also performed in the Engineering Geology Laboratory of the Earth Sciences Department of Sapienza University of Rome on the sampled soil covers to estimate geotechnical parameters and attribute a class of the USCS (
Table 1) to each collected sample. In particular, the following tests were carried out:
Grain-size analysis (both sieving and sedimentation process; ASTM D422–63 standard);
Liquid limit, plastic limit, and plasticity index of soils (ASTM D4318–10 standard);
Specific gravity of soil solids (ASTM D854–10 standard);
Laboratory compaction tests (ASTM D698–07 standard);
Permeability of granular soils (ASTM D2434–68 standard).
Four different geotechnical units for the soils covering the lithologies cropping out (see
Figure 1) were defined based on the resulting geotechnical characteristics:
Sands and conglomerates (1): characterized by calcarenites (Pliocene) with a significant granulometric fraction in the gravel sizing (about 15%);
Clays, marls, and sands (2): mainly represented by the weathered covers developed onto the Plio-Pleistocene blue clays;
Clays and marls (3): related to the marly clays formation (Pliocene), belonging to the piggy-back basins in the hanging wall of the main thrust, and the scaly clays (Miocene-Cretaceous);
Clays, marls, and limestones (4): which includes the Paleogenic flysch (Messinian-Serravallian) and the gypsum formation (Messinian).
Classification and properties of the sampled soil covers are summarized in
Figure 4 and in
Table 1.
3.2. 1D Time-Dependent Model of Rainfall Infiltration
Hydrus-1D is modelling code for the analysis of water flow and solute transport in an unsaturated medium [
42]. In this study, the code was used to estimate the saturation of the alteration covers and the maximum depth of infiltration of the wetting front during the pluviometric event that involved the study area from 14 August 2018 to 17 August 2018 (
Figure 5). The process was carried out for each lithology, whose permeability parameters were constrained by in situ and laboratory tests, detailing the cumulative rainfall over time recorded in each pluviometric station.
In particular, the equation proposed by Van Genuchten [
43] was adopted, which allows the calculation of the soil water retention curve as:
where
is water retention (dimensionless);
is the suction (L);
is the saturated water content (dimensionless);
is the residual water content (dimensionless); and
and
as fitting model function of pore size distribution.
To predict hydraulic parameters for each unit we take advantage of the embedded ROSETTA Lite software [
44], which relies fundamentally on a database of measured water retention and other properties for a wide variety of media. Starting from the soil grain size distribution, the code is able to estimate the
and
parameters, together with saturated (
) and residual (
) water content, for reconstructing the water retention curve (i.e., the relationship between soil water suction and the amount of water remaining in the soil) assuming the van Genuchten’s model [
43] (
Table 2).
3.3. Static, Hydraulic Controlled and Coseismic Cover Stability
To compute the SF in the soil cover and consider landslide processes featured by shape ratio (L/H > 10) (i.e., mainly shallow landslides), we used a modified version of the preceding equation which was proposed by Lu and Godt [
38] to account for partially saturated soil conditions:
where
is the suction stress;
is the internal peak or residual friction angle (°);
is the cohesion (kPa);
is the unit weight of soil cover (kN/m
3); and
is for slope angle (°).
The
values have been estimated for each geo-lithological unit in the laboratory and referred to literature values reported in Giannini et al. [
45] and Forte et al. [
46]. As the value of H, we considered and fixed different depths within the soil cover to assess the in-depth distribution of SF and weight the contribution of capillarity in the slope stability.
Regarding the mechanical parameters, peak strength conditions were respectively assumed for slope portions that did not experience instabilities in the past while residual conditions were adopted within the inventoried landslide perimeters (i.e., already occurred in the past).
In the case of shallow landslides, for the calculation of
, which accounts for the capillary forces that counterbalance the shear stress in unsaturated soils [
47], we assumed the relation formalized in Lu et al. [
48]:
where Sr is the degree of saturation, defined as in Equation (8) and α and
represent the fitting parameters of the Van Genuchten’s [
43] equation reported above (
Table 2).
The value of Sr was obtained according to the following equation:
where
is the water content (-);
is the water content in saturated conditions (-), and
is the residual water content (-).
To calculate the SF values in seismic conditions, the Newmark’s method was applied, based on a yielding acceleration (K
c) threshold calculated firstly by a limit-equilibrium pseudo-static analysis through the relations:
respectively used if acceleration acts parallel to the slope or horizontally.
According to Romeo [
49], the analysis referred to a simplified solution neglecting the vertical component of the ground motion. Coseismic displacement was computed adopting the semiempirical regression model proposed by Hsieh and Lee [
37] for cohesive soils, which upgrades the original formula by Jibson [
50,
51] in the following form by considering a worldwide EqTLs dataset:
where
identifies the Arias intensity (m/s) and
the critical acceleration (m/s
2).
was defined basing on the PGA through the Romeo [
49] empirical equation:
where PGA is the peak ground acceleration (m/s
2).
In agreement with Romeo [
49] and reference therein, the critical displacement threshold defined for assuming a generalized failure for coherent landslides, Dc, was set at 10 cm.
A peak friction angle (peak resistance) was assumed for first-time soil failures, while a residual friction angle (residual resistance) was considered as admissible strength for landslide reactivation.
4. Results
4.1. Modelling of Rainfall Infiltration
The 1D numerical models for the soil infiltration ran for each of the soil columns representatives of the different rainfall inputs and the various lithologies involved. Using the lab-derived parameters and the physical conditions of the considered soils, the degree of saturation of the covers was attributed considering its change over time and simulating the increases over time and the decreased over-depth of water content.
As exemplified in
Figure 5, the results of the modelling by Hydrus-1D, demonstrate that the saturation of the soils increased over time, simultaneously with the progress of the rains, until 16 August 2018 (end of T3), i.e., when the rainfall ceased and the Montecilfone earthquake occurred. At the rainfall peak, the degree of saturation no longer increased due to the absence of rainfall intake. The saturation decreases in the most surficial layers of the soil covers due to the infiltration process that progresses and simultaneously increases the saturation of the deeper layers.
Figure 6 shows the degree of saturation evolution as a function of time and depth, spanning from the beginning of the rainfall event to the day after the Montecilfone earthquake when the rainfall ended.
At a depth of 20 cm b.g.l., the degree of saturation increased until the 16 of August and decreased from the next day, thus showing a reasonably direct correspondence with the pluviometric input suffered by the study area. Considering a depth of 40 cm b.g.l., the degree of saturation began to increase after the second day of rain, while the previous days did not present surficial inflow. This step also corresponds to the end of the rainfall infiltration into the soil covers. Due to infiltration processes, the degree of saturation at this depth continued to increase even after the end of the pluviometric input.
Finally, at a depth of 60 cm b.g.l., the degree of saturation slightly increased in limited portions of the study area, where the most permeable lithotechnical units crop out, only beginning from the day immediately following the end of the rains. At the end of the third rainy day, the Montecilfone earthquake occurred, shaking the soil covers already saturated (with different saturation degrees) by the previous rainfall.
4.2. Slope Stability Analysis in Evolving Saturation Scenarios
The slope stability under static conditions was analyzed considering the different saturation scenarios in varying time windows, i.e., from 14 August 2018 to 17 August 2018, and at different depths, i.e., from 20 up to 60 cm b.g.l., as well as evaluating the overall stability of the entire thickness of alteration covers, ranging over the area between 20 cm and 1.6 m (
Figure 7). The slope stability analysis was performed applying the general formulation proposed by Duncan and Wright [
4] and the Lu and Godt [
38] relation for assessing suction stress, which considers the stabilizing contribution of the negative pore pressure in unsaturated conditions. This contribution appears to be relevant considering that the soil covers rarely and locally reach full saturation conditions in the considered study area.
Three loading conditions were considered in the slope stability analysis within the shallow portion of the cover (20–40–60 cm b.g.l.), which respectively account for: (i) the full saturation conditions with breakage of the capillary fringes; (ii) saturation values as a function of time and space in accordance with the 1D infiltration profiles obtained by Hydrus numerical modelling; and (iii) depth higher than 60 cm, where a minimal or negligible increase in saturation can be assumed and in which only the frictional contribution function of the effective stresses must be considered.
The static SF values obtained by the slope stability analysis at 20–40–60 cm b.g.l. reflect the different saturation profiles and pressure heads in the soil cover (
Figure 7), whose stability is directly dependent on the counterbalancing of suction stabilization and soil weight destabilization. For constant depth and lithotechnical properties, values of the static SF increase as the suction stress increases (more negative), as it results in an effect of “apparent cohesion” in shallow cover stability [
52,
53]. Minimum SFs were maintained where the suction stress values were lower. By increasing the time, the static SF decreases weakly with increasing rainfall from T1 (14 August 2018) to T2 (15 August 2018). After T2, the water content remains unchanged, especially at 20 cm depth, and there are no more clear variations in the value of the static SF. By increasing the depth, the reduction of the static SF become more visible (60 cm depth in
Figure 7).
The slope stability analysis under static conditions was also carried out considering the whole thickness of cover as generalized for the total extension of the study area, assuming the potential mobilization of the entire alteration cover. In this analysis, the suction term was considered unchanged, and the weight volume unit was considered both in saturated, natural, and dry conditions imposing the respective γ of the analyzed soil. From the results of this analysis (
Figure 8), it emerges that the rainfall played an aggravating role in stability conditions, worsening the SF as a function of depth, time, and cumulative rainfall. Minimum SFs were reached at the bottom of the soil cover, where an SF close to 1 was computed.
The slope stability analysis under static conditions highlighted the contribution of the rainfall from two different points of view: at the shallowest levels (i.e., 10 cm b.g.l.) the soil cover thickness is not sufficient to establish slope instability conditions since the destabilizing actions, determined by the lithostatic weight, are too low despite the degree of saturation imposed being close to 100% and the capillary fringes are broken. This fact leads to general stability conditions over the entire area.
At 20 cm b.g.l., the calculated SFs are quite high over the whole area, but not sufficient to justify slope instability. Furthermore, the minimum soil thickness does not contribute much to the decrease in the SF. Combined with a minimum stabilizing contribution given by the residual suction in areas that do not reach a complete saturation, these conditions lead to general slope stability over the whole area. At 40 cm b.g.l., the SFs are, however, still quite high but the contribution of the greater soil cover thickness leads to a decrease of the slope stability conditions. The suction is, in any case, present and stabilizing since the capillary fringes are effective. At 60 cm b.g.l., maximum modelled depth of rainfall infiltration there is little presence of infiltration water and the SF decreases. At 80 cm b.g.l., the soil covers are always unsaturated since infiltration never reaches these depths. At greater depth, the SF noticeably decreases, remaining above the unit and reaching minimum values assuming the mobilization of the entire thickness of the soil covers (
Figure 7 and
Figure 8).
Comparing the SF maps in dry and rainy conditions (
Figure 9), obtained assuming the mobilization of the entire thickness of the soil covers, it is possible to see how the infiltration tends to lower the value of the SF, bringing some areas close to the unit value. The presence of infiltration water would have weighed approximately uniformly over the entire study area, reducing the SF by decimals up to units depending on the wetting front depth reached (
Figure 9). However, the areas of maximum rainfall and greater infiltration depth are found in the zones with a lower PGA (peripheral areas), confirming that the rainfall has played mostly an aggravating role in the reduction of the stability conditions but did not act as a trigger.
4.3. Coseismic Displacement Computation
The worsening of the soil cover stability conditions under the action of the seismic shaking was measured by calculating the cumulative coseismic displacement according to the Newmark method [
54], assuming a threshold of 10 cm to evaluate the exceedance probability for the critical pseudo-static coefficient (
Figure 10). The first iteration involved calculating the cumulative coseismic displacement in conditions of peak shear strength, i.e., assuming as strength parameter, the effective friction angle of the soil covers (φ’). The results of this analysis were compared with the first-time CEDIT landslides, returning a back-analysis accuracy of 27.66% (13/47). The second iteration was carried out considering the residual friction angle (φ
r) and assuming lowered shear strength accounting for the potential of reactivation of existing landslides. The model returned a 75% accuracy prediction (21/28) with respect to the reactivated landslides inventoried in the CEDIT catalogue.
5. Discussion
The obtained results show that the intense rainfall was insufficient to lead to the critical conditions (
Figure 7,
Figure 8 and
Figure 9), showing that SFs at different depths tend to decrease with the advance of precipitation. Saturation involved the first 60 cm of weathered soil covers, not resulting in localized instability. In this wetting front, except for the first couple of centimeters where the full saturation is reached, the suction stress played a stabilizing role against the developed shear stress. At greater depth, below the max depth of infiltration, the SF decreased as the load increased. Minimum SF slightly greater than 1 were reached at depth equal to the whole soil cover column (
Figure 8). Therefore, despite being assumed in hydrostatic conditions, the rainfall was not sufficient to trigger landslides and a mainly predisposing role can be attributed to it.
Regarding the coseismic displacement analysis, it was observed that the earthquake was actually sufficient to induce coseismic displacement greater than the critical threshold of 10 cm (
Figure 10), considered in literature enough to assume instabilities for flow and slide mechanisms in cohesive soils [
49,
55,
56].
Referring to the DB of induced ground effects surveyed in the CEDIT catalogue, the back analysis shows good performance resulting in critical coseismic conditions for 13 out of 47 first-time failure locations and 21 out of 28 points where landslide reactivations were registered (
Figure 10). Most of the landslides predicted in the here modelled scenarios fall into the Keefer distance for coherent slides where the maximum seismic load occurred (
Figure 11—rainy conditions), however several landslides have been activated at greater epicentral distances where the maximum rainfall was registered.
Limitations in the ex-post reconstruction of landslide occurrence can be found in several critical aspects, most of which rely on the large-scale approach adopted and local uncertainties in either the soil cover model, in the extrapolation of rainfall intensity in Thiessen regions, in geotechnical parameterization, or in combination.
The displacement map reported in
Figure 10 clearly shows the good fit with crowns of already inventoried landslide perimeters, marking the initiation area of earth slides and the accuracy of the performed inventory. The diffuse instability conditions retrievable from coseismic analysis with respect to the survey ground effects pointed out a tendential underestimation of the number of EqTL reactivations since the evaluation of reduced coseismic displacement in existing landslide mass is impractical by direct survey over large areas. Alternatively, an underestimation of the residual strength in landslide mass may have been included, however, a deterministic approach for specific cases is not feasible in large scale approaches. The lower back-analysis performance emerges from first-time landslides, whose spatial distribution clusters are in the epicentral area and in the near-fault sector, where site effects cannot be excluded a priori, and transiency of pore pressure is possible. As it regards the role of forcings, the comparison between the effective landslide scenario and the one predicted by Newmark’s approach reported in
Figure 11 is more suitable in cases of rainy conditions with respect to the modelled results under the hypothesis of absence of antecedent rainfall. This proves the non-negligible effect related to the rainfall saturation of soils for the earthquake-induced shallow landslides and further demonstrates the multi-hazard influence in the landslide triggering.
Considering these critical features, the adopted workflow has proved to be useful for the reconstruction of multi-hazard scenarios and the quantification of the role of causative factors in seismically induced landslides in earth slopes. The simplified approach, despite being based on empirical and regression models and supported by the limited available site information, provided a reliable scenario, proving its applicability for the analysis of the impact of a seismic event. This approach integrates with the regional-scale predictive approaches available to date [
45,
57,
58], and are thus useful for supporting civil protection activities and emergency management duties.
6. Conclusions
The results obtained demonstrate the suitability of the applied analytical approach for the reconstruction of seismic landslide scenarios, at basin scale, in hydraulically controlled conditions (multi-hazard conditions).
As the main outcomes, the performed analyses highlighted how concatenated or combined factors can result in more severe landslide scenarios than the one expected assuming the single contribution of the seismic trigger, especially in highly landslide susceptible areas. The antecedent rainfall by itself, despite it lowering the static SF, would not have been able to induce the first activations of landslides in the area, which can thus be considered seismically induced. However, the latter occurrence results in an unexpected event, according to the low magnitude of the mainshock, strengthening the hypothesis that the landslide scenario involved slopes that are predisposed to shallow landslides and prepared by saturation of soil covers after rainfall occurrence.
Low-to-middle magnitude earthquakes can thus reflect increased risk conditions over large areas, therefore, even if such combinations are rare, they must be considered for land management and in regulatory technical products. Combinations of hydraulic conditions and parameters related to ground motion should be adequately considered in large scale evaluation of EqTL scenarios.