Post-Collapse Evolution of a Rapid Landslide from Sequential Analysis with FE and SPH-Based Models

Abstract

Propagation models can study the runout and deposit of potential flow-like landslides only if a reliable estimate of the shape and size of the volumes involved in the phenomenon is available. This aspect becomes critical when a collapse has not yet occurred and the estimation of the unstable volume is not uniquely predictable. This work proposes a strategy to overcome this problem, using two established analysis methods in sequence; first, a Strength Reduction Method (SRM)-based 3D FEM allows the estimate of the instable volume; then, this data becomes an input for a Smoothed Particle Hydrodynamics (SPH)-based model. This strategy is applied to predict the possible evolution of Sant’Andrea landslide (North-Eastern Italian Alps). Such a complex landslide, which affects anhydrite–gypsum rocks and is strongly subject to rainfall triggering, can be considered as a prototype for the use of this procedure. In this case, the FEM–SRM model is adopted, which calibrates using mapping, monitoring, geophysical and geotechnical data to estimate the volume involved in the potential detachment. This volume is subsequently used as the input of the SPH model. In this second phase, a sensitivity analysis is also performed to complete the evaluation of the most reliable final soil deposits. The performed analyses allow a satisfactory prediction of the post-collapse landslide evolution, delivering a reliable estimate of the volumes involved in the collapse and a reliable forecast of the landslide runout.

1. Introduction

In the recent literature, the propagation analysis of real rapid landslides has been utilized considerably [1,2,3,4]. Over the last 20 years, a large number of particle-based or mesh-less methods have been under development using different numerical strategies [5,6] to overcome the limit imposed by the small strain assumption typical of finite element modeling. These new methods still implement the continuity equations, but adopt a discretization constituted by a cloud of material points or particles, which do not have an explicitly defined connectivity. All the physical properties are attached to the particles and not to the mesh. Among them are the Smoothed Particle Hydrodynamics (SPH) methods [7,8,9,10,11,12], the Material Point Method (MPM) [13,14], and the Particle Finite Element Method, (PFEM) [15,16]. Generally, the choice of the rheology to use and the consequent calibration of the soil parameters becomes a tricky and essential phase of the analysis. Various authors have applied different approaches to calibrate the model, thus obtaining the most reliable rheological parameters to be assigned [6,17,18,19,20]. In almost all these cases, the user has to assign the basal topography of the site, which remains unchanged during the simulation; then, the detachment volume has to be defined to clearly distinguish the deforming part of the slope. Clearly, the unstable volume is known only if the failure has already occurred, but, when the model is implemented for a prediction, the assumption about the starting geometry is affected by a great uncertainty degree. Based on the authors’ knowledge, when a propagation model is used for prediction, the authors focus more on the soil model calibration [17,18,19,20,21,22] taking for granted the geometry of the unstable volume. However, an error at this stage can seriously affect the estimation of the area finally involved in the run-out.

In this regard, the present study aims at combining two established analysis methods in sequence: the first, performed at small strains, to define the slope stability conditions and to forecast the unstable volume; the second, at large deformations, estimating the displacement evolution after the trigger, starting from the results of the first phase. The combination of two analyses allows the determination of the most reliable unstable volumes and their possible kinematic behavior after collapse, also considering different hydraulic conditions in a slope characterized by a very complex geological setting.

The first phase of the procedure is based on a FEM-based stability analyses, performed through the shear strength reduction technique [23,24,25,26]. It provides the safety margin (also known as factor of safety or FS) of the studied slope as the number by which the original shear strength parameters must be divided in order to bring the slope to failure. The failure is assumed as the condition in which the analysis is not still able to converge to a numerical solution. At the end of the analysis, among other results, a FEM stability analysis allows the estimation of unstable volume, i.e., the volume that would be affected by significant displacements when, according to the reduction technique, the slope reaches instability conditions.

FEMs are mesh-based methods that traditionally solve the static equations according to the basic hypothesis of “small strain”. When the shear strength reduction is set, the soil locally reaches a state of plasticization, i.e., it accumulates large plastic deformation without a unique stress-strain relationship. The errors introduced in the numerical solution at this state do not permit a correct evaluate of the displacement entity, that become very large. Consequently, considerations about the evolution and propagation after the slope collapse are not possible.

The GeoFlow-SPH model proposed by Pastor et al. [11], specifically developed for studying the flow-like landslide, is adopted to perform the propagation analysis. It is an interesting compromise between a sufficiently good description of the propagation processes on a 3D real morphology and a simplified description of the flow behavior of the propagating mass. Its theoretical formulation is based on the hypothesis that the propagating wave has a height smaller than its longitudinal or transversal sizes. This assumption makes it possible to not consider one dimension of the problem, because velocities and displacements perpendicular to the base are negligible, and to adopt a depth-average technique for the evaluation of many flow characteristics, first of all the mean velocity. As in similar models, the soil is idealized as an equivalent homogeneous one-phase material; the rheological behavior of the water–solid mixture is described with a simple viscous law relating the displacement rate to the shear stress [5].

This procedure is applied to the case study of the Sant’Andrea landslide, a large slope instability that has been active in the North-Eastern Italian Alps for many years which has been subject to continuous topographic monitoring since 2013. The slope presents a particularly complex geological condition, involving anhydrite–gypsum rocks affected by karst groundwater circulations. Thanks to many in situ investigations carried out so far, it is possible to reconstruct a 3D terrain profile on which to perform the stability analysis with a 3D-FEM. In this analysis, indications from the surficial displacement monitoring play a crucial role in improving the in-depth stratigraphy assessment, starting from soil parameters chosen a priori. The unstable volume determined with 3D-FEM becomes the input data in a SPH propagation model [11]. In this second phase, as it is not possible to evaluate with specific tests the rheological characteristics of the involved soils [6], a sensitivity analysis is performed to have an evaluation of the most reliable deposits formed in the valley downhill in the case of a paroxysmal collapse event.

