Crater-rim collapse at Stromboli volcano: understanding the mechanisms leading from the failure of hot rocks to the development of glowing avalanches

The failures of volcanic crater-rims frequently lead to the development of avalanches of glowing rocks, a hybrid phenomenon between rock avalanches and pyroclastic density currents, reaching considerable distances from the eruptive centres and therefore be a serious threat for inhabited areas. The collapse conditions have been here taken in to account by means of slope stability analysis using a Limit Equilibrium Method analysis on the crater-rim of Stromboli volcano (Italy). A Stromboli, crater-rim collapses occurred frequently (at least seven events in the last two decades) and were always associated with high-level of magma within the conduits, testified by the increased eruptive activity and ground-deformation. The more frequent/intense eruptive activity produced a greater accumulation of volcaniclastic material, whereas the high level of magma increased magmastatic thrust on the deposits. Volcaniclastic material was modelled by combining the nonlinear failure envelopes as the Generalized Hoek and Brown criterion, with the addition of the failure's envelope of the rockfill-like material described by the Barton-Kjaernsli criterion, taking into consideration the presence of discontinuities within proximal, partially welded, volcaniclastic masses. In addition to the lithological and morphological characteristics of the crater terrace rim and the magmatic thrust, the effects of the explosions in terms of seismic ground acceleration and disturbance factor (D) of the volcaniclastic material were also considered here. While the ground acceleration compatible with the explosive activity of Stromboli has little influence on the stability of the crater terrace rims, the increase in D increases the proneness for failure.


Introduction
The collapse of the unstable portion of a volcanic crater's rim usually generates glowing avalanches, a mixture of hot volcanic rocks, mainly bombs and spatter agglutinates, and ash-and-lapilli sized volcaniclastic matrix that has intermediate characteristics between pyroclastic density currents and rock avalanches (Davies et al. 1978;Nairn and Self 1978;Cole et al. 2005). Avalanches of incandescent material have been reported for the Fuego volcano Guatemala (Davies et al. 1978;Risica et al., 2022), as well as for the 1975 eruption at Ngauruhoe volcano, New Zealand (Nairn and Self 1978;Lube et al. 2007). Deposits associated with these phenomena are present on the flanks of Arenal, Costa Rica (Cole et al. 2005), Shin-Fuji volcano, Japan (Yamamoto et al. 1993), Asama volcano, Japan (Yasui and Koyaguchi 2004), and Tungurahua volcano, Ecuador (Kelfoun et al. 2009), and have been observed in their genesis on three Italian volcanoes: Etna (Calvari and Pinkerton 2002;Behncke et al. 2008;Andronico et al. 2018), Vesuvius (Hazlett et al. 1991;Arrighi et al. 2001), Stromboli (Pioli et al. 2008;Di Traglia et al. 2014Calvari et al. 2016Calvari et al. , 2020. These are events that involve small volumes (10 4 -10 7 m 3 ) and can travel a few kilometres from the source area but are emplaced at very high temperatures (Davies et al. 1978;Nairn and Self 1978;Hazlett et al. 1991;Cole et al. 2005;Calvari et al. 2016Calvari et al. , 2020, being potentially dangerous for communities close to the volcanoes and tourists. To improve hazard assessments associated with these phenomena, it is necessary to expand the understanding hot rocks failure, considering the different lithological factors, topographic conditions, triggering mechanism. This paper presents the results of some stability analyses obtained by modelling with the limit equilibrium method (LEM), taking into consideration the Stromboli volcano, for which there are many supporting data. For the reconstruction of the conceptual model that leads to the instability of the Stromboli crater-rim, both the mechanical data of the Stromboli rock masses (i.e. Apuani et al. 2005a, b;Tommasi et al. 2007;Boldini et al. 2009;Rotonda et al. 2010) and the geophysical and geomorphological monitoring data (i.e. Calvari et al. 2016;Di Traglia et al., 2014 present in the literature were critically reviewed, to be able to correctly analyse both the predisposing phenomena and the triggering mechanisms of the collapses. The results of this work can be used either to correctly understand the monitoring data (i.e. ground deformation; topographic variation) and can be used as input for propagation models of hot material flows, and therefore for the reconstruction of the impact scenarios connected to these phenomena at Stromboli. More generally, the approach proposed here can be applied to study the same phenomena also in other volcanoes, adapting the stability models to the correct lithological, morphological and trigger conditions (Table 1).

