A Methodologic Approach to Study Large and Complex Landslides: An Application in Central Apennines

Abstract

The evaluation of landslide hazards in seismic areas is based on a deterministic analysis, which is unable to account for various uncertainties in the analysis process. This paper focuses on the probabilistic local seismic hazard analysis and extends the results to the landslide hazard analysis to consider both the uncertainties of the ground deformations and the strengths. The work studies the areas between Nibbiano and Sant’Erasmo hamlets in the Camerino municipality located in central Italy, where all constructions present evidence of damage caused by both the seismic sequence of 2016–2017 and the slope instability. An exhaustive geological and geophysical investigation has clarified the geological, geomorphological, and hydrogeological characteristics of the area, enabling a new characterization of material stress-strain behaviour. The study reveals that the low stiffness of the debris covers, and their fair degree of permeability contribute to potential instability scenarios triggered by both intense rainfall and the effects of strong earthquakes. The goal was to utilize the results to support local urban planning because in-depth knowledge of the possible evolutionary scenarios of the slopes is fundamental to the management of the degree of danger for structures, especially for people. Moreover, it was shown once again how a multi-source approach, with different investigation techniques, cannot be ignored for the study of the evolution of complex landslides.

1. Introduction

Gravitational modelling represents one of the main factors of change on the Earth’s surface to the exclusion of tectonically induced deformations. It is recognized in the literature how this can act at different scales, thus moving from very broad phenomena that are generally slow to more localized ones. The present work is set in a “wide” context (sector of a ridge between two rivers). Therefore, the bibliography on DSGSDs (Deep-seated Gravitational Slope Deformations) and large landslides were taken into consideration [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20]. Gravity-related morphogenetic processes, therefore, also play a major role in shaping the reliefs in the Apennine-Adriatic sector of central Italy [21,22,23,24,25,26,27,28,29,30,31,32,33]. Here, the land is subject to recent uplifts and is widely affected by deepening linear erosion and, thus, high relief energy values [34,35,36]. Under such conditions, mass movements of different types and sizes are triggered, depending on the litho-structural, geomechanical, and hydrogeological characteristics of the bedrock, the nature of the eluvial-colluvial cover, weather and climate conditions, land use, and anthropogenic changes to the area [36]. The erosion that is exerted at the base of slopes by watercourses is one of the main factors in gravitational morphogenesis. Landslide movements, even large ones, have always accompanied the phases of rapid river downcutting linked, from time to time, to the return of biostatic conditions in the general deepening of valley systems. In recent times, this latter process has been particularly intense in the medium-low stretches of the main river courses as a result of intense human excavation in the riverbed. The incision propagated in the proximal sections of the minor hydrographic network has favoured the reactivation of ancient landslide phenomena and/or the activation of news [37,38]. Another important trigger of gravitational movements in the study area is tectonic activity and the recurring earthquakes associated with it [39,40,41,42,43]. Seismic shocks can, in fact, induce unstable conditions on slopes and favour the triggering of mass movements through various mechanisms, summarized as follows. Cyclic variations in interstitial pressures that, in fine sands and water-saturated silts, can give rise to liquefaction phenomena, oriented accelerations, surface tectonic faulting and fracturing, fracturing of rock masses, and disconnection of semi-coherent debris deposits. Also relevant to slope stability may be co-seismic modifications, even minor, of slope angles and underground water tables. Finally, extreme climatic events and, more generally, the recurrence of intense and prolonged rainfall are important. In this context, the increase in summer storms in recent decades, when cracking in clay soils is more intense and favours greater infiltration, can certainly be a favourable factor for landslides [44,45,46].

Gravitational movements can develop over very different time intervals (up to the order of thousands of years), with also very different return times (depending on the recurrence of meteorological or seismic triggering factors, considering apart human activity). Bearing in mind that a phenomenon can be defined as inactive if it is no longer capable of being reactivated in the present morpho-climatic system [47], it can be stated that the phenomena with characteristics such as being identified as inactive are extremely few.

Studies have shown numerous cases of interference, in ancient or recent times and sometimes on several occasions, of mass movements on river dynamics with the production of dams and diversions. Phenomena of this type have been recognized in recent decades in the high hilly stretches of the area, where landslide damming of watercourses has resulted in short-lived lake episodes due to easy threshold erosion. In the narrow valleys of the Apennine chain, the recurrence of river barrages in the Upper Pleistocene and Early Holocene is evidenced by lacustrine deposits and landslide accumulations alternating with alluvial gravels [48,49,50,51]. Geognostic investigations are still fundamental and irreplaceable tools for supporting surface geological surveys and morphological interpretations of landforms. The continuous core drilling technique, for both soils and rocks, is the most suitable direct investigation method for defining the lithologies traversed and the stratigraphic relationships present in the subsurface. The most common limitation of these investigations is obviously represented by the possibility of obtaining information only on the vertical investigated. For this reason, and because it is often technically and economically impossible to drill large numbers of holes, it is increasingly common to integrate them with geophysical techniques. Together with core drilling, in fact, surface geophysical prospecting represents an excellent tool for a continuous view of stratigraphic relationships, the presence of anomalies, and any buried discontinuities. In recent years, the significant development of hardware and software technology has led to a considerable improvement in the interpretation of the propagation of seismic waves, both volume (P and S waves) and surface (Rayleigh and Love waves). Their application in the study of the evolution of large landslide phenomena enabled the identification of mobilized volumes and sliding surfaces even at great depths [52,53]. The main advantage of integrating the various investigation techniques is the correlation of point and/or local information with spatial information of both geometries and potential strength of materials, which are crucial for an accurate definition of landslide models. In the study of landslide movements, the application of limit equilibrium criteria is increasingly an outdated tool. The application of the method provides only partial answers to the problem, as there is no information on the deformations acting in the soil and mainly on the changes in shear resistances around the sliding surface. In this analysis, displacements are not considered, assuming that the acting and resisting forces are independent of deformation. This corresponds to a perfectly plastic behaviour in which the resistance remains unchanged after fracture [54]). Moreover, in most cases, the analyses are constructed based on standard landslide or homogeneous soil models [55]. This is quite different from real landslide conditions where the heterogeneity of the materials involved and their mechanical properties, which often vary even over short distances, are the main predisposing factors for instability. In recent years, the use of advanced numerical codes has become increasingly common to support hypotheses on the gravitational evolution of slopes as well as control factors, particularly for large-scale phenomena. Numerical analysis, in fact, makes it possible to simulate the possible triggering mechanisms and assess their potential evolution by varying the boundary conditions and considering possible reductions in internal resistance. The main limitations of a numerical model are due to the uncertainties of the data entered as input both for the geometric reconstruction of the model itself and for the geotechnical parameterization to be used in the chosen constitutive models. Based on the research paths summarized above, the present study, therefore, combines various methodologies useful for better interpreting large and complex mass movements, such as the one investigated, also with the aim of identifying and fine-tuning an adequate methodological approach for this type of phenomenon.