2. Method: FEM and SPH Models

The FEM analysis is performed using the Midas GTS NX 2019 v2.1 package (MIDAS IT Co., Ltd., Seongnam, Korea). The program allows insertion of the three-dimensional topography of a slope starting from the contour lines of the site. The soil stratigraphy can be reproduced by inserting the data collected in some vertical boreholes and with the surficial geological survey, which are automatically interpolated by the software to obtain the interlayer surfaces. The 3D hydraulic conditions can be assigned imposing a constant pressure in one or more confined aquifers or assigning a defined piezometric surface and performing a previous seepage analysis to evaluate the total head distribution. Several linear or non-linear constitutive models may be assigned to the materials. In order to evaluate the stability conditions of a slope, after the setting up of the model, a non-linear static analysis is firstly performed to provide the at-rest effective stress distribution in relation to the given hydraulic conditions and the soil self-weight. Afterwards, a shear strength reduction procedure is applied (in Midas, called the Strength Reduction Method (SRM)), which gradually decreases the resistance parameters of all the soils by a certain Strength Reduction Factor (SRF). The code continues to increase SRF until there is no numerical convergence to a solution. This situation is considered critical because the number of nodes that exceed their shear resistance (plastic points) is so high as not to ensure the slope stability. The maximum SRF before reaching non-convergence is assumed as the slope FS.

Among the various results that can be analyzed with a FEM analysis thus carried out (distributions of stresses, strains, pore pressure, displacements, etc.), basically two outcomes are extracted for the aim of the present study: the distribution of the surficial displacements and the 3D sliding surface, i.e., the inner surface on which the maximum shear strain is reached before the non-convergence of the SRM procedure. The first information is extremely important for evaluating the reliability of the model geometry. In fact, the orientation and distribution of surficial displacements depend on the position and the shape of the sliding surface, which in turn are controlled by the soil stratigraphy. Consequently, the first attempt soil stratigraphy, automatically determined by the software in the input phase, is slightly corrected to better reproduce the displacement distribution detected with the topographic monitoring. Any correction is imposed always respecting the outcomes of both the surficial geomorphological survey and the geognostic boreholes.

The 3D sliding surface is detected by considering the iso-displacement surfaces resulting from the FE analysis. An iso-displacement surface is a 3D surface that includes all the points characterized by a certain total displacement value. If a threshold value for the total displacement is chosen, the corresponding iso-displacement surface divides the slope in two portions, one stable because it shows displacements below the threshold, and another unstable because it exhibits larger displacements. It is important to note that the threshold selection is arbitrary and it is not possible to define a single absolute value, because models reproducing the same geometry but with different hydraulic conditions exhibit different FS_max and different final ranges of displacements. It appears more correct to express the threshold as a percentage of the maximum displacement reached in the SRM analysis.

After identification of the critical surface, it is possible to obtain the slope portion potentially unstable and its base. These two outputs are transferred as input in the propagation analysis, performed with the GeoFlow-SPH model [11]. This model integrates the Zienkiewicz–Biot mass and momentum balance equations of a continuum media, developed with the Lagrangian formulation. As the SPH method is a particle-based model, the continuum media is described as an assembly of virtual particles which transport all the information about the physical quantities (mainly density, stress, strain, flow rate). The interactions among particles are described through the kernel functions, which are gaussian functions describing the probability of a particle to have contact and momentum exchange with the other neighbor particles, in relation with their spacing. The model adopts a good compromise to model the flow-like landslide propagation, assuming that the flowing mass is distributed along a longer (or wider) base with respect to its depth. A depth-averaging in the z direction is adopted for the evaluation of physical quantities such as the velocity, the density, etc. Thus, one dimension of the flowing mass is lost because the velocity and the displacement perpendicular to the base are not considered. If this type of approach greatly reduces the calculation times, it must however be considered that the simplification of the model can lead to less accurate results, especially when the heights of the moving material increase. However, the wide use of the model in literature allows the results to be considered reliable.

The heterogeneous material that constitutes the moving mass is considered as an equivalent fluid [20], governed by simple rheological relationships, with constant physical and mechanical characteristics. Several models, such as the Mohr–Coulomb, Bingham, Bagnold and Völlmy ones, are included in the code to have the most appropriate description of the rheological behavior of various materials. In the case described here, the soil is a mixture of fine matrix and granular particles, with a size of up to 1 m, and, according to many literature suggestions [21,27,28,29,30], its behavior can be reliably described using the Völlmy model, which is also discussed in Section 4.2.2. Within this limit, it becomes fundamental to choose rheological properties that can simulate the expected bulk behavior of the real landslide [31], as the numerical results directly depend on the chosen rheological model.

3. Perarolo di Cadore Landslide

3.1. Site Description

The Sant’Andrea landslide is located in the municipality of Perarolo di Cadore, Belluno province (North-Eastern Alps, Italy). The landslide position, just upstream of the village center and directly overlooking the left bank of the Boite river (Figure 1), puts this site at high hydro-geological risk. A potential collapse could produce a temporary dam on the river flow, thus causing extensive damage to the downstream urbanized site.

