Sph propagation modelling of an earthflow from southern italy

XII International Conference on Computational Plasticity. Fundamentals and Applications COMPLAS XII E. Oñate, D.R.J. Owen, D. Peric and B. Suárez (Eds) SPH PROPAGATION MODELLING OF AN EARTHFLOW FROM SOUTHERN ITALY L. CASCINI*, S. CUOMO*, M. PASTOR†, M. COPPOLA* * Laboratory of Geotechnics, Department of Civil Engineering University of Salerno (UniSa) Via Giovanni Paolo II, 84084, Fisciano (Salerno), Italy e-mail: [email protected], [email protected], [email protected] † Department of Applied Mathematics, ETS de Ingenieros de Caminos, Universidad Politecnica de Madrid (UPM), Calle del Profesor Aranguren, 28040, Madrid, Spain e-mail: [email protected] Key words: earthflow, modelling, mobilization, propagation, reactivation Abstract. Natural slopes in clayey soils are often affected by failures which may cause the onset of landslides of the flow type travelling large distances and damaging buildings and major infrastructures. Particularly, the so-called earthflows pose challenging tasks for the individuation and forecasting of the remobilized masses; as a consequence, the mathematical modelling of the propagation stage allows enhancing the understanding of earthflows in order to obtain reliable assessments of run-out distances and displaced soil volumes. This paper deals with the reactivations of Montaguto earthflow (Southern Italy) occurred from 1998 to 2009 that are simulated, through the depth-integrated “GeoFlow-SPH” model, thanks to the availability of a detailed data-set. The achieved results provide a satisfactory agreement with the in-situ information and outline how a change of the rheology of the mobilized masses can affect the whole phenomenon. 1 INTRODUCTION Earthflows [1] are a common type of mass movement in mountainous regions with finegrained soil or very weathered bedrock and span from small events of 100 m2 in size to large events encompassing several km2 [1], [2], [3]. Most rapid earth flows [4] occur in areas of highly sensitive clays with low plasticity due to salt loss from pore water; typical cases are documented in Scandinavia, eastern Canada, Alaska, Japan, the former Soviet Union and New Zealand. Particularly, once the initial slide begins moving, the mobilized material often liquefies and begins to flow. The entire failure process is usually completed within several minutes to several hours. High sensitivity and low post-failure shear strength of quick clay often results in very long run-out of rapid earth flows. SPH propagation modelling of an earthflow from Southern Italy 522 L. Cascini, S. Cuomo, M. Pastor and M. Coppola 2 The deposits of earthflow are hummocky, several meters to a few tens of meters lower than the original ground surface, and they slope gently towards the stream channel. The postfailure rheological response of quick clays is controlled by complex interactions among pore water chemistry, clay mineralogy, and grain size distributions. On the other hand, slow earthflows occur in plastic silts or clays, as well as rocky soils that are supported by a plastic silt-clay matrix. Slow earthflows commonly have a teardrop or bulbous, a sinusoidal profile, are elongate in the direction of down slope movement and are several times wider than thick. Coe et al. [4] show that the variations of the hydrologic boundary conditions at the ground surface are closely correlated in time with accelerations of an active earthflow. Once reactivated, the earthflow can attain quite a high velocity (up to meters per hour and more) moving as a flow; later on, its velocity progressively decreases to a complete stop that can occur even tens of years after reactivation, unless a new trigger provokes acceleration. Mathematical modelling of the propagation stage contribute to enhance the understanding of the complex behaviour of an earthflow and also can allow reducing losses inferred by this type of phenomena, as it provides a means for individuating the hazardous areas and defining the best mitigation measures. This paper deals with Montaguto earthflow (Southern Italy) reactivated four times from 1998 to 2009 as outlined by Cascini and Di Nocera [3] During the last two reactivations a strategic transport corridor located at the toe of the hillslope was allowed to be interrupted for several weeks, thus requiring urgent remedial works. On the basis of the advanced geological and geotechnical knowledge acquired by Cascini and Di Nocera [3], this paper investigates the 2006 earthflow reactivation which was characterised by a long run-out distance; based on numerical modelling, some indications are provided for rheological behaviour of mobilised volumes which can be useful for the analysis of future reactivations. 2. MONTAGUTO EARTHFLOW Montaguto landslide is located in the Daunian Apennine (southern Italy) and is one of the largest active earthflows in Europe [3], [5], [6]. The landslide developed in the Cervaro valley following the local morphology and spanning a total length of almost 3 km, with an elevation drop of about 430 m (Fig. 1a). The toe zone of the earthflow has been affecting the National SS90 road and the Roma - Bari railway which are two key corridors of the E-W transport system in Italy. For this reason, to date several studies were carried out to characterize the whole affected area [3], [7] and the main results are hereafter summarized and will be later used as input data for geomechanical modelling of the propagation stage. In the Montaguto earthflow area (Cerrato valley), the geological setting is quite complex with outcropping plastic clays, silty clays, sandy marls, marly calcareous and clayey flysch (Fig. 1b). Slope elevations span from 956 m a.s.l. at “la Montagna” ridge to 401 m a.s.l. in the valley. In the upper zone, the slope aspect is mainly W-SW with the drainage network and secondary ridges SW-NE exposed; in the medium-low portion of the hillslope NNW-SSE is the main exposure. Slope angles are lower than 15° in the valley where the earthflow moves while the adjacent zones are 10° to 30° steep. Due to the main morphological elements whose details are provided by Cascini and Di 523 L. Cascini, S. Cuomo, M. Pastor and M. Coppola 3 Nocera [3], the landslide body is displaced at rates equal to meters/month, while local landslide lobes move at velocities of meters/week or meters/day, respectively, due to retrogression or progression phenomena. Figure 1: Montaguto earthflow: a) ortophoto dated on 2006 (source: Avioriprese); b) geological setting, legend: landslide body (a1); alluvial deposit (b); colluvial deposit (b7); Altavilla Geological Unit: sands (UTAa) and silts (UTAb); Toppo Capuana clayley marls (TPC); Faeto Flysch (FAE). (modified by [6]) Montaguto earthflow suffered four main reactivations between 1998 and 2009 [3], which were strongly related to the superficial water circulation and can be modelled thanks to a comprehensive data-set well described in Cascini and Di Nocera [3] and including: i) several historical and recent topographic maps from 1876 to 2009, ii) orthophotos available for the whole area dated on 1998, 2006 and 2009 (data from Avioriprese), ii) in-situ investigations and geotechnical slope stability analyses iii) satellite images interpretation [7]. Figure 2 shows the propagation path of the earthflow on May 2006 [3] and a reconstruction of the mobilized/deposited soil heights based on 10mx10m Digital Terrain Models (DTMs) dated on 2006 and 1998. In this period, about 10-20 meters of soils were mobilized at the crone zone while deposition heights of 5-10 meters were recorded in the central part of the earthflow (at the neck of propagation path) and deposits of 10-20 meters at the terminal part of the earthflow. The complex dynamics of the earthflow is also evidenced by Casagli et al. [7] who measured - based on SAR interferometry technique - a maximum velocity of 2.9 m/day on June 1st, 2010 at the lower part of the earthflow and provide thermal images which indicate the wet areas corresponding to the highest velocities of the mobilized masses. 524 L. Cascini, S. Cuomo, M. Pastor and M. Coppola 4 a) Difference (m) b) Figure 2: Montaguto earthfow: a) photo of Montaguto earthflow (dated on May 2006), b) Difference between 10x10m DTMs obtained from the topographic maps dated on 1998 and 2006. Negative values correspond to soil mobilization while positive values to deposition (modified by [6]). 3. SPH PROPAGATION MODELLING 3.1. Methods and inputs The “GeoFlow-SPH” model [8], [9] is used to simulate the propagation stage of Montaguto earthflow during one of the last major reactivations. The mathematical model is based on the theoretical framework of Hutchinson [10] and Pastor et al. [11] and schematizes the propagating mass as a mixture of a solid skeleton saturated with water. The governing equations are: i) balance of mass of the mixture combined with the balance of linear momentum of the pore fluid, ii) balance of linear momentum of the mixture, iii) rheological equation relating soil stress tensor to deformation rate tensor and iv) kinematical relations between deformation rate tensor and velocity field. In the case of earthflows, similarly to many other flow-like landslides [8], [9], average depths are small in comparison with their length or width and thereby the governing equations can be integrated along the vertical axis and the resulting 2D depth integrated model presents an excellent combination of accuracy and simplicity. In the “GeoFlow-SPH” model, the Smoothed Particle Hydrodynamics (SPH) method is used which discretizes the propagating mass through a set of moving “particles” or “nodes”. Information, i.e. unknowns and their derivatives, is linked to the particles and the SPH discretization consists on a set of ordinary differential equations whose details are provided by Pastor et al. [8]. The accuracy of the numerical solution and the level of approximation for 525 L. Cascini, S. Cuomo, M. Pastor and M. Coppola 5 engineering purposes depend on how the nodes are spaced, as recently reviewed by Pastor and Crosta [12] and Cuomo et al. [13]. Furthermore, the detail of the available DTM is a crucial issue for a proper description of the ground surface and it entails the global reliability of the numerical results [11]. Based on the description provided by Cascini and Di Nocera [3], four main zones (Fig. 3a) of Montaguto earthflow can be depicted. Zone 1 is the source area of the 2006 reactivation, zone 2 is a deposition area located at the middle part of Montaguto earthflow. In the zone 3 – which includes zone 2 – either propagation of unstable volumes takes place or remobilization of material previously deposited, as occurred during the 2006 reactivation. Finally, zone 4 is the terminal propagation/deposition zone of the whole earthflow (Fig. 3a). Two series of numerical analysis are performed to simulate: i) the propagation down slope the zone 1 occurred during the 1998-2006 reactivation, ii) the remobilization of material from zones 2 and 3 towards the toe of the slope. As input data, a 10mx10m Digital Terrain Model (DTM) is obtained from the topographic information available for 1998 and 2006 (data from [3]), thus defining the source area as well as the lateral and basal limits of the slip surfaces on which the flow moves. It is worth noting that for the upper part of the hillslope (source area) the 2006 DEM was used (post-failure information), while the 1998 DTM was used (pre-failure information) for the lower part of hillslope (depositional area) where the earthflow had not yet arrived. The topography so elaborated is shown in Figure 3b. As for the 2006 reactivation, the total triggered volume is individuated by difference of the digital elevation models available for 2006 and 1998 [3], [5]. An area of about 180,000 m2 is estimated inside the zone 1 with a maximum triggered soil height equal to about 27 m at the upper portions of the hillslope; this volume is used as input data for propagation analysis (Fig 4a). Analogously, it is estimated that the 2006 reactivation almost mobilised the zones 2 and/or zone 3 (whose global extent is about 105 m2) with initial soil heights equal to about 4-5 m (Fig. 4b, 4c). The objective of SPH analyses is twofold: i) investigate the potential of “GeoFlow-SPH” model to adequately simulate this landslide type, and ii) investigate the most adequate rheological parameters to capture the propagation behaviour of the earthflow during one of the last major reactivations (period 1998-2006). To this aim, a simple viscous-type a) b) Figure 3: a) Main zones of Montaguto earthflow, b) 3D-view of the topography used for modelling. 526 L. Cascini, S. Cuomo, M. Pastor and M. Coppola 6 rheological model such as Bingham law is used since it proved as much as effective to backanalyse a wide class of case histories [15]; particularly, the scientific literature outline that mud/debris materials are characterised by a Bingham yield stress (  y) equal to 102÷103 Pa and viscosity (  ) spanning from 20 to 800 Pa·s [15]. Alternatively, the mass could be schematized as mixture of solid particles with a frictional rheology and pore water with pressures changing in time and space due to consolidation [8]; however, such analysis is beyond the scope of this paper. The list of the most significant numerical cases is reported in table 1. It is worth noting that the use of 10 m DTMs has the major advantage to allow the investigation a large series of cases in a reasonable time while has the drawback to ensure a global description of the ground surface which not consider site-specific details [12], [13]. However, this approximation is acceptable for gentle slopes and also consistent with the availability of DTMs which are some years time-spaced each other. a) b) c) P Figure 4: Initial soil heights for numerical simulation of the earthflow reactivation: a) stage 1 and b-c) stage 2. Table 1: List of the main numerical simulations. Stage ID* [Pa.s] [Pa] source zone triggered heights of soil (m) Stage 1 Stg1_4 50 17 zone 1 27 Stg1_3 500 170 Stg1_2 1000 1000 Stg1_6 1000 3000 Stg1_7 1000 5000 Stg1_8 1700 5000 Stg1_1 2000 1700 Stg1_5 2000 2500 Stg1_10 5000 1500 Stg1_9 6000 1700 Stage 2 Stg2_5 500 100 zone 2 5 Stg2_6 500 50 Stg2_11 300 50 Stg2_7 100 50 Stage 2 Stg2_0 1000 500 zone 3 4 Stg2_1 500 100 Stg2_2 500 50 Stg2_3 300 100 * soil unit weight  = 1600 kg/m3 527 L. Cascini, S. Cuomo, M. Pastor and M. Coppola 7 3.2. Numerical results A first set of analyses is aimed at simulating the reactivation of the whole earthflow (zones 1-4), also providing some reference values for Bingham rheological parameters as it concerns the upper part of the earthflow (zone 1). Considering that the whole earthflow travelled a very long run-out distance (about 3 km), the lowest values of rheological parameters provided by the scientific literature are firstly referred [15]. However, most of the numerical tests (Stg1_1 – Stg1_8) are unable to provide a reliable description of the earthflow, because the masses mobilised in the zone 1, after a travelled distance of about 1 km in the “Fosso Nocelle” torrent, fall into the “Tre Confini” valley at the right-hand side of the observed propagation path (Figure 5). Figure 6a shows the simulated heights at point P of figure 4a (in the zone 2) and similar time trends can be noted with different peak and final values. Particularly, at point P, the initial soil volume is mobilized due to gravity and soil heights decrease from about 5.5 m to less than 1÷3 meters in few tens of seconds. Then, a surge of material is recorded with a 2÷9 m peak height. Finally, propagation (Stg1_3, Stg1_4) or deposition (Stg1_6 – Stg1_8) can occur at point P depending on mass rheology. The final height looks much more affected by viscosity () than yield stress (y), as the comparison of cases Stg1_6, Stg1_7 and Stg1_8 outlines. It is worth noting that: i) all the mentioned cases refer to both rheological parameters larger than those from literature, ii) a very large value of yield stress (y) – up to 5000 Pa – does not afford a shear strength sufficient to avoid the material falling into the “Tre Confini” valley. Further analyses (Stg1_9 and Stg1_10) – with a high value of viscosity () and yield stress (y) taken from literature – provide different scenarios which is straightforward discussing (Fig. 6b). Viscosity () regulates the whole duration of the propagation/deposition stage, also modifying the run-out distance and propagation pattern in time and space (Fig. 5e). Particularly, an increasing value of viscosity () and a reduction of yield stress (y) correspond to the propagation of more surges of material at the monitored point P (see test Stg1_9 compared to Stg1_1 – Stg1_8 in figure 6b, as also outlined by arrows). This behaviour, in turn, corresponds to velocity fluctuations of the whole earthflow, as outlined by Leroueil [16] for an active eartflow. Both tests Stg1_9 and Stg1_10 provide deposition heights equal to 5 ÷ 6 m in the zone 2 which well agree those obtained from topographic maps dated on 1998 and 2006 (Fig. 2b). An overview of the best fitcase (Stg1_9) is provided in figure 7 which evidences a proper description of the observed earthflow through the selected rheological parameters. It is also outlined that the earthflow cannot be simply schematized as a unique process which triggers a mass in the zone 1 which later propagates up to the distal part of the zone 4. This is confirmed by further propagation analyses performed referring to the rheology of Stg1_10 and considering as source areas: i) zone 1 and 2, ii) zone 1 and 3. In both cases, some mass falls in the “Tre Confini” valley and some mass stops before reaching the distal part of zone 4. Therefore, multiple stages of the earthflow must be simulated referring to distinct source areas in time. A second set of analyses is aimed at investigating the remobilization of earthflow from the middle (zone 2 or 3) to the terminal part (zone 4) of the path. If the zone 2 is considered as remobilization area (with initial soil heights equal to about 5 meters) an overall mismatch is obtained between of the numerical results and the observed earthflow, for any combination of 528 L. Cascini, S. Cuomo, M. Pastor and M. Coppola 8 a) b) c) d) e) Figure 5: Planview of the initial (a) and final heights of soil for the cases: b) Stg1_4, c) Stg1_6, d) Stg1_8, e) Stg1_9 of Tab. 1. 0 1 2 3 4 5 6 7 8 9 10 0100 200 300 400 500 600 700 800 900 1000 soil height (m) time (s) Stg1_1 Stg1_5 Stg1_9 Stg1_10 9 5 1 10  0 , 0 b) 1 2 3 4 0 1 2 3 4 5 6 7 8 9 10 0100 200 300 400 500 soil height (m) time (s) Stg1_4 Stg1_3 Stg1_6 Stg1_7 Stg1_8 3 6 7 8 4  a) Figure 6: Simulated heights of soil at point “P” compared to in-situ observations (4-6m tick deposit) a) b) c) d) Figure 7: Some simulated steps (a-d) of the best fitcase (test Stg1_9 of Tab. 1) for the stage 1of the earthflow reactivation in the period 1998-2006. 529 L. Cascini, S. Cuomo, M. Pastor and M. Coppola 9 rheological parameters among those listed for tests Stg1_1 – Stg1_10 of table 1. With different rheological parameters (viscosity  equal to 100÷500 and yield stress y equal to 50÷100, tests Stg2_5-6-7-11 of Tab. 1) the numerical modelling is slightly better but still poor (Fig. 8). It is worth noting that rheological parameters either higher or similar to those from literature allow a satisfactory simulation of the stage 1 of the earthflow. It entails that the remobilization did not occurred in the zone 2 during the second major stage of the earthflow. Aimed at simulating the second stage of the analysed earthflow, the zone 3 was assumed as the whole remobilization zone in a set of analyses with different rheological parameters listed in table 1. The obtained results (Fig. 9) outline the possibility to properly describe the earthflow, provided an adequate calibration of the rheology. Particularly, the results of the best fitcase are shown in figure 10 and it can be noted that: i) the maximum run-out distance is well simulated, ii) different surges of material are simulated as typically occurs for earthflow, iii) a partial mobilization of the initial heights is correctly simulated in zone 3. a) b) c) d) e) Figure 8: Planview of the (a) initial height of soil in the zone 2 and final simulated height of soil for case: b) Stg2_5, c) Stg2_6, d) Stg2_7, e) Stg2_11 of table 1. a) b) c) d) e) Figure 9: Planview of the (a) initial height of soil in the zones 2 and 3 and final simulated height of soil for case: b) Stg2_0, c) Stg2_g1, d) Stg2_2, e) Stg2_3 of table 1. 530