More specifically, this work shows the study and detailed characterization of a large gravitational phenomenon located on the eastern side of the Umbria-Marche Apennine between the Chienti and Potenza rivers (Figure 1).

This Apennine part is characterized by deep gravitational deformations and large landslides, which have had a complex geological and geomorphological quaternary evolution [56,57,58,59]. In the areas of overthrust, large landslide phenomena have occurred over time, with mechanisms of collapse, overturning, and sliding (translational and rotational), which have produced powerful landslide deposits. The objective was to analyse the geomorphological evolution of the gravitational phenomena and the relation between the geological-structural setting and the landslide mechanisms. Furthermore, the aim is to understand the correlation between the bedrock and the gravitational phenomena, as well as to consider the high seismicity of the area. The Nibbiano landslide is the one that best suited our purpose, given the emptying of the large Paleo-landslide accumulation and the bedrock that is not too deep. This setting has facilitated the interpretation of the landslide debris/rock contact. The area was affected by the seismic crisis of 2016, which caused the partial reactivation of various gravitational landslides, including the Nibbiano hamlet.

The current work was drafted according to the following description. To reconstruct the current geological and geomorphological conditions of the area, a search and a collection of information from the current bibliographic sources were accomplished. Following this, a geognostic and geophysical investigation was conducted to study the stratigraphic trend and obtain geotechnical information in depth. Finally, numerical simulations were used to identify the activity and the possible evolution of the landslide phenomena. The calculations were realized with 2D and 3D models to study both the static and the dynamic conditions. The results were used to identify the areas to be relocated because the degree of hazard does not allow their future use. Moreover, the main scientific contribution of the research was to support the local administrations for the purposes of urban planning in the area under investigation.

2. Geological Features

The Umbria-Marche Apennines are represented by an arcuate mountain ridge that stretches for about 450 km and has very pronounced reliefs, with altitudes above 2000 m above sea level, compared to the adjacent foothills. The Apennine reliefs are mainly made up of pelagic, Jurassic-paleogenic carbonate sediments, deformed by folds and overthrusts associated with Neogene compressive tectonics. The Apennine structure has been classically interpreted according to a model of pellicular tectonics, with imbrication of sedimentary units that have dislodged along Triassic evaporites on an undeformed and buried crystalline basement [60,61,62]. Quaternary extensional tectonics subsequently disarticulated the edifice in folds and overthrusts with normal faults dipping predominantly towards the WSW that created tectonic basins, even large ones (Figure 2). The examined area is located on the western edge of the Camerino syncline and in the foothills affected by the front of the Valnerina—Monte Igno overthrust (Figure 2). A detailed geological-structural survey was carried out along the Monte Igno overthrust, including a significant area with respect to the Nibbiano landslide phenomenon (Figure 3).

The Monte Igno anticline shows the eastern flank overturned towards the east with predominantly calcareous, calcareous-marly, and marly rock formations (Cretaceous-Paleogenic Group) that overlap through one or more overthrusts, the predominantly marly-calcareous, marly-clayey and arenaceous formations of the Miocene Group. The Cretaceous-Paleogenic Group includes the Marne a Fucoidi (marly limestones and marls), the Scaglia Bianca and Scaglia Rossa (marly limestones and limestones), and the Scaglia Variegata and Scaglia cinerea, the latter consisting mainly of marly-clayey and clayey formations.

The formations of the Miocene Group include the Bisciaro (marls and marly limestones), the Schlier (marls, calcareous marls, and clayey marls), and the Camerino Sandstones, consisting of silicoclastic turbiditic deposits (Figure 4a,b). The geological cross-section in Figure 4 highlights the lithological and structural characteristics of the Nibbiano landslide bed. In the summit part, roughly corresponding with the detachment niche of the landslide, the contact by overthrusting is realized between the Scaglia Rossa and Scaglia variegata formations with the underlying Scaglia cinerea, the latter outcropping for a wide portion. All the formations involved are characterized by overturned layers, particularly fractured and with pervasive shear zones in the more marly lithotypes. Moreover, the geometric analysis of the numerous thrusts in the areas not affected by landslide phenomena highlights slopes of 30–40° degrees plunging towards the west [57]. The covering colters consist mainly of debris deposits and landslide accumulations, characterized by heteromeric materials ranging from medium grain sizes to blocks. They have variable friction angles, very low cohesion values, and often a high water content.

3. Geomorphological Features

On a large scale, the hilly system of the area is narrow between the reliefs that border the Camerino basin and present morphological features fragmented by a dense network of fluvial incisions with colters of varying thickness. The high energy of the reliefs to the west and the mechanical characteristics of the often-saturated soils have contributed to generating important landslide movements with different states of activity, in many cases also connected to the seismic conditions of the area and the extreme events [27,33,46,63].