The unstable mass includes an area of 72,000 m² that extends from about 490 to 580 m a.s.l., involving the down-slope part of a larger old landslide that affected the Southern slope of Mt. Zucco. The Sant’Andrea landslide can be described as an active complex slide characterized by slow movements. It involves detrital deposits, about 30 m thick, overlying anhydrite–gypsum rocks. Its activity has highlighted alternate phases characterized by slow movements and accelerations, with a progressive enlargement of the area involved in the last few years. In particular, the relationship between the landslide activity and rainfall has been defined in response to both heavy and lasting rainfall events.

At the end of the 20th century, the first sign of instability was identified by the formation of cracks. A significant rapid movement of about 3000 m³ was registered in November 2000, involving the southernmost portion of the unstable slope. After a period of heavy rainfall, the landslide evolved in paroxysmal propagation that interested the underlying river with the formation of a temporary dam. Later, the movements affected mainly the lower section of the slope, showing a progressive and irreversible worsening of instability conditions. However, the phenomenon no longer manifested other rapid events, involving only the material present in the degradation escarpment.

Until 1985 the area was crossed by National Road 51, replaced in that year by a new section with viaducts, bridges and a tunnel moving on the other side of the valley. Currently, the hairpin bend, visible on the East landslide side (Figure 1c), is part of the minor local road, the SP 42. Even the Treviso-Calalzo railroad, that until 2001 crossed the area, was replaced by a new section that moves in deep in the mountain in a long tunnel. Since the 1990s, numerous stabilization works have been carried out (Figure 1b,c): a retaining linear structure made up of micropiles and soil anchors was built on the South side of the railroad in 1992; two curbs, one of which rested on a base of micropiles, the other on jet grouting piles, were built in 2005. Finally, a reinforcement composed of several jet-grouting piles was installed in 2007, very close to the first linear retaining structures dating from 1992. A surficial drainage system was realized at the rear of an old retaining concrete wall existing on the North side of the original railroad. All these interventions produced only temporary effects, as the landslide seemed to stop for a period. Then, the movements started again and gradually destroyed all the reinforcements. Only the old wall, even if affected by important fissures, still seems to be effectively operative for the support of the uppermost part of the slope.

3.2. Geological Setting

The study area is located inside the geographic region of the Dolomites, in the Southern-Eastern Italian Alps. This Dolomitic zone is characterized by a complex tectonic setting, also known as Giunzione Cadorina [32], in which the main Cenozoic deformation led to the development of structures with almost perpendicular geometry. The study area presents a Carnian (Upper Triassic) sedimentary succession in which the outcropping lithologies exhibit a strong contrast in geomechanical behavior. Triassic formations are mainly represented by dolomitic and anhydrite limestones, dolomites and marls, locally affected by karstic cavities [33]. From the base to the top, the geological sequence is formed by the: Dolomia Cassiana, characterized by massive dolomites; Heiligkreuz Formation, characterized by arenites; Travenanzes Formation, made up of an alternation of carbonate and siliciclastic lithologies with evaporitic intercalations; Dolomia Principale, characterized by dolomites.

Different surveys were carried out over time (1987, 2003, 2008, 2018–2019) to provide direct/indirect and spatially distributed information about the geometrical characteristics of the geological framework, considering the thickness and the spatial distribution of the units (Figure 2). To improve the knowledge about the landslide behavior, a complete 3D geological model of the area is evaluated, first collecting all the available stratigraphic data and integrating the information from boreholes with the results of distributed geophysical investigations, i.e., ERT and seismic refraction surveys [34,35,36,37]. Moreover, the elastic properties of the material involved are evaluated, considering petro-physical characteristics such P- and S-waves velocity or electrical resistivity [35].

The interpretation of the available data lead to the definition of a reference geological model, in which a debris mass, with variable grain size and thickness, overlies a bedrock with alteration and fracturing characteristics variable with depth. The stratigraphic record presents, from the top to the base:

• Layer A: The uppermost part of the debris mass. It is characterized by heterogeneous cohesionless deposits mainly composed of coarse gravel and a matrix of poorly graded sand with local silty levels. The thickness in the landslide body varies around 5–10 m;

• Layer B: Probably formed by the overlapping of different debris flow events. Within this layer, two different sub-units could be identified: layer B1 composed of silty clay deposits characterized by low permeability; layer B2 represented by a fine matrix of sand and gravel intercalations. The thickness of B1 is variable between 15–20 m while B2, not always present in the boreholes, reaches a maximum thickness of 5 m;

• Layer C: An altered and fractured gypsum layer which can be recognized at the base of the detrital mass and attributed to the bedrock below. This alteration layer is characterized by gypsum and anhydrite lithology mixed with a silty-clay sand. The thickness is variable, around 5–15 m;

• Layer D: Corresponds to the anhydrite–gypsum bedrock, represented by a marly-limestone lithology assigned to the Travenanzes Formation. The rock mass presents the deformed anhydrite–gypsum levels as observed in the geological core analysis in 2019 [36];

• Layer E: Characterized by a dolomitic rock, observed externally to the unstable area, on the East side, but never identified inside the boreholes. Considering the geological setting and the available information, it is attributable to the Dolomia Cassiana formation.

