Assessing flank instability of Stromboli volcano (Italy) by reappraising the 30 December 2002 tsunamigenic landslides

Bearing in mind the destructive potential of tsunamis induced by volcanic landslides, the tsunamigenic event occurring at Stromboli volcano in Italy on 30 December 2002 has been re-examined here, by means of visible images and slope stability analysis. This was one of the few examples in the world of a flank collapse occurring at a volcano that was directly observed. We present the results of stability analyses, together with a sequence of photos collected from a helicopter a few minutes before the collapse. The result of this study is that the sequence of landslides triggering the 2002 Stromboli tsunami can be defined as the final stage of a lateral magma intrusion that exerted a high thrust at high altitude, destabilizing the entire slope. This study allows a more complete understanding of the event that took place on Stromboli on 30 December 2002. Furthermore, the approach used here, if appropriately modified, can be used in other contexts, contributing to the understanding of the condition that leads to tsunamigenic landslides.


Introduction
The instability of volcanic edifices is the result of the complex and mutual relationship between endogenous and exogenous processes (Gorshkov 1959;Francis 1997, 2001;Reid et al. 2001). Volcano slope instability can produce destructive mass-wasting (Voight et al. 1981;Siebert 1984Siebert , 1996Voight and Elsworth 1997), which in turn, in the case of insular or coastal volcanoes, can generate tsunamis (Ward and Day 2003). In the last centuries, several episodes of tsunamis generated by volcanic landslides have been documented: Oshima-Oshima volcano (Japan) in 1741 (Ioki et al. 2019), Unzen volcano (Japan) in 1792 (Sassa et al. 2016), Mt. Augustine (Alaska, USA) in 1883 (Kienle et al. 1987), Ritter Island (Papua New Guinea) in 1888 (Karstens et al. 2020) and Anak Krakatau in 2018(Walter et al. 2019. Here, data on the tsunamigenic landslide sequence on 30 December 2002 on the Stromboli volcano in Italy (Bonaccorso et al. 2003;Tinti et al. 2006;Chiocci et al. 2008;Tommasi et al. 2005Tommasi et al. , 2008 has been re-examined. In addition to reviewing the data collected before the collapse and already published (Bonaccorso et al. 2003;Calvari et al. 2005;Burton et al. 2008;Pioli et al. 2008;Federico et al. 2008), the stability of the NW slope of Stromboli (the socalled Sciara del Fuoco, SdF) was re-analysed considering all the geomechanical/geotechnical data on the SdF infilling (Apuani et al. 2005a, b;Tommasi et al. 2005;Boldini et al. 2005Boldini et al. , 2009Apuani and Corazzato 2009;Rotonda et al. 2010;Verrucci et al. 2019). A reliable stratigraphy is proposed to model the sequence of destabilization phases of the SdF flank. In addition, a series of photos from helicopter surveys collected just before the collapse (Calvari et al. 2005) were analysed. This unique set of data shows the propagation of the fractures along the slope a few minutes before the slope failure and improves our interpretation of the causes triggering the instability. This study can shed light on both the triggering mechanisms of the tsunamigenic phenomena of Stromboli and make a general contribution to the understanding of volcanic edifice failures that can trigger tsunamis.

Stromboli volcano
Stromboli volcano (Aeolian Archipelago) is located in the Tyrrhenian Sea (Central Mediterranean Sea; Fig. 1a). The volcano extends from 2200 m b.s.l. up to 926 m a.s.l. and is renowned for its persistent Strombolian explosive activity (Blackburn et al. 1976), produced from a series of vents located within a crater terrace at 750 m a.s.l. The persistent activity consists of small-scale explosions launching pyroclasts up to a few hundred metres high (Patrick et al. 2007). Sometimes, larger explosions occur (Barberi et al. 1993;Bevilacqua et al. 2020;Calvari et al. 2021). Explosive activity is sometimes accompanied or replaced by effusive activity (Calvari et al. 2005(Calvari et al. , 2010(Calvari et al. , 2014(Calvari et al. , 2020Di Traglia et al. 2014, 2018a. At Stromboli, flank eruptions are described as effusive eruptions with vents located outside the crater terrace, generally within the SdF flank (see Marsella et al. 2012). The volcano has experienced repeated lateral collapses, both on the NW slope (SdF) and on the SE flank (Romagnoli et al. 2009;Francalanci et al. 2013). In particular, the SdF has a slope ranging from 40 to 45° in the subaerial part, remaining > 30° down to 300-400 m b.s.l. (Kokelaar and Romagnoli 1995), and is filled with volcaniclastic deposits and lavas (Kokelaar and Romagnoli 1995;Casalbore et al. 2010Casalbore et al. , 2011Di Traglia et al. 2018b).
Current landslide phenomena at Stromboli vary from rockfall and gravel flows (i.e. < 10 6 m 3 ) to rockslides/rock avalanches (i.e. 10 6 -10 7 m 3 ). The significant steepness of the slope, combined with the considerable heterogeneity of the infilling material and the constant emplacement of volcaniclastic material due to explosive and effusive activity, are all factors that promote the persistent small-to medium-scale mobility of the flank. Larger events include deepseated gravitational slope deformations that evolve into debris/rock avalanches (i.e. > 10 7 m 3 ), typical of sector collapse events (Schaefer et al. 2019). The most recent landslide of notable importance (i.e. 25-30 × 10 6 m 3 ) Chiocci et al. 2008Chiocci et al. ) (2023Chiocci et al. ) 20:1363Chiocci et al. -1380 occurred on 30 December 2002 (Bonaccorso et al. 2003;Acocella et al. 2006) and caused a tsunami sequence that impacted the northern coast of the Stromboli village (maximum run-up heights of about 11 m; Tinti et al. 2005Tinti et al. , 2006. Tsunamis in Stromboli are quite frequent (on average 1 every 20 years), as evidenced by recent events such as in 1879, 1916, 1919, 1930, 1944and 1954(Maramai et al. 2005. Tsunamis in Stromboli mainly impact the northern coastal area of the island (Turchi et al. 2022) but can also reach the Sicilian and Calabrian coasts ~ 60 km away (Tinti et al. 2006). Moreover, tsunami deposits were identified on the north-eastern coast of the island and would be connected to events that took place during the Middle Ages (Rosi et al. 2019;Pistolesi et al. 2020). Extreme run-up deposits probably connected to the Sciara del Fuoco flank collapse that presumably occurred 5 kyr ago have been found on the southeast flank of the island at about 120-m elevation (Tanner and Calvari 2004). Observed tsunami events have always been connected to intense volcano activity and to (i) submarine landslides which, in turn, are triggered by flank destabilization induced by magma intrusions or material accumulations (as in 2002; Acocella et al. 2006) or (ii) pyroclastic density currents (PDCs) produced by paroxysmal explosions (as in 1930;Di Roberto et al. 2014).

The 30 December 2002 landslides
The most significant event of slope instability to occur on Stromboli in recent times was the tsunamigenic landslide sequence 2 days after the onset of the eruption (30 December 2002Bonaccorso et al. 2003;Calvari et al. 2005). The 2002-2003 flank eruption was anticipated by an intensification of the eruptive activity from May 2002, with increasing intensity and frequency of the explosive activity, and lava overflows from the summit craters, due to the rising magma within the shallow conduit (Burton et al. 2008). This activity was accompanied by high soil CO 2 fluxes and Rn concentrations in the summit detected about 2 weeks before the eruption, testifying to the increase in the bulk (gas + magma) input rate from the deep storage towards the surface (Burton et al. 2008;Federico et al. 2008). In the afternoon of 28 December 2002, the opening of a NE-SW trending, 300-m-long eruptive fissure below the NE crater (NEC) occurred (a in Fig. 2a; Calvari et al. 2005;Acocella et al. 2006). A hot avalanche, comprising a mixture of lava overflow and loose volcaniclastic deposits from the NEC, occurred on the same day at 17:30 UTC (Calvari et al. 2005;Pioli et al. 2008). After a 2-h eruptive pause, a new lava flow erupted from the lower end of the eruptive fissure and reached the sea. A helicopter survey carried out the morning of the next day (29 December 2002) using a thermal camera revealed the inactive lava flows and a number of transversals, horizontal cracks developing on their surface and later interpreted as the first signal of slope instability (Bonaccorso et al. 2003;Calvari et al. 2005). New lava flows erupted on the same day after about a 12-h pause from vents at 600-and 500-m elevation along the SdF (b and c in Fig. 2a), lasting again a few hours and reaching the sea (Calvari et al. 2005). A new helicopter survey (30 December 2002 morning) enabled detecting the development of a deep-seated slope deformation called α slide (Fig. 2a, b) affecting the central-northern SdF  and the formation of several fractures developing just before the failure. Another helicopter survey carried out in the afternoon, after the tsunami, revealed two large scars in the middle area of the SdF (Bonaccorso et al. 2003). The area affected by the α movement experienced failures and large deformation between 600 and 700 m a.s.l. and the establishment of longitudinal shear bands probably associated with the development of shear surfaces dipping seawards . The occurrence of a submarine failure (τ1 slide, Fig. 2a, b), from the coast to 340 m b.s.l., with a sub-circular slip surface, removed ~ 9.5 × 10 6 m 3 of material . Retrogressive submarine-subaerial collapses (τ2 slide, β slide and γ slide, Fig. 2a Tommasi et al. 2008), accounted for 11.6 × 10 6 m 3 of material removed, with a total mass involved in the failure of 21.1 × 10 6 m 3 .
Based on broad-band seismic signals, it was possible to affirm that the submarine failure detached on 30 December 2002 at Grande-Le Schicciole), the direction of diking since 100 ka (modified after Tibaldi et al. 2009) and the location of the Crater Terrace, where the present-day persistent activity occurs 12:14:05 UTC, whereas the subaerial failures started to occur at 12:22:38 UTC (Bonaccorso et al. 2003). Seismic signals associated with the tsunami wave hitting the northern side of Panarea island (20 km SW of Stromboli) were recorded at ~ 12:20 UTC and ~ 12:27 UTC, whereas tide-gauge recorded the tsunami arrival at the Panarea harbour between 12:19 and 12:24 UTC (Pino et al. 2004;La Rocca et al. 2004). The maximum wave run-up was just under 10 m high on the northern Stromboli coastline, with a maximum inundation extension of 100 m inland on the NE part of Stromboli, causing significant damage to buildings but no casualties (Tinti et al. 2006).
Understanding what type of instability action triggered the α movement is a starting point to explain the sequence of the 30 December 2002 landslides. The hypothesis that seismic acceleration initiated the deformation process in the upper part of the slope was discarded, firstly because of the small amplitude of the seismograms recorded during those days, secondly because during the 5 April 2003 paroxysm (Calvari et al. 2006), albeit the powerful explosion occurred with associated large deformation (Bonforte et al. 2008), the slope was not affected by any failure. This was also confirmed during the strong explosions taking place in 2007 and 2019 (Ripepe et al. 2021). Various authors proposed that the α movement was probably caused by the thrust exerted by a dyke (Acocella et al. 2006) intruding between the volcaniclastic deposits and a stronger SdF substratum (Apuani et al. 2005b;Tommasi et al. 2005). The lateral thrust exerted by the magma ascending along the dyke would have caused the plasticization of the summit area and mobilized the shear strength along the SdF deposits. According to previous studies Boldini et al. 2009), the progressive shear deformations caused by the α movement activated grain crushing within the shear zones, enabling the transition to an undrained condition in the submarine portion of the slope. The increase in pore pressure promoted the draw-down of shear strength, in a process called static liquefaction. The removal of submarine material deprived the terrestrial slope of its support, leading to the subaerial failures β and γ (Fig. 2a, b).

Photographic surveys
The 30 December 2002 helicopter survey began at Stromboli at 11:19:24 (UTC time) and lasted until 11:43:17, for a total flight of 23 min and 53 s. During the helicopter flight, we collected a total of 266 visible photos, taken at a distance of flight path from the slope between about 800 m and 400 m. The photos were recorded at a frequency of one image every 2-4 s and using a digital camera. We used a CANON A40 digital camera with lens of 50° furnishing images of 2272 × 1704 pixels. Unfortunately, the camera clock was not calibrated before the flight; thus, we lack a precise correspondence between photos and other geophysical data.

Limit equilibrium analysis
The preliminary approach to estimate the slope stability in active volcanic environment was proposed by Voight et al. (1981) who produced an analysis of the 1980 lateral collapse that affected Mount St. Helens through 2D limit equilibrium method (LEM) of slices. Voight and Elsworth (1997) simulated the destabilizing action of magma intrusion, hydrothermal activity and tectonic seismicity through a 3D LEM analysis performed on a simplified model of the volcanic slope. Over the years, several authors adopted LEM analyses to assess instability phenomena affecting volcanic flanks (Reid et al. 2001;Simmons et al. 2005;Battaglia et al. 2011;Borselli et al. 2011Borselli et al. , 2016Dondin et al. 2017;Gonzalez-Santana and Wauthier 2021;Di Traglia et al. 2022), including the SdF depression (Apuani et al. 2005a;Tommasi et al. 2008;Di Traglia et al. 2018c;Schaefer et al. 2019). In this work, the slope stability at Stromboli is evaluated with the Slope Stability Analysis Program (SSAP; Borselli et al. 2011Borselli et al. , 2020, a free software that investigates the factor of safety (FS) of sliding surfaces through bi-dimensional LEM analysis by assigning the shear strength parameters for different failure criterion as: linear Mohr-Coulomb; generalized non-linear Hoek and Brwn (Hoek et al. 2002;Hoek and Brown 2019), non-linear Barton-Bandis (Barton 2013). SSAP also allows setting the water table position, additional loading conditions, already used to estimate the stability of volcanic edifices in different geological contexts (Borselli et al. 2011;Dondin et al. 2017;Di Traglia et al. 2022). In this study, the factor of safety is estimated through the Morgenstern-Price method (Morgenstern and Price 1965). The search of critical sliding surfaces is carried out through the "Random Search" engine that generates possible sliding surfaces (sub-planar, convex, concave-convex, composite) through the Monte Carlo methods (Siegel et al. 1981). SSAP generally adopts a feature called "Dynamic Surface Search Attractor" that progressively reduces the initial search area set by the user or provided automatically by the software, according to the surfaces with lower FS gradually and progressively identified (Borselli 2020). This function is extremely useful because it allows the user to find ever more critical surfaces within a specific area. At the same time, the software tends to discard potentially unstable areas if not identified in the early stage of the research. SSAP gives as output a graphical representation of the surfaces with lowest FS, internal distribution of forces and pressures along the sliding surface and in the internal mass (e.g. pore pressure, effective normal stress, tangential and horizontal and inter-wedge forces), distribution of local FS values along the sliding surface, local distribution of the numerical reliability of general FS numerical solution (namely RHO index (Zhu et al. 2003;Borselli 2020)). At the same time of the FS obtained by LEM, SSAP produces a 2D colour map with distribution of average local FS acquired by local stress distribution obtained by the large set of Monte Carlo generated sliding surfaces. According to the SSAP documentation (Borselli 2020), three types of 2D colour maps can be produced, defined as quasi-finite element analysis (qFEM) (Borselli 2020) derived by classical approach from Schofield and Wroth (1968), Griffiths and Lane (1999) and over stress ratio (OSR) maps (Farias and Naylor 1998). These maps, based on a combination of LEM and FEM procedures, can highlight the presence of local areas within the slope that may display a critical state of stress (high OSR and a low FS). Although the actual state of deformation is not calculated by the software (SSAP does not allow the user to include the strain parameters in the slope characterization), this feature provides a graphical representation of those areas in which progressive failure can originate. Based on Schofield and Wroth (1968) and Griffiths and Lane (1999), the qFEM map describes the local distribution of the FS, including recording the local slope of each portion of all generate surfaces, a representation of the potential direction of possible plasticization for those areas where local FS < 1.0. A local OSR map, based on Farias and Naylor (1998), displays in terms of local mean stress (principal stress and path stress) the areas where the maximum local shear stress is greater than the local shear strength. The areas with OSR > 1.0 are the most likely to develop progressive failure.
The slope's geometry is obtained by digitalizing the contour lines included in Tommasi et al. (2008) and Baldi et al. (2008). Then a slope profile line that passes through the area mostly affected by failures and deformation during the 30 December 2002 landslide sequence was selected for the LEM analysis (maximum number of 100 nodes imposed by the SSAP). In each simulation performed in static and seismic conditions, the water table is assumed to be at the sea level (above it the rocks are in dry condition). This assumption is also supported by meteorological and geophysical observations reported in Revil et al. (2011).

Input parameters and material constitutive models (strength criteria)
Although the volcaniclastic material comprising the slope is an alternation of layers with different grain size and lithology, for modelling purposes it is necessary to define layers with homogeneous physical and geotechnical properties. We assumed that in each layer there are coexisting, at the same time, two different systems: (1) a rock mass fractured system, where the factures and discontinuity are randomly distributed and (2) a system where fractures or discontinuities in stratigraphic alternation or joints are structurally or stratigraphically defined. This coexistence exists in stratovolcano edifices where we can observe rhythmic alternation in volcaniclastic deposits and lava deposits. The interface of volcaniclastic deposits/lava may be a zone where we can consider a condition of potential discontinuity that is preferentially oriented and where we have a potential zone of weakness. This stratigraphic condition requires a change in the classic way we use to parametrize the constitutive volcanic origin geomaterials. In particular, we need to change the way to model local shear strength and account for possible variability in shear strength and in order to search for a potential zone of instability and weakness in this type of peculiar rock mass. The potential zone of weakness at the interface between volcaniclastic deposits and lava deposits is analogue to the interface between bedrock and rockfill introduced and modelized by Barton (2013).
In this case, the SSAP software allows considering for each strata a combination of two non-linear failure criteria that are selected applied in different portions along the analysed sliding surface. The choice of the appropriate failure criteria is done, at run-time, depending on the local slope along of the assumed sliding surface (Lunardi et al. 1994a, b). Some slope units have assigned at the same time two non-linear criteria. In particular, it is assigned the generalized Hoek and Brown criterion (GHB) (Hoek et al. 2002;Hoek and Brown 2019), implemented in SSAP, following Carranza-Torres (2004) and Lei et al. (2016), and the non-linear failure envelope by Barton-Kjaernsli (1981) and Barton (2013) (B-K criterion) for the rockfill-like material in its local discontinuity with lava local substratum.
SSAP software according to the Lunardi et al. (1994a, b) methodology selects t to assign the B-K criteria where the local slope angle of sliding surface coincides with the slope angle of predefined discontinuity existing at alternation between rockfill-like material and lavas. Figure 3 shows a scheme of the SdF, the association of volcaniclastic (rockfill material) and lava-breccia unit (alternated by lava and breccia layers). The mean slope angle (β) influences the mean slope deposition surfaces and interface of alternate volcaniclastic and lava deposits. We propose an ubiquitous shear strength combined criterion that associates the B-K criterion, in case of local collinearity of sliding surface portions with β + Δβ angle, and GHB criteria where the sliding surfaces are outside this interval of slope. The anisotropy exists locally depending on the local slope of considered failure surface. The GHB criteria are not assigned to a preferential set of discontinuities because it is applied at any type of slope of failure surfaces. B-K instead can be applied only to specific families of discontinuities. The GHB must be associated to a rock mass random set of fractures (e.g. fractured lava rock, mass and volcanoclastic rockfill), instead B-K criteria can be associated to the existing discontinuities between lavas and volcanoclastic material. Because these materials are in alternation on the volcano slope and cannot be done a clear separation between them, we assume the coexistence of both strength criteria (GHB, B-K). The selection of strength criteria to be used is made at runtime by the software SSAP, depending on the local geomaterial properties of the sliding surface assumed.
Given the extremely large range of confining pressures considered, the interpolation of HB parameters in MC criterion led to unrealistically high values of cohesion at shallow depth, as demonstrated by Li et al. (2008) and Shen et al. (2012). For this reason, generalized Hoek-Brown criterion (GHB) should always be preferred. In literature can be find the application of non-linear HB failure criterion to non-rock mass geomaterials as overconsolidated clays (Li and Yang 2019). At the same time, we encounter that rockfill must be modelled in terms of non-linear failure criteria (Charles and Watts 1980;Indraratna et al. 1998;Liu 2009;Frossard et al. 2012;Barton 2008Barton , 2013. The SdF behaviour in response to both static (e.g. not seismic and not dyke intrusion) and loading (seismic and/or dyke intrusion pressure) conditions has been reproduced considering a set of geotechnical unit stratigraphy as follows: -A shallow layer in the subaerial part of the slope made up of volcaniclastic material Boldini et al. 2009;Rotonda et al. 2010) that can be classified as rockfill-like material. The GHB failure and the non-linear failure's envelope by Barton-Kjaernsli (1981) are considered for this layer. -Lava-breccia unit (Apuani et al. 2005a, b), described as lava layers (between 35 and 65%) alternated with breccia material, that accurately and efficiently describe the behaviour of the volcanic edifice both at shallow and higher depths. This lithostratigraphic unit was proposed for rock masses with analogous percentage of lava and breccia components, namely as a lower boundary in the geotechnical characterization of the volcanic structure. The GHB failure and the non-linear failure's envelope by Barton-Kjaernsli (1981) are considered for this layer. -A stronger substratum mainly composed of the lava unit (i.e. percentage of lava layers greater than 65%, Apuani et al. 2005a). The Mohr-Coulomb failure criterion is considered for this deeper zone.
The D values are the most problematic parameters in modelling the rock mass in continuous domain following the GHB failure criterion. Hoek and Brown (2019) give only qualitative indication and not definitive criteria of selection D parameter. In case of rockfill-like materials alternate to lava domains is uncertain what value to be used.
We use a value of D = 0 (Table 1) as primary choice because the most relevant problems are the possible weak and persistent discontinuous surface at the interface between rockfill-like materials and lava breccia with lava rock mass. Furthermore, in Di , it has already been demonstrated that the variation of D has little influence on the destabilization of volcaniclastic deposits on Stromboli.

Loading setting
The slope stability analysis is performed in static and dynamic (seismic action) conditions, to assess the reliability of the model. In static conditions, the slope is expected to remain stable, with an FS sufficiently larger than the unit. The capacity of the model to reproduce the shallow and small-scale instability of the slope increases the accuracy of the simulation. For the stability analysis under dynamic conditions, SSAP adopts the pseudo-static method and seismic coefficients. The horizontal force applied to the barycentre of each portion of sliding sling mass is equal to F = K h W where K h is the horizontal seismic coefficient and W is the portion of sliding mass's weight. K h is a function of the maximum horizontal acceleration and the topographical and lithological properties of the location (Baker et al. 2006;Chowdhury et al. 2009). A scenario is implemented in the model's validation, considering the slope being subjected to a seismic acceleration comparable to those generated by the 3 July 2019 paroxysmal explosion (Giudicepietro et al. 2020). Seismic data are collected from the STRG station of the INGV-OV (https:// www. ov. ingv. it/ ov/ it/ monit oragg io-sismo logico-di-strom boli. html), on the northern edge of the SdF, close to the shoreline (De Cesare et al. 2009). Seismic data are converted from "counts" to velocity data and then derived to obtain the acceleration. The horizontal seismic coefficient estimated is K h = 0.008.   Landslides 20 • (2023) The main purpose of the LEM analysis is to identify a sliding surface consistent with those proposed in the scientific literature and characterized by an FS's envelope correlated to the magnitude of the magma thrust. To reproduce the initial deep-seated α movement, the slope's response to a lateral force is investigated, namely the resultant of triangular pressure distribution (magmastatic), based on the magma's unit weight (γ m = 25.5 kN/m 3 ; Barberi et al. 1993). The susceptibility of the slope to a deep-seated landslide, comparable to the lateral collapse that created the SdF depression less than 5 ka (Gillot and Keller 1993), is evaluated by simulating a dyke intruding at the lower boundary of the volcaniclastic infilling material, reaching the topographic surface at 650-660 m a.s.l., which was the site of an effusive vent opening during the 2002-2003 flank eruption (Calvari et al. 2005).
After the preliminary global analysis, in static conditions and without external loading, we carried out an analysis to identify potential failure surfaces with a geometry comparable to the α slide (see Fig. 4). To simulate a tension crack, the identified failure surfaces have been selected to convert the curvilinear upper boundary into a vertical profile coherent with the original geometry. Then a lateral destabilizing force at the top of the single failure surface has been applied to calculate the FS (using the SSAP's tool "Individual Surface Analysis"). The implemented force is the resultant of triangular pressure distribution (magmastatic) and an excess pore pressure (overpressure), according to the equation: where F m is the magmastatic force, F e is the magma thrust due to overpressure, γ m is the magma unit weight, H is the height of the tension crack and P 0 is the overpressure. The reliability of modelling the slope response to dyke intrusion through the application of a lateral force on a vertical tension crack is attested by the numerous applications in scientific literature (Voight and Elsworth 1997;Apuani et al. 2005b;Tommasi et al. 2008;Battaglia et al. 2011).

Fig. 4 LEM model set-up.
In a is shown the entire flank, with the geotechnical model that best passed the test phase (static and dynamic) to which the forces in the vertical tension crack were applied, as reported in b. In b is reproduced the lateral thrust caused by the dyke. The elevation of the topographic profile is expressed in metres The main issue is that this approach is helpful in performing back-analyses of collapsing events, where the geometry of the sliding surface is already known but becomes less convenient for predictive analyses of potential failures. The exploration phase of potential sliding surfaces occurs before the application of the lateral thrust; thus, the surface's search engine is not influenced by the alteration of the stress domain caused by the dyke. That is because after the generation of potential sliding surface, the calculation of final FS associated to it is strongly influenced by the applied loads (static and/or dynamic). This methodology allows to consider in any case the full loading condition applied to ant type and shape of generated randomly possible sliding surface. The most critical surface thus will emerge due to assumed loading condition.
The dyke overpressure in vertical tension cracks can be applied only in a predefined surface. This option enables avoiding that the applied external resultant force is out of range of the physical sense of possible range of values, and at same time it gives the possibility to test a large set of possible values that can be applied at the same surface. This apparent limitation of the stability analysis with dyke overpressure is because the magmatic thrust can be modelled as a single force and not as a pressure distribution along a surface. In terms of lateral thrust magnitude, there is no difference between a pressure distribution and the corresponding resultant force. However, the function that assigns a lateral thrust to the head of the sliding force is designed for triangular pressure distributions, so that the resultant force is applied at 2/3H, where H is the total depth of the tension crack, measured from the topographic surface.

Slope condition and fracture field before the collapse
The 30 December 2002 helicopter survey began at Stromboli at 11:19:24 UTC, thus about 1 h before the first submarine landslide occurrence and lasted for 23 min and 53 s (Bonaccorso et al. 2003). When the survey started, the sea was "clear" and the SdF displayed one effusive vent located at 670 m a.s.l., erupting a small lava flow spreading north. At 11:21:22 UTC and during the following 3 min, the sea surface was still clear also along the coastline, indicating that no submarine landslide had occurred yet. The first image showing a little brown sediment within the seawater and up to a distance of a few metres from the coastline was recorded at 11:25:07 UTC, 3 min and 45 s after the start of the helicopter survey. The seawater was mostly clear during the initial phases of the helicopter survey and was observed for the first time "dirty" (i.e. muddy), presumably because of ash mixed within, from the coastline up to about 50-70 m from the coast at 11:40:48, after about 19.5 min from the start of the flight. This observation indicated that presumably some minor submarine landslide was occurring at that time, well before the submarine failure detached at 12:14:05 UTC, or the subaerial failures started at 12:22:38 UTC (Bonaccorso et al. 2003;La Rocca et al. 2004).
Starting at 11:34:25 UTC, 13 min after the start of the helicopter flight around Stromboli, the eastern fractures bounding the block that was going to fail began to emit ash. It was a sort of "passive smoking", with ash slowly rising from the slope all along the moving fracture, but this ash output was not associated with the release of gas under pressure typical of explosive activity. Conversely, the ash was emitted slowly and was rising almost passively from the slope. The process at the basis of this ash release was evident only after the slope collapsed, when we realized that the rising ash was caused by rock friction and fragmentation due to the sliding of the landslide block against the stable country rock. The ash release was displaying the opening of fractures along the slope, which we have represented in Figs. 5 (map) and 6 (graph). In general, during the helicopter flight we observed an expansion of fractures both upslope and downslope, with the length of the fractures increasing with time and merging until the limits of the subaerial landslides are identified (Figs. 5 and 6). At 11:41:04, for the first time the upper fractures doubled, displaying two parallel segments both releasing ash, and these expanded also downslope at 11:41:14 when another segment between 400-and 450-m elevation appeared on the right (west) side of the landslide scar.

Slope stability condition and scenarios
The static analysis recognized a set of sliding profiles analogous to the α slide (Fig. 7a, b), whereas in the submarine portion of the slope it identified several sliding surfaces comparable to the τ 1 + τ2 slides. Figure 7 shows that the area characterized by low local FS and high OSR values satisfies all geometric characteristics observed after the 30 December 2002 landslides : the head of the submarine landslides is located at the shoreline and the collapse scar has a maximum depth of 40-50 m ; the submarine landslide's toe can be placed at 300 m b.s.l. Similarly, an area affected by low local FS and high OSR values is located just upslope, in the subaerial part, resembling the β and the γ slides. Moreover, the OSR map (Fig. 7b) revealed that the area with the highest OSR is set at the α slide head, justifying the presence of a (magma-filled) tension crack in this area.
To establish the relationship between magmatic thrust and slope instability, the lateral thrust, initially equal to the magmastatic force, is progressively increased by 10 MN/m steps (Fig. 8a,  b). To assess the susceptibility of a flank collapse comparable to the destructive event that generated the SdF less than 5 ka (Lucchi et al. 2019), the dyke intrusion is modelled with the single-surface approach. A deep-seated curvilinear surface is collected from static global analysis. The surface with minimum FS was selected from the entire set of generated surfaces. The head of the slide, originally located in the crest area, is cut with a tension crack that starts at 513 m a.s.l. and rises to 675 m a.s.l., reaching vertically the topographic surface. The toe of the sliding surface is placed at 320 m b.s.l. (Fig. 3). Analyses are performed at increased steps of 10 MN/m, simulating different magma thrust. The outcomes of the stability analyses show that the slope destabilizes for significant magma thrust (3500 MN/m, Fig. 8a).

Discussion
The Sciara del Fuoco geotechnical stratigraphy A coherent lithological and geo-mechanical characterization of a volcanic slope is extremely difficult due to the considerable heterogeneities of the material composition and to the lack of extensive in situ and laboratory analysis. The analysis approach adopted for the Stromboli volcano favours the exploitation of a different distribution of geotechnical parameters within the slope rather than assuming different geotechnical parameters from those present in Landslides 20 • (2023) the scientific literature. The integration of geophysical data makes it possible to apply a dynamic loading that the slope underwent. Stability analyses performed in static and dynamic conditions highlight the role of the lava-breccia unit, a lithostratigraphic unit proposed for rock masses with analogous percentage of lava and breccia components. This unit can be considered as a lower boundary in the geotechnical characterization of the volcanic structure because the possible presence of weak interface between breccia facies (rockfill-like) and lava facies. In this case, the possible presence of long persistent interface with constant slope amplifies the possibility of a zone with low shear strength. Since generally geotechnical properties tend to improve with increasing depth, except for localized weak areas identified by specific geological investigations, it is reasonable to suppose that the internal part of the volcano is constituted by material with a shear strength equal to or greater than the lava-breccia unit. Another important aspect is the problematic conversion of parameters from the Hoek and Brown (HB) non-linear failure envelopes to Mohr-Coulomb (MC) linear failure criterion (Bromhead 1992). Given the extremely large range of confining pressures considered, the interpolation of HB parameters in MC criterion led to unrealistically high values of cohesion at shallow depth, as demonstrated by Li et al. (2008) and Shen et al. (2012). For this reason, generalized Hoek-Brown criterion (GHB) should always be preferred. The applicability of nonlinear HB failure criterion to non-rock mass geomaterials (Li and Yang 2019), as well as the fact that rockfill must be modelled in terms of non-linear failure criteria (Charles and Watts 1980;Indraratna et al. 1998;Liu 2009;Frossard et al. 2012;Barton 2008Barton , 2013, has been demonstrated. The geotechnical model of the slope that provides the more reliable results is constituted by three layers: the area unaffected by the last lateral collapses, defined by high values of friction and cohesion, according to MC criterion; the lava-breccia unit is assigned to most of the SdF infilling materials and a shallow subaerial layer is described with both, at the same time, GHB and the Barton-Kjaernsli strength criterion, in ubiquitous way.

Loading conditions
According to scientific literature (Acocella et al. 2006), the magmatic thrust exerted by a dyke's intrusion triggered the tsunamigenic sequence of landslides occurring on 30 December 2002. The dyke intrusion is modelled with a lateral force applied on the sliding surface tension crack. The main issue of this approach is that it is successful in performing back-analysis, since the geometry of the failure surface is already known, but the investigation of potential sliding surfaces is heavily affected by the choice of the failure contours. This limitation does not affect this study, since the geometry of the 2002 α movement and of the lateral collapse occurring in the past are estimated by previous authors. The SSAP's simulation in static conditions agree with the topo-bathymetric measurements made after the collapse, i.e. the area with the lowest FS and highest OSR coincides with the area that actually collapsed, even if with FS values slightly above 1, confirming that a perturbation to the system (magma thrust, grain crushing, etc.) is capable of inducing the slope itself to collapse in conditions of metastability.
Stability analysis performed on the deep-seated sliding surface comparable to the lateral collapse occurring 5 ka (Gillot and Keller 1993;Lucchi et al. 2019) shows that a large value of magmatic thrust is necessary to destabilize the slope. Some comparisons can be made with previous studies on the deep-seated instability of the Stromboli NW flank. Apuani et al. (2005b) designed a translational deep-seated sliding surface bounded at the top by a 220-m tension crack, with its toe placed at the shoreline and crossing a lava-breccia layer where a disturbance factor (Hoek et al. 2002) D = 0.6 was assigned. They observed that, given that stratigraphy, instability can be reached adding an overpressure equal to 0.5 MPa to the magmastatic pressure. Tommasi et al. (2008) performed a stability analysis, in the 3D domain, to determine the critical geotechnical parameters of the volcaniclastic material within the SdF for which the failure occurs. The magma thrust was modelled with only the magmastatic component. Both lithostratigraphic characterizations proposed by these authors have been tested in the previous paper but neither showed a reliable response in static conditions. Tommasi et al. (2007) andVerrucci et al. (2019) evaluated the consequences of a dyke intrusion on Stromboli with FEM analysis. They observed that even for large values of magmatic thrust, the slope underwent major deformation near the dyke location, but the stability of the entire flank was always ensured. Therefore, they proposed that the magma must also intrude below the SdF volcaniclastic deposits to generate large instability phenomena. These findings are also in agreement with the evidence of volcano tectonics and structural geology studies, which show a propensity for small-sized dykes in the current conditions (e.g. 2002-2003, 2007, 2014 eruptions), which instead change  Fig. 7 The figure shows the two 2D colour maps, defined as quasifinite element analysis (qFEM) and over stress ratio (OSR) maps, considering static conditions along the entire Sciara del Fuoco slope, and thus representing the back-analysis of the 30 December 2002 landslides sequence (pre-intrusion). a qFEM map and b OSR map into larger dykes in the presence of greater in-filling thicknesses in the Sciara del Fuoco ).
Since the magnitude of the lateral force is related to height of the tension crack (H), it is important to carefully choose how to design the tension crack. The correlation between H and the variation of FS for different values of applied pressure makes this approach to the study of the dyke intrusion appropriate for back-analysis where the geometry of the surface is well-known. Instead, the estimation of potential failure in volcanic environments must be carefully carried out, assuming several tension crack geometries based on the volume of material involved and on the dyke location.

The 30 December 2002 landslides: a collapse in four acts
The sequence of events observed before and during the 30 December 2002 collapse, consisting of an increase in eruptive parameters as a consequence of the magma stationing at high altitudes in the volcano (i.e. below the crater terrace), opening of the eruptive fracture that started from the NEC and consequent landslide of the external portion of the NEC, was later observed also in the early stages of the 2007 (Casagli et al. 2009;Marsella et al. 2009)  -Run-up phase (Fig. 9a). A preparatory phase for the eruption, not defined as true unrest because in the case of a persistently active volcano this can be misleading. During this phase, there is an increase in the magma that accumulates under the crater terrace with related phenomena (more intense activity, overflows, small landslides from the edge of the crater terrace, frequent shifting of the activity from the central to the lateral vents, most commonly the NE crater). These phases can evolve into flank eruptions, as happened in   (Fig. 9b). As the magma input rate increases, the magma residing below the crater terrace of Stromboli is no longer able to depressurize itself through eruptive activity (explosions + overflows). Furthermore, the existence of a structural limit to the inflation of the crater terrace has been hypothesized (Di Traglia et al. 2014), and therefore the magma filling the dyke propagates laterally. This propagation triggered the collapses of the NEC talus (or southern talus), produced avalanches of hot rocks and fed the first lava flows with vents located at high elevations. This occurred in (Calvari et al. 2005Burton et al. 2008;Casagli et al. 2009;Di Traglia et al. 2018a;2022). In the run-up and onset phases, magma densification caused by frequent explosive activity may have played a role in the increase in magmastatic pressure and therefore in the instability of the summit areas, as demonstrated by Di . -Intrusion in the Sciara del Fuoco slope (Fig. 9c). This phase creates many deformations and fractures in the upper part of the SdF, feeding lower-elevation vents (500-400 m a.s.l.). Sheet intrusion in the SdF only occurred in (Acocella et al. 2006Casagli et al. 2009;Neri and Lanzafame 2009;Calvari et al. 2010). -Deep-seated deformation and landslides collapse (Fig. 9d).
The intrusion phase in the SdF can stop (for example owing to lithological-structural discontinuity, (Marsella et al. 2009;Intrieri et al. 2013), or evolve with the deformation of the entire slope, with the propagation of fractures towards low altitudes (as highlighted in this work). The deformation of the entire SdF testifies to the triggering of a deep-seated slope deformation (what in 2002 was called α slide), which in turn triggers the detachment of more superficial slides. This involvement of the whole SdF slope can take place without the dyke intruding at low altitudes if the magma thrust is sufficiently high (as evidenced by the results of our stability models).

Conclusions
This critical review of the data published on the sequence of tsunamigenic landslides of Stromboli occurred on December 30, 2002, has allowed us to better understand the mechanisms leading to the collapse of a part of the Sciara del Fuoco, the unstable slope of Stromboli volcano. The images collected from a helicopter survey less than 1 h before the collapse displayed the propagation of the fractures along Stromboli's NW slope just before the failure and improved our interpretation of the causes triggering the instability. The definition of a reliable geotechnical stratigraphy of the unstable flank was based on pre-existing data rearranged and tested using numerical models with the limit equilibrium method including different types of failure criteria. This method to combine different failure criteria for the geomechanical units enabled understanding both the predisposing factors and the triggering mechanism of the landslide sequence. From the analysis of the eruptive activity that anticipated the landslides, it is possible to identify 4 evolutionary phases: (1) a runup phase of the eruptive activity, with the intensification of the Strombolian explosions and the occurrence of lava overflows from the crater terrace; (2) the onset, with the progressive migration of magma from the central areas of the crater terrace, towards the lateral vents (NE or southern vents), triggering a small landslide that was used to evolve as gravel flows (e.g. 2014) or a proper PDC (e.g. 2002); (3) intrusion of the SdF, with the continuation of lateral and the simultaneous downward migration of the magma, with associated slope deformations and the opening of some effusive vents along the volcano's flank; (4) slope deformation and collapse, without the magma propagating at low altitudes. This destabilization, on 30 December 2002, induced the flank movement and the subsequent collapse of more shallow parts (both submarine and subaerial), which were responsible for triggering the tsunami waves. The stability analysis also revealed that for collapses of larger dimensions to occur, a very high magma thrust is necessary, suggesting a complex mechanism for these events.
The approach used in this work, which integrated the geotechnical modelling and analysis of the slope through a stability analysis with LEM and the observations acquired just before the collapse, can be generalized and applied to other instability contexts of volcanic edifices. In particular, the association of volcanoclastic units and lava-breccia units with an adequate failure criterion, such as the combination of generalized Hoek-Brown criterion with Barton-Kjaernsli criterion, used to describe the levels constituted by the stratified alternation of volcaniclastic deposits and lava, such as stratovolcanoes, can also be generalized to other volcanoes characterized by persistent or semipersistent activity. This combination of well-known non-linear criteria in determining the appropriate local shear strength values depending on the stress condition and local slope of the assumed sliding surface is promising for the application to the complex nature of stratovolcano flank deposits. The proposed new approach may overcome the limitation of parametrization and shear strength modelling in a continuous domain only (M-C or HB failure criteria).

Funding
Open access funding provided by Istituto Nazionale di Geofisica e Vulcanologia within the CRUI-CARE Agreement. This research was funded by the "Presidenza del Consiglio dei Ministri-Dipartimento della Protezione Civile", through the UniFi-DPC 2019-2021 agreement (Scientific Responsibility: N.C.). The contents of this paper represent the authors' ideas and do not necessarily correspond to the official opinion and policies of the "Presidenza del Consiglio dei Ministri-Dipartimento della Protezione Civile". This research was also funded by the Project FIRST-ForecastIng eRuptive activity at Stromboli volcano: Timing, eruptive style, size, intensity, and duration, INGV-Progetto Strategico Dipartimento Vulcani 2019 (Delibera n. 144/2020; Scientific Responsibility: S.C.). The SSAP software research and development was funded by CONACYT (Mexico): Proyectos Ciencia Basica: CB-2016/286764.

Conflict of interest The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.