Crater-rim collapses at Stromboli volcano
Stromboli is an island volcano in the southern Tyrrhenian Sea, being part of the Aeolian archipelago. The volcano, morphologically dominated by unstable flanks (NW and SE), is characterized by low to medium intensity explosive activity called Strombolian (Barberi et al. 1993;Rosi et al. 2013). These consists in the ejection of coarse pyroclastic material, which generally falls in its crater terrace (Fig. 1). Stronger explosions usually lead to the ejection of coarse clasts (bombs and blocks) that, following ballistic paths, can impact

Crater-wall weakening
Reduction in the strength parameters Schaefer et al. (2015) Erosion Thinning of crater wall Calvari et al. (2016) Overloading

Tephra accumulation
The accumulation of material (lava or tephra) on the crater-rim or volcano flank produces overloading on the volcano slope Nolesini et al. (2013) Magma thrust

Magmastatic variation
Magma thrust due to the high level in the conduit Di Traglia et al. (2014)

Dike intrusion
Thrust of magma in fluid-filled cracks Voight and Elseworth, (1997) Viscous indentation Magma thrust due the formation of cryptodome structures or massive bulbous intrusions, and outward bulging and flaring when intruding poorly consolidated rocks Cole et al. (2005) 1 3 zones outside the crater terrace. During these events, ash plumes are usually generated from a few hundred meters (major explosions) up to a few kilometres in height (Strombolian paroxysms), leading tephra dispersion on the island (e.g. Andronico and Pistolesi 2010) or, sporadically, on larger areas (e.g. Giordano and De Astis 2021). A recently proposed classification by  sought to combine classifications based on the effects of explosions with classifications based on geophysical signals measured during the events. In particular, the proposed classification is based on the VLP size and the VD parameter, being the former the maximum value of the RSAM of a 30-s sliding window that moves in a 30-min time interval of signal, filtered in VLP frequency band (0.05-0.5 Hz), whereas the latter is the product of the muzzle velocity and the explosion duration, both derived from the analysis of the monitoring camera images. Using these two devices, used to measure the crater terrace displacement, has been also reported parameters, it is possible to identify 4 classes (from 0 to 4 in ascending order), where class 0 includes the ordinary Strombolian activity, class 1 the "intermediate" explosive activity (according to the classification of Andronico et al., 2008), class 2 the "major explosions", class 3 includes the "Strombolian paroxysms" (also called Vulcanian basaltic, see Calvari et al. 2006). The Stromboli explosive activity can be accompanied or replaced by effusive activity, as lava overflows from the crater area (e.g. Di Traglia et al. 2014;Calvari et al. 2016), or as lava flows fed by ephemeral vents (e.g. Di Traglia et al. 2018) within the unstable NW flank (locally called Sciara del Fuoco).
Leaving aside the PDCs column-collapse events, such as those observed during class 3 explosive events (and only sporadically during class 2 explosions), avalanches of incandescent material are generated in Stromboli due to two main mechanisms: • Crater-rim collapse PDCs (sensu Cole et al. 2005), being related to the accumulation of material on the rim of the crater terrace, and are triggered by dike injection (Calvari et al. 2005;Casagli et al. 2009;Di Traglia et al. 2018), magma fingering (Cole et al. 2005;Calvari et al. 2016). • Deposit-derived gravitational collapse PDCs (sensu Sulpizio et al., 2014 andRisica et al., 2022), being generated by the failure of material deposited during the paroxysmal explosions on the volcano flanks (Di Roberto et al. 2014;Salvatici et al. 2016). These deposits, generally spatters that were deposited on loose volcaniclastic material, collapses by simply exceeding the friction angle (Di Roberto et al. 2014;Salvatici et al. 2016). This occurred at least two times at Stromboli in the last century, during the 1930 and 1944 paroxysms (Rittmann 1931;Di Roberto et al. 2014;Salvatici et al. 2016), when small-volume flows of hot pyroclastic material spread on the volcano flanks, reaching the sea and, in the case of the 1930 event, causing casualties and destruction along its path (Rittmann 1931).
This work takes into account the generation of crater-rim collapse PDC for two main reasons: i) glowing avalanches are direct risk factors in the areas where they can flow (the stretch of sea in front of the Sciara del Fuoco, in the case of Stromboli); ii) often associated with overflows and the early stages of flank eruptions, they can be indicators of conditions leading to the transition from explosive to effusive activity. This explosive to effusive transition, in the case of Stromboli, is the condition that can lead to the dikes intrusion with consequent destabilization of the unstable flank, triggering of landslides, and tsunamis generation (Acocella et al. 2006;Di Traglia et al. 2014).
The observations made in the last 30 years of Stromboli's activity, supported by an advanced volcanological, geophysical, geochemical, and geomorphological monitoring system, have allowed us to identify two main phenomena (Table 2): • the collapse, triggered by overflows, of the hornitos that form on the crater-rim during the spattering activity (Calvari et al. 2014(Calvari et al. , 2016(Calvari et al. , 2020Di Traglia et al. 2014. These phases are frequently observed in Stromboli (Coppola et al. 2012), and only in some cases were they anticipatory of effusive flank eruptions (Di Traglia et al. 2018).
During these periods of more frequent/intense activity, the geophysical parameters usually increase (Calvari et al. 2014(Calvari et al. , 2020, especially with deformations of the crater terrace and increased seismic tremor (Di Traglia et al. 2014;Calvari et al. 2020). These phases are also characterized by an increase in surface degassing, highlighted by the increase in SO2 emissions (Calvari et al. 2014); • the collapse, triggered by magma-fed fractures at the onset of flank eruptions, of the volcaniclastic talus that form around the crater terrace due to the accumulation of the material produced by Strombolian activity. This phenomenon was observed in the early stages of the 2007 and 2014 flank eruptions (Casagli et al. 2009;Di Traglia et al. 2018), but a deposit associated with this phenomenon has been described by Pioli et al. (2008) for the 2002-2003 eruption. Glowing avalanches have also been described in the early phase of the 1985 flank eruption (De Fino et al. 1988).
The 2014 effusive eruption revealed that the precursors of flank eruptions can be very similar to those of overflow activity (i.e. more intense / frequent explosive activity, deformation of the summit area, increased tremor) and are associated with ascents of significant magma mass below the terrace crater (Di Traglia et al. 2018). The 2014 eruption was anticipated by 2 months of "anomalous" activity, during which many overflows occurred, in turn often associated with crater-rim collapses (Di Traglia et al. 2018;Casalbore et al. 2021). The failure of a hornito that had formed on the edge of the crater terrace occurred on 6 August 2014 (Di Traglia et al. 2018). This collapse had similar features to others observed at Stromboli, such as the one that occurred on 12 January 2013 (Calvari et al. 2016), or 19 May 2021 (Fig. 2). However, while for these latter events the eruptive activity decreased as a result of the overflows associated with the collapses, after 6 August 2014 failure the activity continued to increase (Di Traglia et al. 2018). This led to the lateral propagation of the magma stored beneath the crater terrace, with the opening of an eruptive fracture early on 7 August 2014, which triggered the collapse of the volcaniclastic debris talus located immediately below the NE area of the crater terrace (Di Tioukov et al. 2022), allowing the magma to feed the 2014 flank eruption (Di Traglia et al. 2018) (Fig. 3).  1 3

Stability analysis: methods and settings
The stability of the volcaniclastic deposits has be analysed using a 2D LEM by means of the Slope Stability Analysis Program (SSAP; Borselli et al. 2011Borselli et al. , 2020, a free software that examines the Factor of Safety (FS) of sliding surfaces through 2D LEM analysis by assigning the geo-mechanical parameters, according to different failure criteria (Mohr-Coulomb linear failure criterion, generalized Hoek and Brown and/or Barton-Bandis nonlinear failure criteria), already used to estimate the stability of volcanic edifices in different geological contexts (Borselli et al. 2011;Dondin et al. 2017). The FS is here calculated using the Morgenstern-Price method (Morgenstern and Price 1965), whereas the research of critical surfaces is carried out through the "Random Search" engine, using the Monte Carlo methods. SSAP implements the "Dynamic Surface Search Attractor", a feature that progressively reduces the initial search area according to the surfaces with lower FS gradually identified, allowing for the identification of more critical surfaces within a specific area, also discarding potentially unstable areas if not identified in the early stage of the search process. SSAP outputs consist in: (i) a graphical representation of the surfaces with lowest FS; (ii) the internal distribution of forces and pressure, distribution of local FS, local distribution of the numerical reliability of general FS numerical solution (namely RHO index), (iii) a 2D colour map with distribution of average local FS obtained by local stress distribution, based on a combination of LEM and Finite Element Analysis (FEM) procedures, defined as quasi-FEM (qFEM; Schofield and Wroth 1968;Griffiths and Lane 1999), describing the local distribution of the FS, including a representation of the potential direction of plasticization for those areas where FS < 1. Although the actual state of deformation is not calculated by the software (SSAP does not allow the user to include the strain's parameters in the slope characterization), these features provide a graphical representation of those areas in which progressive failure can originate. A constitutive model of the slopes is created by using the mechanical parameters present in the bibliography (Apuani et al. 2005a,b;Tommasi et al. 2005Tommasi et al. , 2007Tommasi et al. , 2008Boldini et al. 2009;Casagli et al. 2009;Rotonda et al. 2010;Nolesini et al. 2013;Verrucci et al. 2019a, b). 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 geomechanical/geotechnical properties. Therefore, the physical and mechanical characterization data present in the literature has been critically analysed. In this work, the constitutive model of the slope of the Stromboli crater-rim has been reproduced considering: • nonlinear failure envelopes as 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), accurately and efficiently describe the behaviour of the volcanic edifice both at shallow and higher depths; • failure's envelope of the rockfill material described by the Barton and Kjaernsli (1981) criterion, already implemented in SSAP as described by Barton and Bandis (1991), Barton (2013) and Prassetyo et al. (2017). This method has already been proposed for the filling materials of the Sciara del Fuoco, the unstable slope of Stromboli by Apuani et al. (2005a, b), Tommasi et al. (2007), Boldini et al. (2009), and Rotonda et al. (2010).
Regarding the loading settings, the slope stability analysis has been performed in static and dynamic conditions, to assess the reliability of the model. In static conditions, the slope is expected to remain stable (FS > 1) or metastable (FS ~ 1). 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 wedge is equal to F = K h W where K h is the horizontal seismic coefficient and W is the wedge's weight. K h is a function of the maximum horizontal acceleration and the topographical and lithological properties of the location. A scenario is implemented in the model's validation, considering the slope being subjected to a seismic acceleration from K h = 0 to K h = 0.008, comparable with the ground motion recorded during the 3 July 2019 paroxysmal explosion (Giudicepietro et al. 2020).
Once the constitutive model of the slopes has been validated, the stability of the incandescent rock masses is estimated considering the thrust of the magma. To simulate the conduit suppling the vents on the crater-rim and, as well as a magma-filled tension crack propagating with the volcaniclastic debris talus, the failure surfaces identified by the LEM analysis is 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 will be applied to calculate the FS. 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 and F e is the magma thrust due to overpressure. The reliability of modelling the slope response to dike 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).
The slope's topography is obtained by using available Digital Elevation Models in Di Traglia et al. (2020), Di ). 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 Table 3 Geomechanical parameters used in the LEM simulations. The geomechanical parameters are: Uniaxial Compressive Strength of the intact rock (UCS); Geological Strength Index of the rock mass (GSI); material constant mi (textural and petrographical character of the intact rock); factor of disturbance D; dry bulk volume (γ dry ); saturated bulk volume (γ sat ); joint roughness coefficient (JRC); joint compressive strength (JCS); joint residual friction angle (φ r ); joint length (L); laboratory-scale joint sample length (L 0 ); slope face angle (β) and its variation (Δβ) Geomechanical parameters following the GHB (Hoek et al. 2002) Layer  Barton and Kjaernsli (1981) Criterion following parametrization of Barton and Bandis (1991) and Lunardi et al. (1994) Layer in dry condition). This assumption is also supported by meteorological and geophysical observations reported in Revil et al. (2011). The SdF behaviour in response to both static and loading conditions has been reproduced considering a mechanical stratigraphy as follows (Table 3): • a shallow layer in the subaerial part of the slope made up of volcaniclastic deposits modelled as rockfill material (Tommasi et al. 2005); • Lava-Breccia unit, made by a mixture of breccia and lava layers (described by Apuani et al. 2005a), modelled with the nonlinear failure envelopes as the Generalized Hoek and Brown criterion GHB (Hoek et al. 2002).

Stability analysis: results
The degree of stability of the crater terrace rim has been quantified by means of a 2D stability analysis using LEM on the profiles 1, whereas the stability of the volcaniclastic talus has been quantified by LEM analysis on the profiles 2 (Fig. 4). From the stability analysis carried out on the two profiles it can be deduced that the volcaniclastic accumulations have conditions of diffuse instability, and that in general they are in metastable conditions (FS). The local FS maps (Fig. 5) show how the variation of the factor D, which depends upon the degree of disturbance to which the rock mass has been subjected by explosion damage (Hoek and Brown 2019), do not significantly influence the stability of the slope. In the graph of Fig. 6a it is possible to estimate how the seismic acceleration induced by the explosions has little influence on the variation of FS in all the volcaniclastic accumulations taken into consideration. The estimate of the effect of the Kh on the stability conditions was carried out by varying the factor D, which increasing reduces the value of FS, without however reaching critical FS values.
From the graphs in Fig. 6b, c, it is possible to observe the variation of FS as a function of the variation of the magmatic thrust in the single surfaces analysed on Profile Fig. 4 a localization of the profiles used for the stability analysis. In the background map: difference between the elevations calculated using DEM reported in Di . b profile 1 (crater terrace rim) and c profile 2 (volcaniclastic talus). The areas inside the crater terrace are affected by uncertainty about the DEM, because the DEMs derived from Plèiades-1 satellite data are noisy due to the plume of volcanic gas emitted by the volcano (see Di Traglia et al. 2022 for DEM description) 1 and 2, respectively. In general, the collapse of the crater terrace rim is reached with lower magmatic thrust values, slightly higher than those induced by the magmastatic thrust. For the latter, the dense rock equivalent (DRE) value of 2,750 kg/m 3 for Stromboli magma (Barberi et al. 1993) and value of 2,750 kg/m 3 of "high-density" scoriae reported in Houghton (2005, 2007) were used. These values correspond to 0.9 MN/m and 1.6 MN/m for the crater terrace rim (Fig. 6b), whereas the magmastatic thrust is in the range of 3 MN/m and 5.4 MN/m for the volcaniclastic talus (Fig. 6c). Stability is reduced by the increase of D.

Discussion and conclusive remarks
The collapses of the crater-rim, regardless of their location and the thickness of the volcaniclastic accumulations that are involved in the collapse, are the final phase of a process of increased eruptive activity, in turn following the ascent to more superficial levels of the  (Calvari et al. 2016;Di Traglia et al. 2018). The collapse of the rim of the crater terrace is associated with the phenomena of spatter accumulation and lava overflow (at least 1 overflow per year, see Coppola et al. 2012;Calvari et al. 2014;Di Traglia et al. 2014;Di Traglia et al. 2018;Di Traglia et al. 2020;Di Traglia et al. 2022), which are much more frequent than flank eruptions (only 3 flank eruptions since 2002). These processes have a fourfold effect: 1. more material is accumulated due to more frequent explosions. This changes the topography (especially the slopes) and the material available. The greater accumulation of material was simulated considering volcaniclastic accumulations comparable to the "pre-collapse" ones ( Fig. 5). Furthermore, the choice of using rockfill-like material, i.e. considering the discontinuities as actually observed in Stromboli (Tommasi et al. 2007), well reproduces unstable areas are those of material accumulated during the eruptions, and not the older material (evidently more stiff, as proposed by Apuani et al. 2005a); 2. the erupted material is disturbed by the explosive activity itself (increasing D). The disturbance effects of the explosive activity were considered by varying the parameters D and the seismic load Kh (Fig. 6), which are factors that add up to the others, but which alone are not significant for triggering collapses. This is consistent with the observation that crater-rim collapses have not been reported during ordinary Strombolian explosive events; a value D > 0 was considered also in previous studies on active volcanoes (Apuani 2005a, b, Borselli 2011). 3. the rise of the magma level in the conduits increases the magmastatic thrust. This effect was simulated (Figs. 6b, c), consistent with the observations made in Stromboli, as well as in other volcanoes (e.g. Arenal, Cole et al. 2005); 4. the magma in the conduits increases in density due to the loss of gases emitted during high passive degassing and frequent explosive activity. The effect of "densification" of magma due to gas loss was considered (Fig. 6b, c), which is consistent with the magma's inability to continue explosive activity and migrate laterally, triggering overflows or, if the magma supply from the depths is greater, leading to the development of eruptive fractures, and therefore to the onset of flank eruptions.
What appears from the analyses carried out is that the rim of the crater terrace is in more critical conditions than the volcaniclastic talus (Fig. 5), especially if we consider surfaces such as those analysed in profile 1 (i.e. the collapse of 19 May 2021) where magmatic thrust values slightly higher than the magmastic one are required, compatible with magmatic overpressure much lower than 1 MPa. On the contrary, magmatic overpressure values of ≈ 1 MPa are necessary to induce the collapse of the volcaniclastic talus, in accordance with what is calculated for the magma lateral migration from the crater terrace to the ephemeral vents during the initial phases of the flank eruption (Di Traglia et al. 2014). Fig. 6 a Relationship between Kh and FS variation along the selected surfaces for both the crater terrace rim and the volcaniclastic talus, which suggests a minimal effect of the ground motion induced by the explosions on the stability of volcaniclastic accumulations. b Relationship between magma thrust and FS along the a surface located on the crater terrace rim (h = 11 m). Magma densification due frequent explosions account for change in magma thrust from 0.9 to 1.6 MN/m. c Relationship between magma thrust and FS along the a surface located in the volcaniclastic talus (h = 20 m). Magma densification due frequent explosions account for change in magma thrust from 3 to 5.4 MN/m 1 3 The results of the stability analysis of the crater-rim of Stromboli can be considered a very important contribution also for the study of instability phenomena in other volcanoes, characterized by similar phenomena. An example comes from the Arenal Volcano in Costa Rica. Cole et al. (2005) had hypothesized that the initiation of crater-rim collapse occurred during the waning phase of the explosive activity, with the indentation of the crater-rim by degassed magma which fed lava flows and sporadically triggered crater-rims collapse PDCs.
Gravitational deformations of volcaniclastic deposits in summit areas triggered by dike intrusions have been measured by Bonforte and Guglielmino (2015) and Schaefer et al (2015Schaefer et al ( , 2019 on Etna and Pacaya, respectively. These cases should be analysed with an approach similar to the one proposed here. At Stromboli, the volumes involved in the crater-rim collapse are in the order of magnitude of 10 5 -10 6 m 3 , and largely depend on the availability of new material accumulated on the rim of the crater terrace. These kind of analyses can be generalized to other contexts, such as the collapses of material deposited on the flanks of volcanoes (e.g. Salvatici et al., 2016;Risica et al., 2022). In this case, stability analyses based on multi-temporal topographic measures (e.g. Casalbore et al., 2022) can be fundamental to estimate the actual volume of the material prone to collapse, and therefore to make adequate and repeated hazard assessments over time.