Deep hydration processes, caused by the infiltration of meteoric water from the upstream zone of the slope, affect the gypsum-anhydrite bedrock, especially in Layer C. This process leads to a plastic rheology of this lithology [38], involving the bedrock and the overlying deposits in a general gravitative process. Therefore, the hydro-geological condition of the site becomes a relevant predisposing and triggering factor for the landslide movement. The complex hydro-geological condition, variable in space and time, depends both on the stratigraphic setting of the site and the rainfall events. The heterogeneous material that constitutes the 30 m thick surficial layer, composed of soils A, B and C, is characterized by hydro-geological and geotechnical parameters that vary with depth. The ERT investigation performed in 2018 in the upstream area of the slope allow identification of two circulation systems within the moving mass (Figure 3): the upper one in the debris mass, between the permeable Layer A and the silty-clay Layer B1, and a deeper one through the fractured gypsum Layer C.

3.3. Geotechnical Parameters

Different in situ and laboratory tests were performed in various periods to provide an accurate geological-geotechnical description and structural setting of the site [35,36,39].

Table 1 reports some of the data directly obtained from these investigations and used to describe the material behavior in the following analysis. In particular, the shear strength parameters obtained by direct shear tests on layers A, B1, B2 and C and the Young modulus determined with in situ dilatometer tests in layer D are listed.

3.4. Monitoring System

To control the kinematic evolution of the Sant’Andrea landslide, various instruments were employed over time, including inclinometers and both long and short extensometers located in the central part of the landslide area, and a pluviometer for rainfall analysis. Geodetic monitoring was set up in 2013 with a Robotic Total Station (RTS), which records the positions of several topographic targets (actually 40). The RTS is a TM30 Leica instrument located in the city center considered stable.

In the last few years, RTS data have allowed the spatial distribution of surficial movements to be established. Progressive involvement in the instability of upper portions along the slope has been noted, with a consequent worsening of the stability conditions. Currently the slope can be divided into two different kinematic zones, called Area 1 and Area 2 in Figure 4. Area 1 is the uppermost part of the slope, which experiences low displacement; Area 2 is the lower part, where the displacement rates are larger. The displacement data measured in 2019 show a significant increase of slope mobility with respect to 2018; in particular, as shown in Table 2, which reports the annual average displacements for two targets chosen as representative of the two areas, velocities in Area 2 significantly increased in 2019, while Area 1 remained practically unmoved. The landslide evolution in 2019 showed a rapid expansion of the unstable area to the upstream zone of the slope, with the activation of previously stable parts.

Figure 4 also shows the distribution of the total surficial displacements measured in the period from 28 February 2019 to 28 February 2020 and obtained by interpolating the displacement data relating to the target positions using the IDW tool (Inverse Distance Weighted function) in the ArcGIS software. In the same figure, the horizontal velocity rates are reported as arrows, directly in the target positions. The portion Area 2 moves as a rigid block with a quite homogenously distributed velocity with the exception of the central part that shows a greater displacement. The maximum displacement measured in the selected period reaches 226 cm. Of course, larger displacements develop on the front of the slope, where, for safety reasons, it is impossible to install targets, but where, systematically, some small portions of the slope detached and rapidly run up to the Boite river as a consequence of heavy rains [37].

Finally, it is important to note that velocities are strongly affected by meteorological conditions. The landslide exhibits rapid accelerations in consequence of heavy rainfalls 1–2 days long but also of long periods of low or medium entity precipitations. Relaxation periods following the rainy periods can last some months [37]. A procedure for short term forecasting the behavior of a landslide, based on the recognition of the meteorological and kinematical conditions of the site, was developed and it is now under review.

4. Results

4.1. FEM-Based Slope Stability Analysis

4.1.1. Geometry and Soil Properties

The 3D FEM model is built by means of Midas GTS NX 2019 v2.1 [40]. The model geometry and mesh (Figure 5) are created on the basis of information collected from topographic surveys of the area and from the geological and geotechnical investigations previously described.

The slope contour lines are derived from the 0.5 × 0.5 m Digital Terrain Model (DTM), generated from LiDAR data in 2009, provided by the Soil Protection Department of Veneto Region. These contour lines have an elevation step of 1 m and are used to generate the outer surface of the slope.

As previously stated, the software reconstructs a first attempt 3D soil stratigraphy, automatically creating a geometry that interpolates the data of all the available vertical boreholes. Subsequently, this first 3D stratigraphy is slightly corrected in the portions where, due to lack of direct measurements, the code performs inaccurate reconstructions. These slight adjustments aim at obtaining a better qualitative reproduction of orientation and distribution of the monitored displacements.

In order to more clearly visualize the model geometry, some reference elements, such as the hairpin bend of the road and the two concrete curbs are drawn above the model in Figure 5; these elements are not modeled, having no effect on the stability of the slope. The old retaining wall, also visible in the same figure, is instead modeled due to its still existing stabilizing effect. It should be noted that the reinforcements installed in the past interventions, such as micropiles and jet-grouting piles, are not considered in the model because they have been completely destroyed as a result of the large displacements accumulated by the landslide, and they do not currently exert any stabilizing effect. As highlighted in Figure 5c, the layer of fractured gypsum is divided into two parts to better represent the hydrology of the slope which, following rainy events, experiences the presence of deep-water seepage.

All the volumetric elements in the model are simulated using Bricks. The following material constitutive laws are adopted:

• Soils: isotropic elasto-plastic model according to Mohr–Coulomb criterion;

• Retaining wall: isotropic-elastic model.

The geotechnical parameters obtained from standard penetration tests (SPT) and laboratory tests (Table 1) are initially considered for carrying out a back-analysis, as described in the next paragraph. Subsequently, they are slightly modified, in order to determine the resistance values which, most likely, characterize the soils in a collapsing condition (i.e., FS = 1) in the worst hydrological situation. Table 3 lists their final values obtained at the end of this refinement process and definitively adopted for the model.