In particular, the study area extends along the eastern slopes of the high-hill alignment between the peaks of Mountains Igno (south) and Primo (north) and, respectively, between the Chienti and Potenza watercourses (Figure 1).

Along this side of the ridge, even from the middle of the slope, there is widespread debris with thicknesses of up to 30–40 m, which in many places mask the outcrop of the thrust plane. This mountain strip is also characterized by several springs, some of considerable discharge, which attest to the emergence of important aquifers.

As a result of the above, the following scenario strongly predisposing to instability is outlined.

• The steep eastern flank of the anticline has the presence of thrust planes, steeply inclined and overturned strata in the upper part, and gently dipping stratification in the lower part, sometimes out of the slope.
• Presence of extensive and thick debris accumulations resulting from degradation and local fall-type landslides.
• The situation of local highly saturated deposits.

To better understand the entire deformation phenomenon, as well as the complex relationship between the lithological-structural structure and the morphodynamics of the extended slope, the Nibbiano-Sant’Erasmo area was identified. Here, geomorphologically, the territory appears as emptied, testifying to the considerable activity of mass movements, past and present, as also noted by the authors in comparable geomorphological contexts [27,31,63].

Going into detail, both settlements are in the medium-low portion of the slope of the ridge, at altitudes of 646 m a.s.l. and 560 m a.s.l., respectively. They are located on the aforementioned detrital deposits, with varying thicknesses also in relation to the lithological and stratigraphic characteristics that have significantly influenced the morphological structure of the landscape and the development of the settlements.

A morphometric analysis of the system was indispensable in the study of landforms. The aim is to superimpose a quantitative evaluation of their characteristics on the description of the forms. Through morphometry, therefore, it is possible to develop mathematical models that effectively contribute to outlining the evolution of the relief [64,65,66,67]. To increase the accuracy of the analysis, a detailed morphometric study was carried out by acquiring LIDAR data (MASE. Creative Commons—Attribution—4.0 International) with a ground resolution of 1.0 × 1.0 m. The data was appropriately processed with specific algorithms to obtain a general view of the landforms, even in those parts currently occupied by residential buildings and/or road infrastructure.

Among the fundamental tools of analysis is certainly the use of the acclivity map, which expresses an initial classification of the territory according to the average slope of the slopes by analysing geometric factors of length and height. The analysis represents a set of techniques useful for quantitatively describing the morphology of places by calculating the average slope of the grid of listed points extrapolating them according to various degrees of inclination using the geometric information of the shapes. Slope inclination represents a predisposing factor of considerable importance for slope stability since it is directly related to the inclination of possible failure planes and/or horizons and, therefore, correlates with the distribution of landslide phenomena.

Using the digital terrain model (DEM), an analysis of the acclivity and energy of the relief can be carried out by calculating the gradient of the plane tangent to the surface in the direction of maximum slope. This is equivalent to the first derivative of the function expressing the elevation change along the same direction obtained for each cell or pixel.

The result is that which corresponds to the maximum gradient of the 3D surface section analysed and is therefore considered as the gradient of the maximum value according to the following formula.

$$Slope=\sqrt{{\left({\displaystyle \frac{\partial z}{\partial x}}\right)}^2+{\left({\displaystyle \frac{\partial z}{\partial y}}\right)}^2}$$

(1)

$$\left({\displaystyle \frac{\partial z}{\partial x}}\right)=\left[\left(z_{i+1,j+1}+{2z}_{i+1,j}+z_{i+1,j-1}\right)-\left(z_{i-1,j+1}+{2z}_{i-1,j}+z_{i-1,j-1}\right)\right]/8\partial x$$

(2)

$$\left({\displaystyle \frac{\partial z}{\partial y}}\right)=\left[\left(z_{i+1,j+1}+{2z}_{i,j+1}+z_{i-1,j+1}\right)-\left(z_{i+1,j-1}+{2z}_{i,j-1}+z_{i-1,j-1}\right)\right]/8\partial y$$

(3)

The formula number two represents the gradient in the East-West direction, while the formula number three represents the gradient in the north-south direction. The Figure 5, an overlay is shown at the locations of the largest gradients calculated for the areas. On the raster image of the slopes, a high-pass filter was applied to highlight slopes greater than 35° and 40°.

4. Site Investigations

The area was investigated with four boreholes. The S1 (depth of 60.0 m) more upstream in the Nibbiano area, the S3 (depth of 51.2 m) and the S4 (depth of 60.0 m) more downstream in the Sant’Erasmo area, and the S2 (depth of 60.0 m) in the intermediate part. In addition, refraction seismic surveys were carried out. The L1 line between Nibbiano and Sant’Erasmo hamlet and the L2, L3, and L4 in the transversal direction. SH waves were also recorded to estimate the vs. velocities in the L2 line. Finally, an electrical tomography (along the L1 line) was carried out to improve the accuracy of the interpretation, and n.12 ambient noise stations were recorded (HVSR). Figure 6 shows the location of the investigations. The position of each test was taken from a preventive study of the major criticalities found by the surface survey in the area.

4.1. Geotechnical Characterization of the Materials

The stratigraphy shows a considerable thickness of detrital materials in the most superficial part of the boreholes. The granular nature, due to the component of calcareous debris, did not allow the taking of undisturbed soil samples. Therefore, the SPTs were performed to characterize the mechanical behaviour of the superficial layer.

The processing was implemented using the characteristic value concept. Eurocode 7 defines the characteristic value as a precautionary estimate of the parameter that influences the start of the limit state. The estimates were based on the Student’s t distribution (more suitable for geotechnical problems) using the following formula.

$$X_k=\underset{\_}{X}\pm t_{n-1}^{0.95}\left[{\displaystyle \frac{s}{\sqrt{n-1}}}\right]$$

(4)

