Reconstruction of the evolution phases of a landslide by using multi-layer back-analysis methods

Back analysis is the most common method to study landslide movements after the event, and it allows us to understand how a landslide evolved along the slope. This paper presents the back-analysis of the Pomarico landslide (Basilicata, Italy) that occurred on January 25th, 2019, on the southwestern slope of the Pomarico hill. The landslide, of rotational clayey retrogressive type—planar sliding, evolved in different phases until it caused a paroxysmal movement in the early afternoon on January 29th, 2019. The landslide caused the collapse of a bulkhead (built at the end of the twentieth century) and of some buildings along the village’s main road. In this paper, a multi-layer back-analysis study is presented, based on the limit equilibrium model (LEM), applying the solution proposed by Morgenstern and Price in Geotechnique 15(1):79–93zh, (1965) and implemented in the freeware software SSAP 2010. The analysis allowed the reconstruction of the entire landslide evolution, using geotechnical parameters obtained from both laboratory and in situ tests, and data from the literature. The application of multilayer back-analysis made it possible to avoid the homogenisation of the layers, modelling the event according to the real conditions present on the slope. The use of the SSAP software made it possible to curb the problem related to the theoretical limitation of the shape of the rupture surfaces, by evaluating independently the friction angle locally and by discarding all those surfaces, which, due to this problem, presented a non-reliable factor of safety (FS) value. The modelling revealed a slope that is highly unstable as the height of the water table changes. The FS calculated under water table conditions close to ground level was less than 1 (FS = 0.98), simulating the first landslide movement (November 2018). The subsequent model reconstructed the critical surface responsible for the January 2019 movement and calculated the FS present on the slope (FS = 1.01). Eventually, the paroxysmal event on January 29th, 2019, was modelled, returning an FS of 0.83, and a sliding surface that sets below the bulkhead, causing its failure. Furthermore, the modelling of the slope in the presence of adequate retaining structures demonstrated the (non-) effectiveness of the retaining wall system represented by the bulkhead. The proposed method of analysis suggests further applications in similar complex multi-layer soil-structure interaction scenarios.