Note that the fractured gypsum is considered divided in two parts characterized by two sets of parameters. This is for considering the effect of degradation due to hydration processes occurring when the water arrives in contact with gypsum. The decrease of resistance consequently to hydration is specifically studied in a work [41] which evidenced that, with saturation, gypsum cohesion might undergo a 30% reduction.

4.1.2. Hydraulic Conditions

The stability conditions are investigated considering a wet configuration. The presence of two separate seepage systems as previously described and illustrated in Figure 3 are included in the model.

The upper perched aquifer involves only the gravel layer (Layer A) with its base located at the interface between the latter and the silty clay layer below. The piezometric height is arbitrary supposed 2 m above the interface surface.

The second one, placed above bedrock, is an unconfined aquifer that interests only the altered gypsum rock, the gypsum bedrock and the Dolomia formation. The water table of this second aquifer corresponds to the top of the fractured dry/wet gypsum layer (Figure 5c), which may be considered relatively more permeable due to its cavities and fissures. The pore pressure in the stress analysis is included in steady state condition, thus considering a drained behavior of the different soils in saturated conditions.

An SRM analysis is also completed for a dry configuration, with the main purpose of defining the maximum entity of the volumes characterizing the landslide cluster. As previously mentioned, the RTS data evidence strong acceleration of the slope during rain events. This aspect is quite obvious considering that the wet configuration represents the worst hydrological condition, for which the slope would arrive very close to instability. As a consequence, the soil shear strength parameters in the SRM analysis are slightly adapted to obtain the safety factor very close to the unit in wet configuration, as normally done in back-analysis [42].

4.1.3. FEM Results

The SRM analyses show that the critical value of the Strength Reduction Factor (SRFmax), assumed here as FS, is equal to 1.17 for the dry model and drops to 1.002 in wet conditions. The FS in wet configuration is assumed as reliable because it is observed that the landslide movements accelerated subsequently to rains and the water pressure increases. In addition to the SRFmax value, the analysis also provides the values of the incremental displacements, accrued at the failure condition and always referring to the initial stable situation. These displacements, however, must not be interpreted in quantitative terms, given that their size has little meaning, but in terms of distribution that is fundamental to define the volume of the unstable mass.

The distribution of the displacement vectors for the wet case is shown in Figure 6a; it shows a good correspondence with the measurements provided by the topographic station (Figure 4). Figure 6b shows the solid shear strain distribution for the wet scenarios, thus localizing the of shear failure for the slope.

Figure 7 compares the distribution of total displacements at the ground surface and along a central vertical section, obtained at instability condition respectively in wet and dry configurations. This comparison shows that in the dry model the sliding surface is deeper inside the slope and affects a larger area in plan. In total, the unstable volume determined in dry state is greater than the one mobilized in wet configuration. However, the mobilization of this greater volume has a minor possibility of occurring because the wet state exhibits a smaller safety factor. In fact, the presence of groundwater has a destabilizing effect given both by the variation of the effective stresses in the soil, and the assumed reduction of the fractured gypsum resistance. This aspect leads to an instability condition with a lower SRF, highlighting larger displacements but involving a smaller volume than that in dry conditions.

As mentioned earlier, to identify the unstable volumes it may be reasonable to define a threshold value expressed as a percentage of the maximum displacement achieved in the SRM analysis. Several attempts are made to determine an acceptable threshold for both dry and wet scenarios. As shown in Figure 8a for the wet case, increasing the threshold reduces the volume considered unstable, but, on the contrary reducing too much the threshold the unstable portion enlarge exponentially. Values between 7.5% and 15% can be considered reliable; among these, the choice of 10% promotes safety for both scenarios. The choice is made also evaluating the distribution of the solid shear strains in the two models; the localization of relevant values of shear strains confirms the reliability of the selected threshold value. This approach allows estimation of the amount of the unstable volumes, which are about 123,000 and 98,000 m³, respectively in dry and wet models. Removing those volumes from the FEM geometry, it is possible to highlight the unstable areas for the two conditions (Figure 8b). Identification of the sliding surfaces and, therefore, of the unstable mass volumes, allowed the realization of the subsequent modeling phase, which is intended to predict the probable kinematic motion of the landslide. This result is consistent with the trends of the solid shear strains (Figure 6b). The choice of using displacements as a parameter on which to identify the unstable volume is made because the strains do not allow to uniquely define the entire detachment volume, but only the area in which the shear strains are maximum.

4.2. SPH-Based Landslide Propagation Analysis

4.2.1. Geometry of the SPH Model

The geometry of the SPH model contains information about the topography of the studied area, the geometry of the unstable volume (source mass) and the boundary conditions. As input data, the basal topography is represented by a terrain mesh of 2 m resolution, which extends also outside of the landslide area. As for the FE model, also in this case the DTM is used, but enlarging the considered area. The mesh is discretized through a set of nodes of x, y and z coordinates which represent the points related to the basal topography, i.e., the original slope corrected by subtracting the geometry of the unstable mass derived from the FEM outputs (Figure 8b). As the sliding mass does not change during the simulation, a fixed mesh is used.