where $X_k$ is the characteristic value of the parameter; $\underset{\_}{X}$ is the average value; $t$ is the value of the “Student” distribution with n − 1 degrees of freedom; $s$ is the standard deviation of the sample; $n$ is the number of data used.

Table 1. Geotechnical debris properties.

Parameters	Records	Min	Max	Average	DevST	Characteristic
Friction (°)	9	24.49	48.54	34.29	±7.38992	29.43
Young Mod. (MPa)	9	3.60	24.30	10.33	±6.73515	5.91

shows the mechanical parameters only for the superficial detrital material.

The bedrock involves marly limestones with a medium degree of fracturing. For each meter of perforation in the rock, the quality of the cluster was estimated using the RQD value. The Rock-Quality-Designation is the ratio (in percentage) of the sum of the healthy core segments with length greater than 0.1 m and the total length of the section in which the estimate is performed. The method [68,69,70,71] was used to define the quality and degree of fracturing. Figure 7 shows the individual values for each survey and their range of distribution.

The values show a homogeneous bedrock between Nibbiano and Sant’Erasmo hamlet. The RQD values vary between 0.0 and 95.8%, with an average value of 56.4% and a standard deviation of 24.64. The data (Figure 8) shows that at the top of the rock mass, RQD values are lower than 50%. This aspect is more marked in the S1, S2, and S3 boreholes for the first 3/5 m, while at greater depths, the quality of the cluster appears superior. Similar conditions for the S4 borehole, where the poor thickness is reduced by about 2 m. This aspect is due to the presence of discrete fracturing in the highest part of the rock, which gives the mass moderate resistance and poor stiffness.

For the geomechanical characterization, the failure criterion developed by [72] and updated by [73] was used. The method defines the resistance in terms of major and minor principal stresses, starting from the characteristics of the intact rock and deriving the representative properties of the mass based on the characteristics of the joints. Among the main classification systems of rock masses, the Geological Strength Index value is more representative. In the evaluation of GSI, there are various interpolation matrices for various rock types [74] or criteria for the evaluation of the description of the structure and conditions of the discontinuities. In this work, it was considered more appropriate to use the equation proposed by [75], which indicates calculating GSI as a function of the RQD parameter as:

$$GSI=18.7 \mathrm{*} exp exp\left(0.0125 \mathrm{*} RQD\right)=44$$

(5)

Furthermore, the uniaxial compressive strength (Uniaxial Compressive Strength) was calculated. The ${\sigma }_{ci}$ is defined as the maximum stress supported before failure. In this case, it was estimated between 10 ÷ 11 MPa derived from the results of the uniaxial compression tests performed on intact rock samples. Finally, a disturbance degree of the rock mass ranging from 0.7 (Good) to 1.0 (Poor) was considered. Since the uniaxial compressive strength ${\sigma }_{ci}$ is less than 100 Mpa, the deformation modulus of the rock was calculated by [73] according to the following formula:

$$E_m(\mathrm{G}\mathrm{P}\mathrm{a})=\left(1-{\displaystyle \frac{D}{2}}\right)\sqrt{{\displaystyle \frac{{\sigma }_{ci}}{100}}}·{10}^{({\displaystyle \frac{GSI-10}{40}})}$$

(6)

Considering a variable degree of disturbance, a range of elastic modulus between 0.18 and 0.83 GPa was obtained. In the analysis, the bedrock was considered as an equivalent continuum model. For this, it was necessary to determine the failure criterion in terms of Mohr-Coulomb and define a stress range, a friction angle, and an equivalent cohesion for the rock mass. The approximation of the non-linear Hoek-Brown failure envelope with the linear Mohr-Coulomb failure envelope was done by fitting the curve generated by the following equation for a range of minor principal stresses defined by ${\sigma }_t<{\sigma }_3<{\sigma }_{3max}$:

$${\sigma \prime }_1={\sigma \prime }_3+{\sigma }_{ci}{\left(m_b{\displaystyle \frac{{\sigma \prime }_3}{{\sigma }_{ci}}}+s\right)}^a\to \tau =c\prime +\sigma ·tan\phi \prime$$

(7)

where, as stated by [73], $m_b$, $s$ and $a$ are constants of the material. The fitting process involves balancing the areas above and below the Mohr–Coulomb diagram. The parameter values for the bedrock are shown in

Table 2. Geomechanical rock properties.

Hoek-Brown Classification					Hoek-Brown Criterion			Mohr-Coulomb Fit
sigci	GSI	mi	D	Ei	$m_b$	$s$	$a$	$c\prime$	$\phi \prime$
MPa				MPa				MPa	°
10.5	44	13	0.7	48.0	0.599	0.0003	0.509	0.174	30.72

where sigci is the intact uniaxial compressive strength value; GSI is the geological strength index; mi is a material constant for the intact rock; D is the disturbance factor; Ei is the elastic modulus of the intact rock; $m_b$ is the reduced value of the material constant mi; $s$ and $a$ are material constants which depend upon the characteristics of the rock mass; $c\prime$ is the cohesion and $\phi \prime$ is the friction angle.

.

4.2. Surface Geophysical Surveys

To obtain a continuous vision of the subsoil to more depth, a serious of surface geophysical surveys were carried out. Seismic refraction surveys [76,77] were carried out with the roll-along technique, according to acquisition sequences with overlapping of 24 geophones at a time.

Table 3. Seismic acquisition parameters used in geophysical surveys.

Parameters.	Line 1	Line 2	Line 3	Line 4
Length (m)	835	355	360	350
Type	Wave P	Wave P-SH	Wave P	Wave P
Intergeophone distance (m)	5.0	5.0	5.0	5.0
Number of geophones	168	72	72	72

shows the main acquisition characteristics.