Introduction
Landslides are natural phenomena in which large or small masses of soil and/or rock, making up the slope, slide downwards, sometimes causing human losses, and considerable damages (Petley 2012;Hungr et al. 2014;Perera et al. 2018). These movements affect more or less all geological materials (natural rocks/soil, artificial fill, or their combinations), occurring and developing in a large variety of volumes and shapes (Hungr et al. 2014). They can be triggered by various factors, such as rainfalls, earthquakes, and anthropogenic activities (Song et al. 2021). Knowing their evolution, the properties of the soils involved and their mechanical behaviour, can allow to understand the stability or not of the slope (Memon 2018;Pazzi et al. 2019). Nevertheless, to reach this knowledge, it is required the availability of a wide range of observations, measurements, and data, and the definition of geologic and hydrologic conditions related to the phenomena occurrences .
To study the stability of a slope and to define the parameters that characterize its displacement, is commonly used the stability analysis (Memon 2018). Its main objectives are (1) the calculation of the factor of safety or safety factor (FS), defined as the ratio of the soil shear strength to the shear stress (Jiang and Yamagami 2008) and (2) the determination of the critical sliding surface/s. The stability analysis involves the study of the forces that oppose the movement (mainly defined by the mechanical characteristics of the ground, e.g., the friction angle φ' and the cohesion c') and the forces that cause the movement. In the field of geotechnics, this methodology is mainly applied to study the maximum load that can be supported by a slope or a geo-structure without causing subsidence (Sloan 2013), while in the landslide studies it is usually applied after the event to perform a back analysis that enables to determine the conditions that had caused the slope failure (Fidolini et al. 2015;Garduño-Monroy et al. 2020;Morelli et al. 2020). Slope failure implies that the FS at the time of landslide is equal to the unit (FS = 1.0) (Wesley and Leelaratnam 2001;Jiang and Yamagami 2008;Popescu and Shaefer 2008;Zhang et al. 2010). Values of 1 or higher mean that the slope is stable, while values lower than 1 mean that the slope is unstable (Duncan et al. 2014). By knowing the FS and adopting the appropriate method of analysis (see later in this section), back analysis allows to build a model that represents the slope at the time of collapse, reconstructing the equilibrium conditions just before the collapse, and allowing to estimate the residual resistance (Duncan et al. 2014). This analysis is based on the geometry and hydrogeological conditions of the slope before the failure (Duncan et al. 2014). This method is widely used to characterize the geotechnical parameters (φ' and c') of the slope and to design landslide slope rebuilding works (Jiang and Yamagami 2008;Popescu and Shaefer 2008).
Back analysis is a good procedure to estimate the shear strength of soil mobilised by slope failure (Tang et al. 1999;Duncan et al. 2014). A back-analysis study involves knowledge of the moving mass (before and after rupture), geometry, and actual profile of the sliding surface. However, the geometry of the rupture surface or the mass weight parameters is not always known. Estimation of the slip surface can be performed by applying various stability analysis approaches among which the most used are the limit equilibrium method (LEM) (Zhou et al. 2019) and finite element method (FEM) which analyses use the minimum FS, but the estimated surface/s may be different from the real one (Brunetti et al. 2014;Saeidi et al. 2016;Jiang et al. 2020).
The FEM methods are widely used in the geotechnical field for landslide prediction and slope stability prediction (Griffiths and Lane 1998;Sloan 2013;Duncan et al. 2014;Carlà et al. 2021). In FEM methods, the slope is assumed to be a continuous elasto-plastic medium. It is subdivided into a grid composed by many small elements in which the tension and displacement components are identified by the values they assume at the vertices of the elements. It is therefore possible to know the state of tension and deformation at each point. In this type of analysis, no assumptions are made about the shape and position of the sliding surface, but failure occurs when the tensions exceed the limit stresses. The analysis is conducted incrementally, i.e., by increasing the external load or decreasing the soil strength, until failure is reached. This involves a rather complex calculation with many parameters comprised in the calculation.
The LEM is the most widely used method in the field of stability analysis and posterior analysis to estimate shear strength parameters (Krahn 2003;Jiang and Yamagami 2008;Duncan et al. 2014;Deng et al. 2017;Liu et al. 2019;Liu et al. 2020;Beyabanaki 2021;Firincioglu and Ercanoglu 2021;Abdulai and Sharifzadeh 2021). The main advantage of using LEM is dictated by its simplicity, formulation complexity, shorter analysis time (Azarafza et al. 2021), and wide use which has led to the development of different software that increased its reliability. The main disadvantage is the assumption of the rupture surface along which the movement will develop (Duncan et al. 2014). The LEM is a calculation methodology that relies on the Mohr-Coulomb rupture criterion to evaluate the FS (Duncan 1996;Liu et al. 2015;Memon 2018). This method assumes a rigid and perfectly plastic behaviour of the soil, i.e. it assumes that the soil does not deform until the failure condition is reached.
The LEM method, based on the slice analysis approach, divide the slope into slices, bounded by vertical surfaces, above the assumed slip surface (Krahn 2003;Duncan et al. 2014;Azarafza et al. 2021;Firincioglu and Ercanoglu 2021). The theory of the slice method is based on solving static equilibrium by considering the acting and resisting forces on each slice (Duncan 1996;Krahn 2003;Firincioglu and Ercanoglu 2021).
The equilibrium equation is solved on each slice by summing the forces and moments and is calculated for the entire moving mass. If all assumptions in the equation are solved and thus satisfied, the mass is in equilibrium, and the analysis is satisfied (Azarfza et al. 2021;Firincioglu and Ercanoglu 2021). A unique factor (FS) is obtained by calculating the ratio between resistance forces and moments acting on each slice (Firincioglu and Ercanoglu 2021).
Various methods for solving the equation at equilibrium have been proposed in literature, and most of these assume that the point of application of the normal force acting at the base of the slice is known. Methods that consider a circular sliding surface assume that the equilibrium of moments is around the centre of the circle for the entire mass divided into slices. In contrast, solutions that study non-circular sliding surfaces consider the equilibrium for each of the individual slices (Duncan et al. 2014). These methods differ mainly in the number and type of assumptions adopted to solve the equilibrium equation and can be divided into two main groups: non-rigorous and rigorous solutions. In non-rigorous solutions (Fellenius 1936;Bishop 1955;Janbu 1954), the number of assumptions used to solve the problem is greater than those required, thus making the problem overdetermined. In rigorous solutions (Morgenstern and Price 1965;Spencer 1973;Sarma 1973Sarma , 1979Chen and Morgenstern, 1983;Janbu 1973), all equilibrium conditions (equilibrium of vertical and horizontal forces and equilibrium of moments) are satisfied and are therefore based on a strictly necessary number of assumptions (Duncan et al. 2014;Firincioglu and Ercanoglu 2021). In the "Methods" section, a comparative table to summarize these points is presented.
In agreement with Fredlund and Krahn (1977), the abovementioned solutions provide similar values in FS. The solution proposed by Morgenstern and Price (1965) satisfies all equilibrium conditions, and it can be applied in all cases where interface forces play a significant role in slope stability (Duncan et al. 2014). Furthermore, the Morgenstern and Price method has been implemented in numerous free (such as SSAP by Borselli 2020) and widely used commercial software (such as SLOPE/W by Geostudio 2020; ASPEN 2000 by Newsoft 2019).
There are numerous studies in literature that report the application of back analysis for the determination of geotechnical parameters and the sliding surface of a landslide (Jiang and Yamagami 2008;Akin 2013;Saeidi et al. 2016;Jiang et al. 2020;Rana and Babu 2022).
This study applies the classical back-analysis method, using the LEM method and the solution proposed by Morgenstern and Price (1965), on a multi-layer slope to identify the phases that led to slope failure. The multilayer back analysis, presented here, involves the non-homogenisation of the layers, i.e., it considers each individual layer with its own geotechnical parameters, supported by both laboratory and in situ tests, because the layers can have very different parameters and there often are situations in which the upper or lower layers have worse geotechnical characteristics than the surrounding ones. One of the main problems encountered in the multi-layer slope stability analysis is the theoretical limitation in the shape of the failure surfaces, due to the theory of soil thrust and the computational problems of the FS calculation (Ching and Fredlund 1983;Chowdhury and Zhang 1990;Stianson et al. 2015;Borselli 2020). The inclination of the failure surface identified by analysis must locally respect the Mohr-Coulomb failure criterion in areas where active and passive thrusts prevail. In the presence of high lithological variability (multi-layer system), the values of the limiting angles of inclination can change strongly within the slope. Therefore, it is not suitable to provide a unique fixed limit value for the inclination angles in the passive and active thrust zones. The freeware software SSAP, used in this study, allows to assess the compatibility in terms of geometry and friction angle of the upper and lower layers, as it not only calculates the FS and identifies the potential slip surface, but also evaluates the friction angle locally, considering the variation of the layers, that is an essential parameter for the success of a multi-layer back-analysis study of a landslide movement evolved in several phases. Moreover, before proceeding in calculating the FS, SSAP software automatically verifies each surface, discarding those that locally violate the limiting Landslides slope angles in the strata where the potential sliding surface is present, as they would present unreliable FS values (U.S. Army Corps of Engineers 2003; Duncan et al. 2014;Borselli 2017Borselli , 2020.
In the "Study area" section, the landslide is described, while more details on the method and software are presented in the "Methods" section. The results of the retrospective analysis of the four

Technical Note
identified landslide phases are presented in the "Results" section. A discussion of the results is provided in the "Discussion" section.

Study area
Pomarico village (MT, Italy, Fig. 1) is located at 454 m a.s.l. on an elongated ridge in a NW-SE direction between the valleys of the Bradano and Basento rivers. It is in the south-eastern part of the Bradanian depression, between the Apennine mountains and the Apulian foreland. Since the Lower Pleistocene, this area has been subject to general tectonic uplift (Balduzzi et al. 1982;Bozzano et al. 2002). During the Plio-Pleistocene, the Bradanian depression was filled by a powerful sedimentary succession of clastic origin (Gasperi 1995).
The Pomarico hill presents the complete regressive sequence of the Bradanic cycle (Cherubini et al. 1985). These deposits are represented by a regressive sequence of age between the lowermiddle Pliocene and middle Pleistocene (Guerricchio and Melidoro 1979;Doglioni et al. 2020). They are made up at the base by grey-blue clays (Sub-Apennine Clay Formation, light green in Fig. 1) (Genevois et al. 1984), transitioning upward to sands (Monte Marano Sand Formation, brown and yellow area in Fig. 1) (Cherubini and Walsh 1982) and conglomerates (Irsina Conglomerate Formation, orange area in Fig. 1) (Balduzzi et al. 1982;Perrone et al. 2021). The conglomerates represent the closing deposits of the Bradanic sedimentary cycle (Cherubini and Walsh 1982) (Fig. 1).
Pomarico's slopes and valleys are characterized by large outcrops of detrital deposits originated by landslides and erosive processes of sandy-clay layers (De Marco and Di Pierro 1981; light brown in Fig. 1). The upper part of the slopes is characterized by In this solution, the scaling coefficient acts as an overall scaling factor sandy debris and the central part by clayey debris. These debris bodies are characterized by a variable thickness (between a few meters and 10-15 m.), and because of the direct infiltration of rainwater and the multilayer groundwater flow, they are involved in slope movements (Bozzano et al. 2002;Perrone et al. 2021;Doglioni et al. 2020). The sides of the hill are deeply incised, almost up to the height of the built-up area, by badlands embedded in the grey-blue clays (Guerricchio and Valentini 1979;Genevois et al. 1984) and by furrows and niches of detachment of ancient landslides (Perrone et al. 2021). The Pomarico hill is asymmetrical with irregular profiles caused by the arrangement of strata (dip slope in the southwestern sector and anti-dip slope on the opposite side) and the different erosive power of the two main watercourses, namely Fosso Pezzillo (SW side) and Canale Santa Croce (NE side) ( Fig. 1) and their tributaries that affect the slopes with a dense network of incised and narrow ditches (Perrone et al. 2021).
The slopes of Pomarico are affected by numerous quiescent or evolving landslides that surround much of the village. They are primarily composite retrogressive roto-translative landslides sensu Cruden and Varnes (1996)  The landslide, which affected the southwestern slope of Pomarico on January 25th and 29th, 2019, is part of deeper slope instability. This deep slope is characterized by ancient and recent landslides (mostly quiescent) and a debris cover that has variable thickness and capacity to accumulate water on multiple levels (Bozzano et al. 2002;Perrone et al. 2021). The landslide in 2019 was triggered by heavy rainfall during autumn 2018 and by the long rainfall event that occurred during the night on January 24th-25th, 2019 (Doglioni et al. 2020). It is a retrogressive clay rotational -planar slides (sensu Hungr et al. 2014) that evolved in multiple stages (Hungr et al. 2014;Doglioni et al. 2020;Perrone et al. 2021).  Table 1 to understand the  table under "LEM method  selection"   Table 2 Parameters relating to the ties and piles used for the stability analysis (β = angle formed by the tie rod with the horizontal; L = length of tie rod; T = design load; L c = cemented length; D = diameter of the piles making up the piling; D1 = spacing between piles; D2 = opening length between piles; Cu = equivalent undrained shear strength value of the bulkhead; γ = equivalent unit weight to be applied to the bulkhead)

Technical Note
The first movements occurred at the end of November 2018 and resulted in a ground slip (soil slide sensu Hungr et al. 2014) in front of the bulkhead piles built at the end of the twentieth century downstream of the main street of the village (Corso Vittorio Emanuele II) (Doglioni et al. 2020). The movements then resumed on January 25th, 2019, likely triggered by activation of the clay level placed at the transition between sands and sub-Apennine clays (Doglioni et al. 2020). It was arrested by the bulkhead, which was partially denuded and began to deform at the base. Subsequently, further movement occurred on January 29th, 2019, causing the collapse of both the bulkhead and the buildings constructed along Corso Vittorio Emanuele II (Fig. 2) (Doglioni et al. 2020;Perrone et al. 2021).

Methods
To perform the back-analysis study presented in this work, the SSAP 2010 (Slope Stability Analysis program) software version 5.0.2 (2021) (https:// www. ssap. eu), was employed. SSAP 2010 is a Fig. 4 a Results of the stability analysis. In yellow are shown the 10 critical surfaces identified by the analysis, while in red is highlighted the critical surface ([1], [2], [3], and [4] are layer 1, layer 2, layer 3, and layer 4, described in Table 4); b map of local FS (blue to green and yellow colours show stable areas with a FS equal or higher than 1, orange to red and purple colours show unstable areas with a FS lower than 1) Landslides comprehensive freeware software for verifying the stability of natural, artificial, or reinforced slopes (Borselli, 2020). It allows verifying the stability of slopes in loose soils and/or with fractured rock masses, and/or under liquefaction conditions, using LEM (Borselli, 2020). As inputs, the SSAP software needs a slope section been identified and the geotechnical characteristics of each soils present and involved been assigned. Moreover, the SSAP software allows to model anthropic/structural elements like piles, bulkhead, road, and buildings as overload.
SSAP allows performing in-depth stability checks, choosing among 7 different rigorous computational methods based on LEM. The seven LEM methods are Janbu (1973); Spencer (1967); Sarma (1973); Morgenstern and Price (1965); Chen and Morgestern (1983); Sarma (1979); Borselli (2020). Table 1 shows the characteristics of the different LEM methods, found in literature. Only the rigorous solutions are available in the SSAP software. Figure 3 shows the flowchart of the analysis process based on the following three main steps: (a) identification of the input parameters, (b) choice of the most suitable LEM method based on the available data, (c) limit equilibrium stability analysis to identify the FS and the failure surface/s. Among the computational methods available in SSAP, the solution chosen to carry out the back analysis of the Pomarico landslide was the one proposed by Morgenstern and Price (1965), in the absence of seismic shocks, due to its high reliability. In fact, it uses an arbitrary function to define the direction of the resultant of the interface forces. This solution represents a fully balanced method (Duncan et al. 2014;Borselli 2020). Morgenstern and Price (1965) developed a solution in which the relationship between vertical and horizontal interface forces is expressed through an unknown function, whose parameters are derived from the global mass balance. The solution is fully balanced and can be used to study any failure surface with any shape. Compared to previously described methods, the solution proposed by Morgenstern and Price (1965) does not have a simple final formula, but a system of non-linear equations to be solved by iterations. The method is characterised by high accuracy, reliability, and numerical stability as it satisfies both the balance of moments and forces, and allows for the inclusion of complex stratigraphy and the presence of reinforcement works (Morgenstern and Price 1965;Duncan et al. 2014;Firincioglu and Ercanoglu 2021). The LEM is also used to calculate the FS according to the iterative algorithm devised by Zhu et al. (2005). The calculation of the FS is generalised for any LEM method by means of a new algorithm (Borselli, 2020). In addition, the Sniff Random Search surface tool was used. Sniff Random Search engine tends to maximize, during random surface generation of the surface, the crossing into the layers that have the poorest strength characteristics. Regarding the presence of structural elements as stabilizing piles, the SSAP software uses the formulation by Ito and Matsui (1981) for calculating the maximum mobilized reaction force offered by the pile under conditions of plasticization of the soil-pile interface. In this case, the formulation corrected by Kumar and Hall (2006) for very closely spaced piles was used, as it is more conservative (Borselli, 2020). The anchor was included in the model according to the data given in the project for its construction. Table 2 shows the parameters for the piles and the tie rods used in the modelling. The tie rods, according to the project, were not pretensioned (passive anchor). Table 3 Tests carried out on the four samples collected from the Pomarico landslide. For their locations see Fig. 1 Tests Samples S1 S2 S3 S4 Property index X X X Particle size analysis X X X X Atterberg limits X X X X Direct cut X X X X Annulare cut X   Concerning the presence of anthropic works upstream of the bulkhead (roads, walls, and buildings), the SSAP program allows the insertion of overloads, as mentioned before, to increase the loads acting on a potential sliding surface. In this case, the insertion of a vertical overload allowed considering all the conditions present on the slope.
Considering the type of construction, the overload was assigned equal to 50 kPa, equivalent to about 2.5 m of extra soil.
The height of the water table was inferred from piezometric measurements reported in Bozzano et al. (2002).
As a result of the modelling, SSAP produces a series of graphs. For the purposes of the back analysis presented in this work, only two of them have been selected and an example is shown in Fig. 4, carried out with the parameters given in Table 2 and the geotechnical parameters given in the following section (the "Lithological model" section),. In the upper panel, 10 critical surfaces identified by the analysis are shown in yellow, and the absolute minimum FS surface is highlighted in red. At the end of the global stability check, information about the recalculated local FS value is recorded for each surface on time. Through a specific algorithm, a 2D colour map of the local FS is returned (Borselli 2020), and it is shown in panel b (Fig. 4). The 2D colour map of local FS distribution allows highlighting areas where rupture may occur, i.e., areas that exhibit an unfavourable combination of stress states and therefore have a lower average local FS (warm colours) (Borselli, 2020).

Lithological model
To characterize (i.e., determining soil physical properties and geotechnical parameters) the soils involved in the movement, and to define the input parameters needed by SSAP, geotechnical investigations were carried out. In detail, laboratory tests (property index, particle size analysis, Atterberg limits, and direct and annular cuts) were done on samples taken at representative points of the landslide body and nearby (green dots marked by S in Fig. 1). In addition, in situ permeability tests were performed at the same locations. Table 3 lists the tests performed on each sample, while the results are summarized in Table 4.
The samples analysed are representative of soils characterized by heterogeneous grain size. They consist mainly of sand with similar percentages of silt (S2 and S3), and silt with similar percentages of sand and clay (S1 and S4.) The sand-silt fraction constitutes more than 90% of the material taken outside the landslide body (S2 and S3), while the silt-sand/clay fraction constitutes more than 90% of the material taken along the landslide body. This leads to classify the soils taken in the median part and at the foot of the landslide in silty clays and in clayey sands (SM-SC), while the soils taken outside the landslide body (near the middle-high zone of the latter) as silty clays (CL). Data from the consolidated and drained direct shear tests, in agreement with Bozzano et al. (2002), indicated average friction angle and cohesion values of 22.5° and 22.5 kPa for the silty clays and average friction angle values of 36° for the clayey sands.
The section (black line in Fig. 1) subjected to verification includes the middle-high part of the slope, i.e., the most representative area of the landslide movement (Fig. 1). Thus, based on (a) the inspections, (b) the geotechnical parameters obtained from the laboratory tests, (c) data available in literature, and (d) the stratigraphies obtained from the boreholes, the geotechnical model for stability analysis was created and reconstructed according to the specifications of the SSAP program (Borselli 2020). It is a schematic model composed of 4 lithologies: layer 1-consistent sand (Monte Marano sands formation), layer 2-blue-grey clayey silts, layer 3-remodelled sand, and layer 4-remodelled silts. Geotechnical parameters were assigned to each lithology/layer, and they are shown in Table 5. They were obtained from laboratory tests, and they agree with the parameters reported in Bozzano et al. (2002). The clay substrate was considered impermeable, and the height of the water table was inferred from piezometric measurements reported in Bozzano et al. (2002).

Results
The back-analysis study was divided into 4 different phases with the objective of reconstructing the main phases of the landslide movement. The first analysis (in the following first phase) allowed to identify the critical conditions present on the slope (Ausilio et al. 2001). The second phase made it possible to reconstruct the instability conditions of the slope, considering the presence of structural works. The third phase allowed to study the gravitational movement occurred at the end of November 2018, while the fourth phase allowed studying the movement of 25th-29th January 2019.

First phase (phase I)
The first phase was aimed at identifying the critical conditions present on the slope before the 2018 event. For this reason, the topographic profile was retrieved from Bozzano et al. (2002), and the modelling Table 6 Parameters relating to the ties and piles used for the stability analysis (β = angle that the tie rod forms with the horizontal; L = length of tie rod; T = design load; L c = cemented length; D = diameter of the piles making up the piling; D1 = spacing between piles; D2 = opening length between piles; Cu = equivalent undrained shear strength value of the bulkhead; γ = equivalent unit weight to be applied to the bulkhead)  was carried out in absence of structural elements and considering only the presence of the overload, representing the built-up area of Pomarico. This choice is based on the consideration that the triggering factor of the landslide can be traced back to the heavy rainfall, occurred in the previous months, on a slope already made strongly unstable by its geomorphological configuration and by the presence of areas with high water content (Perrone et al. 2021;Doglioni et al. 2020). The present phase involved performing a series of iterations, varying the height of the water table. This allowed the reconstruction of the FS under in deep and shallow water table conditions. Two representative models of this analysis, both without structural elements, are shown below: a model with deep water-table conditions and one with shallow water-table conditions.

Model without structural elements with deep water-table conditions
The model in the absence of structural elements and with the water table placed at greater depths (Fig. 5) produced a minimum FS of 1.36, with a range between 1.36 and 1.37. These values are quite high and indicate a slope in stable condition. Analysis of the results shows that, the 10 most critical slip surfaces range in length from 138 to 197 m, with a maximum depth of 18 m above ground level, and involve layers 3 and 4 (Fig. 5a). The local FS distribution map (Fig. 5c) shows that the areas with lowest FS are concentrated entirely at the downstream surface of layer 3, near the contact with layer 4 (A in Fig. 5c), and in contact with layer 2 (B in Fig. 5c). The first ones (A) indicate small surface gravitational movements typical of Calanque areas. The second ones (B) are caused by the tensional forces that develop in contact with layers having different geotechnical characteristics, in absence of interstitial pressures in the base layer. These indicate strong variations in FS values as a function of pore pressure.

Model without structural elements with shallow water-table conditions
In the analysis with the water table positioned near the ground level, the model was made with the same configurations as the previous one (the "Model without structural elements with deep water-table conditions" section). The water table height, in this case, was between 3.40 and 5.20 m from the ground level, with a minimum value of 0.60 m (obtained from piezometric measurement reported in Bozzano et al. (2002). The model (Fig. 6) shows that increasing the height of the water table and bringing it closer to the ground level, generates a reduction in FS that falls to 0.98 (meaning unstable slope), with a range between 0.98 and 0.99. The failure surfaces are in the upper part of the slope, affecting layer 1, where they are generated, and layer 3 where they end in contact with layer 4 (Fig. 6a). The 10 most critical sliding surfaces have a length between 134 and 173 m and a maximum depth from the ground level of 20 m.
In the two-dimensional distribution map of FS (Fig. 6c), it is possible to note how the areas with lower FS are located in the zone of contact between layer 1 and 3 (A) and along the contact between layers 3 and 2 (B). In these areas, a FS lower than 1 is recorded; therefore, a strongly unstable slope emerges, especially in the zone of contact between layers.
The stability analysis conducted in the absence of structural elements shows the transition from stable to unstable conditions of the slope, as the height of the water table changes. A FS value lower than 1 is obtained in shallow water-table conditions with the same geotechnical values of the materials used in deep watertable conditions. The above analysis made it possible to obtain a series of critical surfaces representative of the conditions of the slope before the construction of the consolidation work and made it possible to know the factors and mechanisms regulating stability/instability.

Second phase (phase II)
Having determined the most critical slope conditions and having identified the corresponding FS through the analysis of the first phase, a model was constructed to consider the presence of the structural work. The model was built with the configurations (geotechnical parameters and height of the water table) shown in the "Model without structural elements with shallow water-table conditions" section, and the FS values obtained are 0.98 with a range between 0.98 and 0.99. The failure surfaces are located downstream of the bulkhead (Fig. 7a), and from the FS map (Fig. 7c), the areas with lowest FS are all concentrated in the area downstream of the bulkhead (A). The 10 most critical slip surfaces range in length from 123 to 134 m and have a maximum depth from ground level of 17 m.
The model in Fig. 7 shows the slope in its most critical condition with the structural interventions performed in the late twentieth century (Doglioni et al. 2020). Despite the structural intervention, aimed to improve slope stability, the FS is less than 1, and the slope downstream of the bulkhead remains unstable in these given conditions of height water table.

Third phase (phase III)
Subsequently, a representative model of the landslide movement that occurred at the end of November 2018 was produced, when heavy rainfall in October made the slope below the bulkhead highly unstable and triggered fractures in the ground downstream of the piles.
The model was built by varying the topographic profile at the bulkhead, according to the critical surfaces that emerged from the analysis reported in the "Second phase (phase II)" section, thus obtaining a representative profile of the situation at the end of November 2018, when the landslide started to affect the area, directly downstream of the bulkhead and caused a first denudation of the piles. The model was constructed according to the configurations reported in the "Second phase (phase II)" section.
The critical surfaces extend beyond the piling to the upstream buildings (Fig. 8a). The generated FS is 1.01, with a range of FS between 1.01 and 1.09.
The local FS map (Fig. 8c) clearly shows that the lowest FS zones are located in the zone of contact between layers 1 and 3 (A) and between layers 2 and 3 (B).
The whole area around the bulkhead has therefore a rather low local FS, which would explain the detachment of soil, the denudation of the piles and the beginning of the downward rotation of the piles (Doglioni et al. 2020).

Technical Note Fourth phase (phase IV)
Following the landslide at the end of November 2018, a subsequent landslide movement occurred on January 25th, 2019, resulting in the collapse of a further portion of the ground downstream and a deep denudation of the piles. By reconstructing the topographic surface from the results of the analysis shown in the "Third phase (phase III)" section, it is possible to model the January 25th event and explain the collapse that occurred on January 29th, 2019. The detachment of the ground in front of the bulkhead on January 25th, 2019, made the slope even more unstable, and on January 29th, 2019, at 3:15 pm, a paroxysmal slip occurred, leading to the failure of the bulkhead with the overturning of the bulkhead itself and the collapse of the buildings present along Corso Vittorio Emanuele II.
The model was created assuming a change in morphology caused by the movement of the ground downstream of the bulkhead. The FS is even lower with a value of 0.83 and a range between 0.83 and 0.89. The most critical slip surfaces pass underneath the piles, causing them to overturn (Fig. 9a). The local FS map (Fig. 9c) clearly shows the presence of high buckling zones (FS < 1) within layers 3 and 1 (A), especially in the area in front of and behind the bulkhead.

Discussion
This work analyses and describes the landslide events of Pomarico in 2018 and 2019 through the back analysis carried out with the software SSAP, which allowed to model a multi-layer slope and reconstruct the various phases in which the landslide evolved. The multilayer backanalysis, in fact, considers the layers non-homogeneity, attributing to each layer its own geotechnical parameters. This kind of analysis, therefore, must be supported by laboratory and in situ tests to obtain these parameters and cannot be performed using literature values. The higher the number of samples/in situ tests, the better the mean values derived to perform the analysis.
The preliminary back-analysis phases (phase I, the "Model without structural elements with deep water-table conditions" and "Model without structural elements with shallow water-table conditions" sections) allowed the identification of the critical conditions of the slope and therefore closest to those present during the landslide event. Initially, several iterations were performed by varying the height of the water table, considering only the slope with its geotechnical configuration and overburden (in the absence of structural elements). These models made it possible to reconstruct the stability and instability conditions of the slope as the height of the water table varied, obtaining an FS value < 1 in the case of a water table close to the surface with the same geotechnical values of the materials.
Once the slope instability conditions were determined, the slope stability in the presence of the bulkhead was considered. The results show an unstable slope with an FS of 0.98 (phase II). However, despite a rather low FS, the structural works were able to protect the settlement upstream of them. The failure surfaces stopped at the bulkhead, not progressing further, but causing an increase in stresses downstream of the bulkhead that resulted in the formation of fractures in the soil at the bulkhead.
The model used to simulate the late November 2018 event (phase III) shows a slope with an FS of 1.01, with the failure surfaces no longer stopping at the piles but advancing beyond the bulkhead towards the lower portion of the village. The local FS map (Fig. 8c) clearly shows the presence of an FS value below 1 in the portion in front of the piles, where the ground began to slide downwards and cause the progressive denudation of the piles, as well as the beginning of the downward rotation of the pile itself.
The sliding on January 25th (phase IV) led to a further collapse of the ground in front of the piles, causing further denudation. At this point, the piles began to increasingly deform and rotate downwards (Doglioni et al. 2020). The slope maintained this stability until January 29th, 2019, when the collapse of the bulkhead and upstream buildings occurred.
It is important to remember that stability analyses, and the back-testing procedure, are carried out to restore landslide slopes with reinforcement works, such as piles. As reported in Popescu and Schaefer (2008), for the stability of the pile bulkhead, it is necessary to consider both the driving force and the resisting force acting on each pile. An accurate estimation of the lateral force acting on the pile is crucial since the stability action of the bulkhead and the slope are in conflict. This means that a safe stability analysis for the slope does not in turn imply a certain stability for the bulkhead. For a bulkhead to perform the role for which it has been designed, each row of piles must be driven into the slope at an adequate depth across the potential sliding surface. Each pile must be inserted firmly and deeply so that there is adequate mechanical strength, leading to an increase in the FS of the slope, thus ensuring its stability (Xiao et al. 2018;Borselli 2020). For these reasons, a model was created with the same geomechanical characteristics as the previous ones but with a greater extension of the piles. These were lengthened to a total length of 30 m. The bulkhead was built to a total length of 30 m, the upper 24 m of which were placed within layer 1 for passive reaction and the remaining lower 6 m were placed in layer 2 with a blocking function, because it has better geotechnical characteristics than layer 1. Table 6 shows the values assigned to the piles in this modified configuration.
Results are shown in Fig. 10. The resulting FS is 1.34, and the slope surfaces have a length ranging from 1.44 to 1.51 m. Therefore, it is possible to assess that under these conditions, i.e., by extending the pile length up to 30 m, the slope will result very stable. This result means that the reinforcement construction projected was correct in terms of type and shape, but it was under dimensioned in terms of depth of piles and total length.
As stated before and in the "Introduction" section, the main limit of performing a back analysis on a multi-layer slope to identify the phases that led to slope failure, is the need of accurate geotechnical parameters, because the layers can have very different parameters. Therefore, a back-analysis study is successful if it is preceded by a thorough terrain study, which justifies the parameterisation included in the analysis software. Moreover, as mentioned earlier, another limitation of multi-layer back analysis lies in the lack of a detailed slope geometry. This drawback can be overcome using geophysical techniques, like electrical resistivity topographies, Landslides active and passive seismic, and ground penetrating radar just to mention some, in addition to traditional intrusive methods, since they allow reconstructing the subsoil landslide geometry (Pazzi et al. 2017.

Conclusion
In this work, a multi-layer back analysis was successfully employed to study the landslide movement affecting the village of Pomarico. The input parameters of the model, i.e., the mechanical properties of the soils (obtained through both in situ and laboratory tests), the pre-event topographical surface and the water table obtained from the bibliography, made it possible to study the slope, returning a solid multi-layer back analysis.
The use of the SSAP 2010 software made it possible to identify the real sliding surface and the FS of the slope in each evolutionary phase. The proposed approach includes four phases: phase I, identification of the critical conditions presents on the slope by varying the height of the water table; phase II, evaluation of slope stability in presence of structural works; phase III, modelling of the first landslide movement; phase IV, reconstruction of the landslide paroxysmal event. An additional simulation was performed varying the supporting works, in terms of depth and total length, bearing in mind that to perform the function for which the bulkhead is designed, each row of piles must be anchored to the slope at an adequate depth across the potential sliding surface.
The results obtained from each model show that the sequence adopted in the a posteriori analysis is necessary to carry out an adequate study of a slope characterised by a complex sequence of clayey and sandy layers, the movement of which is triggered by variations in the water table, which plasticise the clays and create a sliding surface. The multi-layer back-analysis method implemented using SSAP software allows to reduce the oversimplification of a classical back-analysis based on single strata that may produce, in complex stratigraphic situations, possibly unreliable results. The constraints on of local slope values of sliding surface, associated to each stratum, allow a reliable reconstruction of multi phases instability process. The proposed analysis method suggests further applications in similar complex multi-layer soil-structure interaction scenarios.