The source mass is also discretized through a mesh with 1 m resolution. The coordinates of mesh nodes are in this case x, y and h_mass, where h_mass indicates the heights of the unstable mass. The unstable mass volumes come from the FEM outputs in the two cases, and they become here the source masses, which will be involved in subsequent propagation. The basal topography and the source mass colored according to the local thickness of the unstable volume in wet state are shown in Figure 9. Some reference elements (e.g., the old retaining wall, the hairpin bend) are also overlapped to the mesh in the same figure, to help with the visualization of the geometry. The Boite river is also shown in the figure to aid interpretation, but it is not possible to include its presence in the propagation model. Moreover, as the collapse occurs very rapidly, we can suppose that the presence of 1 or 2 m of flowing water in the river cannot have an effect on the collapsing mass or on its propagation.

4.2.2. Soil Model in SPH

From a geo-mechanical point of view, the moving mass has a complex rheology because the granular component can develop frictional behavior while the finer component a viscous one (Figure 10). To consider an equivalent fluid relationship, it is supposed that the sliding mass can be characterized by a frictional regime affected by turbulence in a collisional dynamic. The two-parameter model of Völlmy [43], originally proposed for snow avalanche analysis, is deemed adequate also for modeling granular flows and rock avalanches [44,45,46]. It is consequently chosen here to describe the Sant’Andrea landslide materials. In this rheological model the basal resistance is a result of frictional and turbulent contributions given by Equation (1):

(1)$\tau =\mu \cos \theta +\frac{v^2}{\xi h}$

where μ is the coefficient of kinematic friction, θ the slope angle, v the velocity of the moving mass, ξ the turbulent factor and h the local thickness of the propagation wave. This model derives from the concept of frictional model, largely adopted when the mass movement is slow enough that the particles are kept in almost permanent contact and the moment exchange is mainly due to frictional mechanisms. In order to include the collisional nature of debris flows in the analysis, it was added a term considering energy dissipation due to inner turbulences. Thus, the two parameters to be defined in the Völlmy model are μ and ξ. Due to the complex and heterogeneous nature of the soil involved, specific traditional tests cannot be performed to uniformly describe the material composing the landslide (Figure 10). It is not even possible to consider the parameters included in Table 1 as representative in the propagation analysis. In the SPH method, the mathematical formulation considers the mass constituted by an equivalent material with average properties of all the materials involved, which flows in conditions of large deformations developed very quickly after the overcome of “quasi-static” failure state. Many authors, e.g., [6,11,14,18,19,20,21,22,31], evidenced that the constitutive models implemented for analyzing the “quasi-static” condition and the flow state cannot use the same numerical values for the parameters, even if, sometimes, the adopted model has the same mathematical formulation (e.g., the Mohr–Coulomb model). Consequently, the sets of parameters for the two propagation scenarios are obtained on the basis of the results of a sensitivity analysis performed to values taken from a literature research.

4.2.3. Sensitivity Analysis

The Völlmy parameters are generally chosen on the basis of a calibration obtained with a trial and error technique by reproducing a really occurred event [6,18]. In the case of the Sant’Andrea landslide, only displacements of the topographic targets are available (See Section 3.4), while no data about significant collapse events are available. Consequently, a back analysis is not possible. For this reason, a sensitivity analysis is carried out to better explore the role of each rheological parameter and to evaluate the range of area involved in the propagation of the mass. The ranges of the Völlmy parameters are chosen on the basis of indications provided by the scientific literature [21,27,28,29,30] and according to the grain-size composition and geotechnical properties of the soil involved in the possible collapse. In particular, μ varies in the range 0.36–0.62, corresponding to a friction angle in kinematic conditions varying between 20° and 32°, which may be considered a realistic range for a soil composed of coarse particles immersed in a silty clay matrix, with highly variable relative percentages. For the same materials, ξ ranges from 200 ms⁻² to 800 ms⁻², according to [21].

The sensitivity analysis is carried out considering the volume of the wet FEM case and varying the values of µ and ξ in two separate sets of analyses. In the first analysis set, µ is varied between 0.36 and 0.62, while ξ is maintained constantly equal to 500 ms⁻², which is the intermediate value in the range mentioned above (set 1). In the second set, ξ ranges from 200 ms⁻² to 800 ms⁻², keeping µ equal to a value of 0.45, corresponding to a friction angle in kinematic condition equal to 24° (set 2). Seven simulations are developed for each set, as reported in Table 4.

The height profiles obtained in all the analyses and evaluated along cross-section C–C, indicated in Figure 9, are shown in Figure 11a,b. It is noticeable that the runout distance is more sensitive to the friction angle value, while the turbulence parameter does not strongly affect mass propagation.

In all cases, the sliding mass moves downslope creating a dam in the Boite river. Such a mass generally remains confined in the river; only with the lowest values of the coefficient of kinematic friction (0.36 and 0.40) the mass run to a longer distance, overflowing the river embankment which protects Perarolo village. In these two cases, the sliding mass may cause damages to some houses near the river.

4.2.4. SPH Results

In this paragraph, the propagation of volumes obtained for the wet and dry configurations from 3D-FEM is analyzed with the SPH model. The adopted constitutive parameters are summarized in Table 5. To consider a greater resistance in dry conditions, μ = 0.58 (which corresponds to a 30° friction angle in kinematic condition) and ξ = 200 ms⁻² are chosen. On the contrary, to consider the destabilizing effect due to the presence of groundwater, which leads to a reduction in resistance, in wet conditions μ = 0.45 (24° friction angle in kinematic condition) and ξ = 500 ms⁻² are chosen instead.