The seismic recordings were analysed with tomographic techniques by progressively reducing the convergence error in an iterative way to define the best distribution of velocities by comparing the real arrival times with the theoretical ones. The method allows us to calculate, for each source-receiver pair, the optimal path of the seismic rays (ray-tracing) and to derive the theoretical dromochrones. The process has allowed us to define a two-dimensional model of the soil with a detailed reconstruction of the buried morphometry and any discontinuities. In fact, seismic waves are correlated with stiffness and, therefore, represent a fundamental parameter for the definition of bedrock and the stratigraphic relationships between it and the sediments covering it [78].

Figure 9 shows the 2D model. The length profile (L1 line) is the one that shows the stratigraphic in the most detailed and complete way. The superficial part has a Vp variability range between 500 and 2000 m/s. There is an improvement of the stiffness until reaching the deep bedrock with P-wave velocities between 2200 and 2850 m/s. Even on the L2, L3 and L4 lines, the two-layer interpretation was possible.

A geo-electrical investigation [79] was performed to improve the interpretation of the cover/bedrock interface. An ERT profile (Electrical Resistivity Tomography) was recorded along the L1 Line, and the main acquisition characteristics are shown in

Table 4. Geoelectric line acquisition parameters.

Parameters	Line 1
Length (m)	720
Type	Array Wenner-Schlumberger
Interelectrode distance (m)	10
Number of electrodes	72

.

Geoelectrical surveys are widely used methods for the reconstruction of stratigraphic buried soils [80,81]. They are also used for the study of landslides because they provide excellent results both in the horizontal and vertical sense. The geoelectrical line was carried out with a 720 m extension and 10 m electrode-electrode spacing, which allowed it to reach considerable depths with good data reliability (Wenner-Schlumberger configuration). The profile shows a superficial horizon of about 20.0 m characterized by resistivity values > 40 Ohm.m and a transition to more conductive materials at depth with resistivity values < 40 Ohm.m. Figure 10 shows a reconstruction of the lithological contacts and of the electrical resistivity anomaly identified shortly before the Sant’Erasmo hamlet.

Before the geological and geophysical investigations, some passive recordings of environmental noise were carried out using the HVSR technique [82,83,84]. Studies in the literature highlight the usefulness of the HVSR test as a support for defining the volumes of potentially unstable masses [85,86,87]. In this work, the HVSR surveys have had a dual purpose. The first is to provide preliminary indications of the depth of contact characterized by greater contrast in soil stiffness, allowing the preparation of an optimized geognostic survey for typology and positioning. Secondly, the HVSR investigations were allowed to evaluate the reliability of transferring seismic stratigraphic models from areas well characterized by geognostic investigations to uninvestigated areas based on similar HVSR curves [88,89]. Twelve HVSR tests were carried out between the Nibbiano and Sant’Erasmo hamlet using a triaxial surface sensor with a natural frequency of 2 Hz. The guidelines indicated in the InterPACIFIC project [38] were followed for positioning, acquisition, and data processing. The tests were characterized by a signal acquisition time of 20 min with a sampling frequency equal to 250 Hz. For all the analyses, a window length of 20 s was used, and good signal coverage was guaranteed in the windows selection process (signal coverage min = 58%—mean = 85%—max = 93%). The post-processing was carried out with the same parameters for all tests; ‘Tapering’ (Enabled—bandwidth = 10%) and ‘Smoothing’ (Konno-Ohmachi—bandwidth = 40) [90]. Figure 11 shows the clustered results of the environmental noise study.

5. Numerical Analysis

5.1. Seismic Analysis of the Slope

The seismic stress, due to the cyclical actions, generates an increase in destabilizing actions in the slopes due to the generation of inertial forces proportional to the seismic acceleration. At the same time, in the soil, a reduction of the resistant actions is generated due to fatigue phenomena. In general, the reduction of the shear strength increases with the load cycle number and, therefore, with the duration of the seismic event. Furthermore, these behaviours are closely linked to the seismo-mechanical characteristics of the soils. The criterion used in this work to estimate the seismic conditions of the slopes is the pseudo-static method. The margin of safety was evaluated with respect to the conditions of limit equilibrium or incipient collapse; therefore, it was evaluated with respect to an ultimate limit state. Since the effects of the inertial forces produced by the seismic motion are constant over time, in modulus and direction, the dynamic effects of the seismic load were represented by an equivalent static action. The analysis consists of introducing a seismic coefficient proportion to the maximum seismic acceleration expected at the site based on the characteristics of the reference seismic event. The selection of an appropriate coefficient is the most important and difficult aspect of a pseudo-static stability analysis. The seismic coefficient controls the force applied to the mass of each individual element, so its value must be related to the measure of the magnitude of the induced inertial force. Typically, the seismic coefficient value is a function of the peak horizontal ground acceleration $a_{max}$ [91]. In this work, the most appropriate seismic coefficient value was estimated by the 2D local seismic response analysis. The seismic amplification is the result of multiple physical phenomena (multiple reflections, diffraction, focusing, resonances, etc.) that waves undergo in correspondence with the heterogeneities and discontinuities of the most superficial layers and topographic irregularities. Local seismic response studies [92] allow us to evaluate the actual modifications that the seismic signal undergoes in its path from the seismic bedrock to the topographic surface [93,94]. The use of subsoil categories relating only to the vs. parameter does not consider the seismic stiffness contrast defined by the ratio between the values of shear wave propagation velocity. Therefore, in this work, the seismic amplification phenomena were improved by the 2D numerical analyses capable of modelling the topographic and stratigraphic conditions of the area.

The analysis was divided into various work phases, one preparatory to the other:

• Estimation of the base acceleration with probabilistic approach (PSHA), definition of the disaggregation data, and identification of magnitude-distance pairs.
• Selection of the natural input accelerograms.
• Reconstruction of the stratigraphic model and seismic parameterization.
• 2D local seismic response analysis in Linear Equivalent (LEQ) conditions.
• Estimation of the maximum accelerations and calculation of the horizontal seismic coefficients.