Figure 12a,b provides the comparison of the computed evolution of the landslide propagation obtained from the wet and dry simulations at different time intervals. The numerical results are represented considering the same time intervals, i.e., t = 0 s, 10 s and 60 s during propagation. They allow to easy comparison of the runout distance of the sliding mass and the depth of the final deposit over the valley. In wet conditions, the unstable source mass is smaller than in dry conditions, but the landslide propagation involves the whole section of the Boite river with deposit heights of about 15 m. In the source area, just a small quantity of the unstable soil remains because the mass has almost totally collapsed. In this case, even if the sliding mass fills the whole Boite bed, the Perarolo village is only marginally affected by the landslide runout. In dry conditions, due to the higher shear resistance of the mass, more material remains in the detachment zone. The flowing mass stopped quite soon, all accumulating in the lower part of the slope, where there is a progressive increase of the deposit heights with time. The collapsed volume only partially involves the Boite section with deposit heights of about 14–20 m. The final profiles of the deposits along the C–C section for both the dry and wet cases (Figure 13) confirm what can be seen in Figure 12.

5. Discussion and Conclusions

The study of the Sant’Andrea landslide highlights the advantages given by the coupled use of a 3D-FEM analysis and a particle-based model to define the kinematical evolution of the landslides and to predict the area potentially involved in the collapse. The results of the FEM models condition the SPH analyzes and therefore the coupling is intended as a subsequent application of the two methods.

In this study, the information provided by mapping and monitoring of surficial displacements play a fundamental role in constructing a reliable 3D-FEM model. In fact, as normally happens in landslide investigation, even after comprehensive geological and geotechnical surveys, many uncertainties about the 3D setting of deep layers remain. The 3D stratigraphy reconstructed on the basis of borehole data, ERT and seismic surveys is slightly adjusted in the analysis stage to find a better fitting with the orientation and areal distribution of the displacement vectors obtained from the analyses of data acquired by the robotic total station. The kinematic behavior of the model is therefore constrained to describe the observed behavior of the real system. This careful and accurate FEM simulation returns in output fundamental information for subsequent analyses at large deformations. Here, a sensitivity analysis allows the recognition of the area which could be involved in slope collapse, even with the lack of rigorous calibration for rheological parameters.

The hydraulic boundary conditions set on the wet models allow the analysis of several factors. First of all, the presence of a deep aquifer, applied in steady state conditions, changes the stress state of the slope, thus modifying the FS and the unstable volume. Furthermore, a reduction in the cohesion of the portion of fractured wet gypsum is considered to simulate in the FEM analysis the degradation of this material in saturated conditions. The wet propagation SPH model includes instead a reduction of the rheological parameters assigned to the homogeneous equivalent soil, so as to reproduce the greater fluidity of the material in this condition.

The study presented here remains a forecasting of a potential event, not really occurred; in fact, to date, there have been no cases of collapse of significant quantities of material. The collapses occurred so far were in fact very small, with involved volumes not exceeding a few tens of cubic meters. However, this condition is typical of any slope-scale model before reaching the critical conditions of the specific case. In that sense, it is important to underline that the strategy adopted allows reliable forecasting of the evolution of the collapse, especially for what concerns the establishing of the expected volume of the collapsible material and the path that such material can reasonably travel.

The results indicate that the most probable and dangerous collapse can occur when the layer of fractured gypsum located above the anhydrite–gypsum bedrock is saturated, situation that cause a reduction of its geotechnical properties [41]. In this case the detachment of a mass volume is about 20% smaller than the volume mobilized in a completely dry state, but the triggering probability is quite higher. Nevertheless, as the material involved in the collapse in wet conditions is probably less viscous, due to its saturation, its higher mobility increases its runout. This condition allows the moving mass to fill the entire section of the river flowing at the base of the slope. Even if the downstream village is only marginally directly involved in the landslide collapse and subsequent runout, the final situation becomes extremely hazardous for the inhabitants, because the river flow would be surely deviated, inevitably causing the village to flood. It is important to underline that the kinematics of the landslide is strongly influenced by the rainfall, as typically happens in gypsum-rich landslides. In particular, all the cases of significant acceleration of the landslide, with velocity that passes from 0.2 to 5 cm/d, are due to either heavy rainfall or long periods of average rainfall. Furthermore, the analysis of the cross-correlation between rainfall and velocity time series shows that the kinematic response of the slope is very rapid, from a few hours to 1.5 d. For this reason, the wet scenario must be considered as a reference scenario. The dry case is here studied for completeness only.

The results obtained, although related to a specific unstable slope, can be considered as significant for gypsum-rich landslides with the presence of aquifers. In any case, they highlight the usefulness of the proposed approach and provided indications for the treatment of similar cases.

Author Contributions

Conceptualization, L.B.; Methodology, L.B. and D.P.; Software, L.B.; Formal Analysis, L.B., E.C. and D.P.; Validation, L.B., D.P., G.T., S.C. and A.G.; Data Curation, L.B., E.C. and D.P.; Visualization, L.B. and G.T.; Writing—Original Draft Preparation, L.B., E.C., D.P. and S.C.; Writing—Review and Editing, L.B., E.C., G.T. and S.C.; Supervision, S.C. and A.G.; Project Administration, A.G.; Funding Acquisition, A.G. All authors have read and agreed to the published version of the manuscript.

Funding

This research activity was developed thanks to the financial support of the Department of Soil Protection, Veneto Region, Venice.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available on request from the corresponding author.

Acknowledgments

The authors wish to thank: Marco Puiatti (Head of the Department of Soil Protection, Veneto Region, Venice), for the financial support of this research activity, Anna Galuppo (Department of Soil Protection, Veneto Region) for the technical and administrative management of the landslide monitoring system, Alvise Lucchetta and Natalino Zoggia (Department of Soil Protection, Veneto Region), Rocco Mariani (Department of Civil Defense, Veneto Region), Osvaldo Carniel and Cosimo Martinelli (Clios Srl, Belluno) for the technical and logistical support. All the used data are property of the Veneto Region; the authors have obtained permission for their use in this work.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. (a) Location of Perarolo di Cadore (NE Italy); (b) frontal view of the landslide area; (c) orthophoto of Sant’Andrea landslide site and its neighborhood. In the orthophoto the boundary of landslide as well as the rivers Boite and Piave, Perarolo village, the main human infrastructures existing near the landslides and the main reinforcement structures built in the past to stabilize the slope are indicated.

Figure 2. Boreholes and geophysical investigations performed on the landslide area [34,35,36,37]. The stratigraphy of the site is derived from field surveys. The different geological units are reported in the stratigraphic log: sandy gravel deposits (Layer A); debris flow deposits in its silty-clay level (Layer B1) and fine matrix sand and gravel level (Layer B2); the bedrock is represented by the fractured gypsum layer (Layer C) and by the anhydrite–gypsum rock mass (Layer D); dolomitic bedrock (Layer E) [36].

Figure 3. Cross-section A–A of the Sant’Andrea landslide with indications about the hydrogeological setting.

Figure 4. Distribution of surficial displacements of Sant’Andrea landslide, identified in the time span from February 2019 to February 2020 obtained from monitoring of topographic targets [37]. The vectors represent the horizontal displacements while the contours refer to the value of total displacements.

Figure 5. (a) Three-dimensional (3D) FEM model geometry; (b) planimetric view with the position of some reference points; (c) section A–A of the slope.

Figure 6. (a) Displacements vectors for the wet scenario (note that the scale used for the vector color is the same used in Figure 7a,b); (b) solid shear strain for the wet scenarios with FSmax (Note that in the contour bar, the percentages represent the frequency of model nodes for each range of values).

Figure 7. (a) Contour plot of the total displacements in wet configuration (FS = 1.002); (b) section B–B; (c) contour plot in dry conditions (FS = 1.17); (d) section B–B.

Figure 8. (a) Unstable volume in relation with the percentage of maximum displacements for the wet condition; (b) landslide sliding surface in wet model obtained removing the unstable mass. The external limit of unstable area in dry state is also overlapped for sake of comparison.

Figure 9. Three-dimensional (3D) visualization of the SPH model in which the basal topography and the mesh of the unstable mass are represented. The figure shows the model geometry and the unstable volume, with the position of the most important details.

Figure 10. Frontal view of the landslide showing the complex composition of the soil involved.

Figure 11. Height profiles of the sliding mass along cross-section C–C from the sensitivity analysis: (a) simulation with different values of friction angle (set 1); (b) simulation with different values of turbulence factor (set 2).

Figure 12. Some simulated steps (t = 0 s; t = 10 s; t = 60 s) of the landslide propagation in (a) wet and (b) dry conditions.

Figure 13. Height profiles of the sliding mass along cross-section C–C for the simulations in wet and dry conditions.

Table 1

Soil parameters obtained from direct shear tests (layer A, B1, B2 and C) and from dilatometer tests (layer D) performed on the materials involved in the landslide [35,36,39].

Unit	Borehole	Depth (m)	φ (°)	c’ (kPa)
Layer A	C2	9.0–9.2	32	1
	S1	10.5–10.7	29	23
Layer B1	C2	20.8–21.0	29	3
	S1	25.7–25.9	26	30
	S1	29.1–29.4	24	25
	S2	18.4–18.8	27	19
Layer B2	S3	28.0–28.5	30	20
Layer C	C2	41.5–42.0	28	0.2
Unit	Borehole	Depth (m)	E (MPa)
Layer D	PS1	67.5	4353
	PS3	69.5	3929

Table 2

Displacement rates in 2018 and 2019 for two reference targets P1 and P4, where V_H, V_Z and V_T state for horizontal, vertical and total displacement rates, respectively.

Area	Target	Annual Average Displacement (cm/y)
		VH		Vz		VT
		2018	2019	2018	2019	2018	2019
1	P1	1.8	1.8	−1.4	−1.5	2.2	2.4
2	P4	45	80	−44	−63	63	106

Table 3

Soil parameters for FEM analyses, where E, φ and c’ state for Young modulus, friction angle and cohesion, respectively.

Material	E (MPa)	φ (°)	c’ (kPa)
Gravel deposit (A)	100	35	5
Silty Clay (B1)	40	26	30
Fine matrix sand/gravel (B2)	75	30	5
Fractured dry gypsum (C1)	25	24	15
Fractured wet gypsum (C2)	25	24	10
Anhydrite–gypsum bedrock (D)	4000	24	100
Dolomia (E)	8000	40	200

Table 4

Rheological parameter values assumed in the SPH sensitivity analysis.

Set 1	φ (°)	20	22	24	26	28	30	32	ξ (ms−2)	500
Set 2	ξ (ms−2)	200	300	400	500	600	700	800	φ (°)	24

Table 5

Rheological and numerical parameters adopted in the propagation simulation in wet and dry conditions.

Source Mass	Wet	Dry
Volume (m3)	97,800	122,600
Rheological parameters
Friction angle φ (°)	24	30
μ = tanφ (-)	0.45	0.58
Turbulence coefficient ξ (ms−2)	500	200