From the seismic hazard map provided by INGV for the Italian territory, the basic seismic parameters (ag, Fo, T*c) were defined, and the disaggregation of the seismic hazard was estimated by evaluating the relative contributions of different seismic sources for different magnitude-distance pairs. The definition of these parameters allowed the selection of the spectral accelerograms compatible with the target spectrum (Cat. A—T1) of the area relating to the 475-year return period. Figure 12 shows the waveforms of the earthquakes used as input in the analyses.

For each accelerogram, an analysis was performed, obtaining the final PGA (Peak Ground Acceleration) of the area as the average of the individual analysis. The simulations carried out with an equivalent linear approach consider the variation of the mechanical characteristics of the soils in terms of stiffness and damping as a function of the stress condition induced by the seismic force. Since the 1980s, numerous experimental approaches have been conducted for the analysis of the dynamic properties of granular materials to define the variation of the shear modulus G and the damping ratio D with the shear deformation γ. Among the most important and interesting tests is the study by Rollins et al. [95], who, starting from the hyperbolic relation proposed by Hardin and Drnevich [96] for the sands, defined the mean curves for the evaluation of G/Gmax—γ e D—γ for the gravels. Rollins defines these curves as independent of both the sample disturbance and the fine content. For the present work, given the granular nature of the cover soils, the [95] relations were used to consider the non-linearity behaviour of the materials under cyclic earthquake actions.

The bedrock, instead, was considered to have a perfectly elastic behaviour.

Table 5. Stress function G/Gmax e Damping used in the analysis.

Layer	G/G0 and D Function
Cover layer	Rollins et al. [95]
Bedrock	Elastic behaviour

shows the variation curves used for each layer.

The 2D numerical model was assembled using the seismic line L1 results between the Nibbiano and Sant’Erasmo hamlet. The seismic amplifications were analysed with the Quad4M code [97] implemented with VisualQ4M [98]. The VisualQ4M is a pre/post-processor capable of discretizing the geo-lithologic model with a 2D finite element mesh and reordering the numbering of the nodes to reduce matrix bandwidth and, consequently, the calculation times. An important phase of the assembling of the model is the discretization. Excessively coarse elements can produce filter effects for the high-frequency components that were avoided since the small wavelengths cannot be adequately modelled by nodes that are too distant from each other. To eliminate the problem, it was imposed an element height $h_{max}$ equal to:

$$h_{max}\le \left({\displaystyle \frac{1}{5}}÷{\displaystyle \frac{1}{8}}\right){\displaystyle \frac{V_s}{f_{max}}}$$

(8)

where $V_s$ is the vs. wave velocity, and $f_{max}$ is the maximum frequency considered in the analysis equal to 20 Hz.

The geometric model was built by discretizing only the cover material, while the bedrock was considered as an infinite half-space with a transmissible boundary, therefore able to absorb part of the incident seismic waves. Each element was assembled in a concentrated mass scheme using springs and viscous dampers as connections. The iterative resolution, with direct integration in the time domain according to the Newmark scheme, minimizes the Rayleigh coefficients based on the fundamental frequency of the model. Figure 13 shows the average distribution of the calculated maximum horizontal accelerations. There are two main focuses (max PGA), one near the Nibbiano hamlet and the other a little further downstream.

5.2. Three-Dimensional Slope Stability

The study of slope stability examines the relationships between the mechanical characteristics of the ground (which oppose each other) and the external forces (agents) that destabilize it. It is fundamental to predict the degree of safety associated with the slope along which significant slides could occur, as these phenomena have important consequences, both from a human and economic point of view [99]. For complex geometries, the use of finite element approaches is more appropriate as it does not impose the a priori hypothesis of the position and/or shape of the sliding surface. In this work, slope stability was studied with the FLAC3D v 7.00 software, an explicit finite difference code for continuum mechanics calculations.

The code simulates the behaviour of continuous systems, which undergo plastic flow when the yield limits of the materials themselves are reached. Given the extension of the area between the Nibbiano and Sant’Erasmo hamlet, a model of 800 × 800 m was built. To obtain consistent results, a mesh was realized with elements of variable dimensions from 10 m in the most superficial part to 30 m in the deepest part.

In the model, the most superficial part has been schematized with a variable thickness (grey colour in Figure 14) function for both the geognostic investigations and the geophysical interpretations.

The cover material was modelled according to the Mohr-Coulomb theory using an effective cohesion equal to 5.0 kPa (conservative value) and an internal friction angle equal to 29.5°. The bedrock, instead, was simulated according to an equivalent continuous model, applying a stress-strain approach according to the Ubiquitous-Joint Plasticity criterion. This approach considers some planes of weakness within a Mohr-Coulomb stress-strain relationship. In the specific case, an effective cohesion equal to 171 kPa and an internal friction angle equal to 32.7° was used, while for the planes of weakness, an effective cohesion of zero and an internal friction angle equal to 16° was used.

The analysis was performed to evaluate the general deformation state and the possible evolution of the system according to two hypothetical levels of the piezometric level and different states of dynamic load. In this paper, we show the results for:

• Static analysis with a water table equal to −1.0 m from ground level.
• Pseudo-static analysis with a water table equal to −4.0 m from ground level.

The most widespread method for dynamic slope stability is certainly the pseudo-static one. In this approach, the seismic action is assimilated to an equivalent static force of an entity equal to the product of the seismic coefficient $k$ and the weight of the potentially unstable ground. To obtain an analysis consistent with the real conditions of the slope, the value of the seismic coefficient was correlated to the seismic performance of the slope under the cycling actions. The seismic coefficient was expressed according to:

$$k_h={\beta }_s·{\displaystyle \frac{a_{max}}{g}}=0.112$$

(9)

where $a_{max}=0.4 g$ is the maximum horizontal acceleration resulting from the analysis of the local seismic response 2D, $g$ is the acceleration of gravity, and ${\beta }_s=0.28$ is the reduction coefficient of the maximum acceleration at the site by Italian building code. In the analyses, the seismic coefficient, considering the horizontal component, is equal for both X and Y directions according to the combination reported in

Table 6. Planar combination (ref. Figure 15) of the seismic coefficients.

Direction	X	Y
Seismic load	$+K_h$	$+K_h$
	$+K_h$	$-K_h$

.

The slope stability analyses were performed using the Strength Reduction Factor procedure. The method is based on the repetition of the same stability test in the same general conditions but introduces systematically and proportionally reduced values in the resistance parameters only. Numerically, collapse occurs when it is no longer possible to obtain a convergent solution of the stress-strain relationship and, therefore, of the global equilibrium. Upon reaching the generalized collapse, the global safety factor of the slope was calculated, assuming the last value of SRF that verifies the stability condition of the system. The Strength Reduction Factor (SRF) is defined as:

$$SRF=\left({\displaystyle \frac{tan{\phi }^{\prime }}{{tan{\phi }^{\prime }}_f}}\right)=\left({\displaystyle \frac{c^{\prime }}{{c\prime }_f}}\right)$$

(10)

where ${{\phi }^{\prime }}_f$ and ${c\prime }_f$ are the frictional and cohesive parameters of breaking resistance. Dawson [100] concluded that the value of the elasticity parameters, therefore Young’s modulus (E) and Poisson’s ratio (ν), have little influence on the result of the safety factor; therefore, the effects of the elastic parameters, in this work, were not a consideration. Figure 15 shows the results obtained by the static calculation and with the piezometric level equal to 1.0 m from the elevation of the ground plane.

Figure 16a shows the maximum values of the shear deformations. The values were concentrated in the southern/south-eastern portion of the Nibbiano hamlet. Typically, the application of the SRF method produces a single safety factor that corresponds to a global minimum stability state. However, for this work, a study of multiple minimum states is more interesting to generate a distribution of the safety factors within the model itself. Figure 16b shows the colour-map trend of the safety factors distribution.

The water table was gradually lowered to −4.0 m depth from ground level, and the pseudo-static load was applied with a force proportional to $K_h$. Also, in this case, the results identify a condition of global instability with safety factors equal to one. Figure 16 shows the trend of the maximum shear deformations.

6. Discussion

This study allowed us to identify some geological particularity that was previously unknown. Indeed, in the area, unpublished geological investigations were conducted to reconstruct a real geometric 3D model and to understand the dynamic soil behaviour with the same conditions as the 2016–2017 seismic sequence. Furthermore, the results were used to simulate by the 3D numerical modelling the possible evolution of the slopes.

According to the geomorphological study, the Nibbiano hamlet is the one most affected by steep slopes (>40°) recognized as scraps of old paleo-landslides. These shapes appear elongated on the north-south alignment with typical arched shapes. The landslide extends along the eastern slope of Monte Campalto, completely encompassing the Nibbiano hamlet. It is possible to state that the entire landslide is separated in two parts by the geometric rise on which both Nibbiano and Sant’Erasmo hamlets stand. The northeastern part of the area presents lower energy, with slopes about 10° maximum; the southwestern portion, instead, presents slopes with greater energy, various counter-slopes, and multiple scarps (Figure 17). In general, the low stiffness of the debris covers, and their fair degree of permeability both contribute to generating possible scenarios of instability due both to intense and/or exceptional rainfall and to induced effects by strong earthquakes. The Nibbiano Hamlet shows both scenarios. These areas present conditions of instability generated along the slopes with significant sliding depths. These instabilities are generated by both the increase of the piezometric level, already considerable given the presence of water emergencies found in the territory, and by possible seismic accelerations such as those already occurred in the Central Italy crisis of 2016–2017. The Sant’Erasmo hamlet, instead, has no relevant morphometric shapes because it is located on a very low-slope plateau with no movement energy. In this area, the damages are due only to recent seismic activity. Furthermore, it is possible to exclude the presence of deep gravitational deformations between the two hamlets.

By the borehole stratigraphies, it was possible to identify a more superficial part with granular behaviour due to the abundant component of calcareous debris, as well as a deep bedrock with altered marly levels in the clay matrix. The superficial part has variable thicknesses from 10.0 m in the S4 borehole to 19.5/20.0 m in the S1, S2, and S3 boreholes. The granular nature of the superficial deposits derives from the accumulation of debris along the eastern slope of Monte Campalto due to both the action of gravity and the effect of washing. The drilling of the bedrock was carried out using a double-core barrel with a thin-walled crown and, therefore, with a minimum cutting surface. This method allowed zero disturbance to the material and a good vision of the extracted core. The bedrock shows some fractures with inclinations varying from 45° to 60° (Figure 18a–c) with more frequency at the top of the cluster. These slickensides are tectonic lineaments clearly visible on the cores and evident signs of fault activity. Instead, in the deeper portions, the mass rock appears completely intact (Figure 18d).

By the interpretation of the seismic profiles together with the stratigraphy of the boreholes, a bi-layered model was defined with a gradual increase of the Vp wave velocity in depth from 250 m/s to 3000 m/s. The superficial part can be referred to as sediments composed mainly of calcareous breccias in a silty-clayey matrix, while the highest velocities were recorded in the deep bedrock consisting of marly limestones. The seismic profiles show a good geometrical correlation, too; in fact, the top of the bedrock was recorded at a depth of 20.0/25.0 m from the ground level with a non-linear trend. The irregularities of the debris/bedrock contact are due to a portion of about 4/5 m of rock with a high degree of fracturing, which is probably linked to the tectonic nature of the places and to the thrust planes. The section L1-L1′, between the Nibbiano and the Sant’Erasmo hamlet, was interpreted in Figure 19 and allows us to define the geometry of thrust planes in the slope. It is evident that they are not arranged according to the “natural” 30–40° inclination but have smaller slopes. These can be interpreted as involvement in the gravitational deformation of the slope, which would tend to rotate the loose materials and some portions of the bedrock together with the planes during the movement towards the valley.

The L4 seismic profile is very interesting in particular. The results highlight a depression with a pronounced concavity towards the top in correspondence with the Sant’Erasmo hamlet. This geometry was interpreted as a small, buried valley that, in seismic conditions, generates amplification and focusing of seismic waves on the surface. In fact, following the earthquake between August and October 2016 and January 2017, the Sant’Erasmo hamlet reported damages to all existing buildings.

The geoelectric profile (ERT) L1 is according to the seismic interpretations and the stratigraphic results of the boreholes. In fact, the profile shows a superficial layer with a thickness of about 20.0 m and resistivity values >40 Ohm.m. This layer was interpreted as remobilized and aerated material attributable to the cover detrital deposits. The results show, in-depth, a transition to more conductive materials typical of the calcareous marls with resistivity values and more homogeneity.

The results of the HVSR tests allowed us to identify the max frequency F0 to define the depth of contact between the cover layer and the bedrock. The tests 5, 7, 11, and 12 show very high fundamental F0 frequencies. These can be interpreted as reduced thicknesses of the cover layer, and they can be considered in some cases as the limit of closure of the volumes. All the HVSR tests with lower frequencies are aligned along the L1 line. These interpretations were essential to define the survey program and for resource optimization.

The results of the geognostic and geophysical interpretation show possible scenarios of instability due both to intense and/or exceptional rainfall and to induced effects by strong earthquakes. Numerical stability analyses were performed to define the stability conditions and explain the possible evolutions of the slope. The data from the investigations were used to reconstruct the 3D model, which was subsequently analysed using Flac code (Itasca Inc., Minneapolis, MN, USA). Particular attention was given to the geometric model. The mesh was realized with tetrahedron elements with variable side dimensions from 10 m in the most superficial part to 30 m in the deepest part. This aspect was important to reduce the calculation time and, at the same time, obtain a good resolution in the upper part. The analyses were elaborated with the strength reduction procedure to obtain an estimate of the safety factor for each condition. Usually, the strength reduction method produces only one factor of safety corresponding to one global minimum stability state. In this work, multiple minimum states were calculated because this is the most interest along the complex slope profile. Therefore, a safety map was constructed for both static and pseudo-static analysis for each model analysed. The slope’s seismic performances were evaluated by the numerical modelling of a 2D section between Nibbiano and San’Erasmo hamlet. The results show a focus under the Nibbiano hamlet and in the adjacent areas due to the multiple reflections and focusing that waves undergo in correspondence to the topographic irregularities and to the bedrock trend. Moreover, the maximum horizontal acceleration value equal to 0.4 g was calculated. An important aspect of this work is the selection of seismic forces consistent with the real conditions of the slope for performing a 3D stability analysis in pseudo-static conditions. To obtain an appropriate seismic coefficient, the value was correlated to the seismic performance of the slope under the cycling actions. Therefore, the seismic coefficient $k_h$ equal to 0.112 was calculated. The results show that the slope area is very vulnerable to pore pressure variations. In the static analysis, using a piezometric level of −1.0 m from ground level, the safety factors are equal to the unit corresponding to an instability condition. The contours <= 1.0 m in the safety map include all built of the Nibbiano hamlet and other portions to the further southern model but exclude the entire Sant’Erasmo hamlet. The results of the pseudo-static analysis are very similar. In the model, a piezometric level of −4.0 m from ground level was used. Even in this case, the Nibbiano Hamlet is entirely involved in the contours of instability, unlike the Sant’Erasmo Hamlet.

7. Conclusions

This work highlighted how essential it is to integrate different approaches to define the 3D geological model. The data were used to define the geometry of what is thought to be the volume of soil potentially affected by gravitational movements. Specifically, it was possible to:

• check the information on the areas with detailed studies aimed at confirming the current geological, geomorphological, and structural knowledge;
• define the buried geometries for an areal extension capable of involving the two hamlets;
• estimate of the geometric variations of the debris covers and of the bedrock in depth;
• estimate the mechanical characterization of the debris along the slope and of the bedrock;
• estimate, by the 2D local seismic response, of the maximum seismic acceleration for the two hamlets;
• identify, by the 3D numerical modelling, the main deformations and instabilities distinct for the two hamlets;
• define a real deformation model that can be extended to neighbouring areas with similar lithological-structural, lithotechnical, and geomorphological characteristics, which is useful for the creation of monitoring systems and possible remediation.

The use of geognostic and geophysical investigations demonstrates once again how a multi-source approach cannot be ignored for the definition of complex geo-morpho-stratigraphic models. The stratigraphic configuration resulting from the area and the very nature of the materials meant that it was possible to highlight a good convergence of the results between the direct, geotechnical tests (boreholes, SPT, laboratory analysis) and indirect, geophysical tests (Seismic, Geoelectric, HVSR). Environmental noise records processed with the HVSR technique also proved to be extremely efficient, providing indispensable support to the definition of the investigation program. The data obtained allowed the reconstruction of a three-dimensional detailed model of the cover/bedrock interaction between the Nibbiano and Sant’Erasmo hamlet to interpret the local geology of the slope and to understand the ongoing landslide process. Finally, 3D numerical analyses were useful in defining the extension of the unstable areas, the maximum possible deformations, and the maximum depth of sliding surfaces. Furthermore, the results were used to understand the future evolution of the slopes.

This research has shown that the use of adequate investigative techniques is fundamental to studying complex landslides. Indeed, an in-depth knowledge of both the intrinsic conditions of the slopes and possible evolutionary scenarios can ensure greater attention in the management of dangers for structures and especially for the lives of the people.