# Explosive eruptions at Stromboli volcano (Italy): a comprehensive geochemical view on magma sources and intensity range

Marija Voloschina, Nicole Métrich, Antonella Bertagnini, Paola Marianelli, Alessandro Aiuppa, Maurizio Ripepe, Marco Pistolesi

► To cite this version:

Marija Voloschina, Nicole Métrich, Antonella Bertagnini, Paola Marianelli, Alessandro Aiuppa, et al. Explosive eruptions at Stromboli volcano (Italy): a comprehensive geochemical view on magma sources and intensity range. Bulletin of Volcanology, 2023, 85, 10.1007/s00445-023-01647-y. insu-04155709

HAL Id: insu-04155709

https://insu.hal.science/insu-04155709

Submitted on 8 Jul 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Explosive eruptions at Stromboli volcano (Italy): a comprehensive geochemical view on magma sources and intensity range

Marija Voloschina<fnref n="1" /> Ⓓ · Nicole Métrich<fnref n="2" /> · Antonella Bertagnini<fnref n="3" /> · Paola Marianelli<fnref n="1" /> · Alessandro Aiuppa<fnref n="4" /> · Maurizio Ripepe<fnref n="5" /> · Marco Pistolesi<fnref n="1" />

Received: 28 September 2022 / Accepted: 3 April 2023 / Published online: 5 May 2023
© The Author(s) 2023

## Abstract

A comprehensive understanding of the processes that occur during magmatic storage and pre-eruptive ascent—and of their associated timescales—is critical to identifying potential precursory signals, and to developing robust volcano early-warning systems. Stromboli's persistent activity comprises continuous degassing and explosive activity that ranges from hourly, low-intensity “normal” activity to occasional, more violent, paroxysmal activity. While the magma source processes that drive normal and paroxysmal activity are reasonably constrained, eruptive activity intermediate in magnitude and intensity (i.e., major explosions) remains elusive in terms of classification, source region, and pre-eruptive timescales. Here, we investigate the 19 July 2020 major explosion that geophysical parameters place at the upper limit of the major explosions field, close to small-scale paroxysms such as the 2003 and 2007 events. The geochemical signatures of matrix glass, olivine, melt inclusions, and embayments—integrated with gas measurements—highlight important differences in eruption source, ascent behaviour, and pre-eruptive timescales of the studied event when compared to paroxysms. Melt inclusion volatile contents identify that magma rise begins from a slightly shallower source (~9.5 km below sea level, b.s.l.) than for paroxysms (11.4 km b.s.l.), with the activation of a shallower ponding zone at 5–6 km b.s.l.. This, in combination with intermediate matrix glass compositions, suggests complex ascent behaviour, characterised by CO₂ buffering in the deep ponding region and magma self-mixing in the shallower zone. Fe–Mg-diffusion modelling in olivine indicates a system perturbation starting ~20–25 days before eruption onset, in agreement with the timescale of volcanic gas CO₂/SO₂ ratio changes observed in the plume, and significantly shorter than that observed prior to paroxysms (~4 months). The geochemical dataset provides insights into the processes controlling the steady-state conditions and the broad spectra in eruption magnitude and intensity at Stromboli and bears important implications for eruption forecasting.

## Riassunto

Un approccio integrato alla comprensione dei processi che avvengono durante lo stazionamento del magma e la sua risalita, incluse le tempistiche associate, è fondamentale per identificare eventuali segnali precursori e per lo sviluppo di sistemi di allertamento affidabili. L'attività persistente del vulcano Stromboli è caratterizzata da un degassamento continuo e da attività esplosiva che varia dall'attività normale, con frequenza oraria e di bassa intensità, ad attività parossistica, meno frequente ma più violenta. I processi sorgente che alimentano l'attività normale e quella parossistica sono ben vincolati, mentre classificazione, sorgente magmatica e tempistiche pre-eruttive sono meno costrette per l'attività esplosiva intermedia in magnitudo e intensità (esplosioni maggiori). Il presente lavoro indaga l'esplosione maggiore del 19 luglio 2020 che, sulla base dei parametri geofisici, ricade al limite superiore del campo delle esplosioni maggiori, prossima ai parossismi del 2003 e del

Editorial responsibility: J.E. Hammer

<figure><img src="image_2.png" /></figure>

<footnote>Marija Voloschina
marija.voloschina@dst.unipi.it</footnote>

<footnote n="1">¹ Dipartimento di Scienze della Terra, Università di Pisa, Pisa, Italy</footnote>

<footnote n="2">² Université Paris-Cité, Institut de Physique du Globe de Paris, CNRS, Paris, France</footnote>

<footnote n="3">³ Istituto Nazionale di Geofisica e Vulcanologia, Sezione di Pisa, Pisa, Italy</footnote>

<footnote n="4">⁴ Dipartimento di Scienze della Terra e del Mare, Università di Palermo, Palermo, Italy</footnote>

<footnote n="5">⁵ Dipartimento di Scienze della Terra, Università di Firenze, Florence, Italy</footnote>
2007. I dati geochimici ottenuti su matrici vetrose, olivine, inclusioni vetrose ed embayment, e la loro integrazione con le misure dei gas vulcanici, suggeriscono differenze sistematiche fra l'esplosione maggiore studiata ed i parossismi in termini di sorgente magmatica, risalita e tempistiche pre-eruttive. Il contenuto di volatili nelle inclusioni vetrose mostra che la risalita del magma inizia da una sorgente leggermente più superficiale (~9.5 km sotto il livello del mare, s.l.m.) rispetto a quella dei parossismi (11.4 km s.l.m.) ed evidenziano l'attivazione di un'ulteriore zona di stazionamento a 5-6 km s.l.m. In accordo con le composizioni chimiche intermedie osservate nelle matrici vetrose, si suggerisce una risalita magmatica complessa, caratterizzata da condizioni di CO₂ tamponata nella sorgente profonda e da mescolamenti nella zona più superficiale. Modelli di diffusione Fe-Mg nelle olivine suggeriscono una perturbazione del sistema che inizia ~20-25 giorni prima dell'eruzione, in accordo con le tempistiche dei cambiamenti nel rapporto dei gas vulcanici CO₂/SO₂ osservati dal degassamento del plume. Queste tempistiche sono più brevi rispetto a quelle che precedono i parossismi (~4 mesi). L'insieme dei dati geochimici approfondisce le conoscenze sui processi che controllano le condizioni magmatiche stazionarie di Stromboli, così come sulla variabilità nell'intensità delle eruzioni, e ha delle implicazioni fondamentali per la previsione delle eruzioni vulcaniche.

**Keywords** Open-conduit volcano · Major explosions · Melt inclusions · Olivine · CO₂ flux · Stromboli

Introduction

In spite of intensive study, the controls on explosion inten-
sity, magma source and ascent conditions that determine
unexpected changes towards more violent events at open-
conduit systems, remain incompletely understood and
the subject of debate (cf. Pioli et al. 2014; Aiuppa et al.
2021, Vergniolle & Métrich 2022; Edmonds et al. 2022).
To aid in our understanding, we consider the triggering of
different eruption styles at open-conduit systems, using
Stromboli as our target, and focusing on the case of inter-
mediate intensity as associated with major explosions at
Stromboli. Stromboli is recognised as the archetype of an
open-conduit volcano due to its persistent activity over
the past ~1200 years. Its volcanic edifice forms the east-
ernmost island of the Aeolian archipelago, ~60 km N of
Sicily and 50 km W off the mainland coasts of Southern
Italy (Fig. 1a), rising 924 m above sea level (a.s.l.), while
the largest part of the edifice continues below sea level
(b.s.l.) to 1500–2600 m depth (Marani et al. 2004; Bosman
et al. 2009). Stromboli is one of the best-monitored vol-
canoes worldwide with a dense monitoring network that
has gradually expanded following the 2002–2003 eruption
crisis (Barberi et al. 2009; Ripepe et al. 2009; Rosi et al.
2013; Di Traglia et al. 2014; Giudicepietro et al. 2019;
Coppola et al. 2020). The monitoring network includes
multiple instruments deployed by the Istituto Nazionale di
Geofisica e Vulcanologia (INGV), by the Laboratorio di
Geofisica Sperimentale (LGS, Università di Firenze), and
by several Italian universities, allowing multi-parametric
and multi-disciplinary approaches to assess the state of the
volcano, and to eventually identify precursor signals and
improve early-warning systems (Ripepe et al. 2009, 2017,
2021; Giudicepietro et al. 2019, 2020). Data streamed by
the monitoring network have been employed for the clas-
sification of eruptive activity. In addition, the amplitude of

ground deformation measured by tiltmeters has been found
to scale with explosion intensity and has proved useful in
short-term eruption forecasting (Ripepe et al. 2021).

Stromboli's persistent activity comprises quiescent (passive) degassing, active degassing ("puffing"; Harris and Ripepe 2007; Ripepe et al. 2008; Aiuppa et al. 2010) and periodic (10 events/hour), low-intensity explosions (Ripepe et al. 2002, 2008; Delle Donne et al. 2006; Patrick et al. 2007; Rosi et al. 2000). These explosions, which typify the classical Strombolian activity, last a few seconds and eject bombs (10–10³ kg) and small tephra volumes (1–10 m³) to heights of 50–400 m above the active vents (Fig. 2a; Delle Donne et al. 2006; Patrick et al. 2007; Ripepe et al. 2008; Harris et al. 2013; Rosi et al. 2000; 2013; Bombrun et al. 2015). Before each of these "normal" explosions, a ground deformation of ~0.1 µrad is systematically detected by tiltmeters, typically starting ~150 s beforehand (Fig. 2b; Ripepe et al. 2021). More violent explosive eruptions occur occasionally and pose a significant threat to inhabitants and tourists on the island (Barberi et al. 1993; Rosi et al. 2006, 2013; Bertagnini et al. 2011; Giordano and De Astis 2020). These more violent eruptions span a wide range in magnitude and intensity, with the intermediate-scale events being termed "major explosions," while high-intensity events being referred to as "paroxysms." Major explosions (ME) occur on a yearly basis (Bevilacqua et al. 2020) and exhibit eruption intensities and characteristics intermediate between normal activity and paroxysmal activity; ground deformation measurable by tiltmeters starts 300 s prior to eruption onset and has an average amplitude of 0.8 µrad (Fig. 2b; Ripepe et al. 2021). Tephra volumes range between 10² and 10³ m³ and mass discharge rates range from 10⁴ to 10⁶ kg/s (Fig. 2a; Bertagnini et al. 1999; Andronico and Pistolesi 2010; La Felice and Landi 2011; Gurioli et al. 2013; Rosi et al. 2013; Pioli et al. 2014). Associated eruption plumes reach heights of few km with the erupted material affecting

<figure><img src="image_3.png" /></figure>
<figure><img src="image_4.png" /><figcaption>**Fig. 1** a Hillshade map of Stromboli volcano highlighting the position of the two populated villages to the NE and the SW. Inset shows position of the island relative to the Italian mainland and Sicily. Black arrow refers to direction of the photo in (b). b Eruption column of the 19 July 2020 major explosion as seen from the village of Strom-</figcaption></figure>

boli (photo courtesy M. Bitetto). c Scattered spatter bombs around the crater area (persons for scale in the background). Pumiceous, LP-type bombs are highlighted by red circles. d A mingled bomb, showing both HP and LP magma types with the scoriaceous HP component being marked in red

larger areas well beyond the crater area (Fig. 2a; Bertagnini et al. 1999; Andronico and Pistolesi 2010; Pioli et al. 2014).

The violent, high-intensity paroxysmal events have decadal frequencies, show a longer duration (few minutes), and often involve several vents simultaneously (Rosi et al. 2013; Bevilacqua et al. 2020). Associated mass discharge rates are several orders of magnitude larger than for normal activity (> $10^6$ kg/s; Fig. 2a; Rosi et al. 2006, 2013; Ripepe and Harris 2008; Pistolesi et al. 2011; Pioli et al. 2014). Associated ground deformation shows an average amplitude of ~10 µrad and can be registered ~600 s prior to eruption onset (Fig. 2b, Ripepe et al. 2021). Explosions are associated

with cannon-like detonations, leading to shockwaves and the ejection of meter-sized ballistic blocks to low elevations, where they pose a significant threat to the two villages of the island (Rosi et al. 2013; Fig. 1a). During the most recent paroxysms in 2003, 2007, and 2019, convective plumes with heights of up to 8 km were observed (Fig. 2a; Rosi et al. 2006, 2013; Pistolesi et al. 2011; Giordano and De Astis 2020). Tephra volumes range from $10^4–10^5$ m³ for the 2003, 2007, and 2019 paroxysms ("small-scale paroxysms (SSP)"; Bertagnini et al. 2011; Pistolesi et al. 2011; Rosi et al. 2013; Métrich et al. 2021), while paroxysms from previous centuries (e.g., 1930 CE, 1456 CE, sixteenth century) have tephra
<figure><img src="image_5.png" /></figure>

**Fig. 2** **a** Log-log plot showing typical plume heights (m) and erupted mass (kg) for normal Strombolian activity (black), major explosions (blue), and paroxysms (red). The range for the 19 July 2020 is shown by pink dashed lines given the observed plume height and the uncertainty relative to ejected mass. The large dots mark average values. Paroxysm data refers to 1919, 1930, 2003, 2007, and 2019 eruptions (Ponte 1919; Rosi et al. 2006; Bertagnini et al. 2011; Pistolesi et al. 2011; Pioli et al. 2014; Giordano and De Astis, 2020, Métrich et al. 2021). Modified from Métrich et al. (2021). **b** Semilog plot of

ground deformation with the stacking deformations of the 2003, 2007 and Jul and Aug 2019 paroxysms (red shaded area) and their average (bold red line) is represented together with the stacked deformations of major explosions (blue shaded area) and their average (bold blue line) and with regular (black line) Strombolian explosions (plot redrawn from and with data from Ripepe et al. 2021 and Genco and Ripepe 2010). The ground deformation for the 19 July 2020 eruption is shown in pink and is from this study

volumes at least one order of magnitude larger (“large-scale paroxysms (LSP)”; Rittmann 1931; Métrich et al. 2010, 2021; Bertagnini et al. 2011).

Numerous in-depth studies focused on the petrological and mineralogical characterisation of eruption ejecta from both normal and paroxysmal activity to constrain their magma ponding regions, and to improve the understanding of the timescales of pre-eruptive and syn-eruptive processes (Bertagnini et al. 2003; Landi et al. 2004; Francalanci et al. 2004, 2005; Petrone et al. 2018; Métrich et al. 2010, 2021; Pichavant et al. 2022). Geochemical characterisation of the ejecta of the Strombolian regime shows that normal activity is fed by a shallow

(<2 km b.s.l.), partially degassed ($H_2O < 0.6$ wt%; $CO_2 < 100$ ppm; $S < 1300$ ppm; $Cl < 2700$ ppm), high-potassium, shoshonitic melt (high $SiO_2 = 48.5–51.5$ wt%, low $CaO/Al_2O_3$ of ~0.48, and high $K_2O/Na_2O$ of ~1.2) that is commonly referred to as high porphyritic (HP), as it contains 50 vol% crystals (plagioclase $An_{60-70}$, clinopyroxene $Mg#_{0.70-0.90}$; $Fs_{12-15}$ and olivine $Fo_{70-73}$, Métrich et al. 2001, 2010; Landi et al. 2004, 2008; Francalanci et al. 2008; Bertagnini et al. 2011, 2013). Paroxysmal activity is characterised by the emission of highly vesicular, poorly crystallised, “golden” pumice, often mingled with the highly crystallised HP member (Métrich et al. 2001; Francalanci et al. 2004; Landi et al. 2004;

<figure><img src="image_6.png" /></figure>
Pistolesi et al. 2011; Pioli et al. 2014). The associated, low porphyritic (LP; < 5 vol% of crystals with mostly olivine Fo$_{82-91}$ and subordinate clinopyroxene Mg#$_{0.90-0.92}$; Fs$_{5-8}$) melt is a volatile-rich basalt (H$_2$O = 1.8–3.4 wt%; CO$_2$ = 890–1890 ppm; S = 1660–2250 ppm; Cl = 1660–2030 ppm) with lower SiO$_2$ (SiO$_2$ = 48–49.5 wt%) and K$_2$O/Na$_2$O (K$_2$O/Na$_2$O~0.7) and higher CaO/Al$_2$O$_3$ (0.68–0.81) compared to the HP magma (Métrich et al. 2001, 2005, 2010; Francalanci et al. 2004; Landi et al. 2004, 2009; Bertagnini et al. 2011). This bimodality in pyroclastic ejecta has been interpreted to sampling distinct portions of a vertically extended magmatic system, with the HP magma stored at shallow depths of ~2 km b.s.l., and the source of the more vesicular magma being placed at depths of 6–9 km b.s.l. down to the mantle-crust interface (Bertagnini et al. 2003; Métrich et al. 2010; Aiuppa et al. 2010; Pichavant et al. 2011). The bulk magma composition of HP and LP melts is similar, suggesting that the shallower magma is related to the deeper region by crystallisation through volatile loss and decompression (Métrich et al. 2001). While melt composition appears unrelated to eruption magnitude/intensity, a relationship between olivine composition and eruption intensity exists (Landi et al. 2004; Bertagnini et al. 2008; Métrich et al. 2010, 2021; Pichavant et al. 2022): olivine Fo$_{>85}$ is associated with paroxysmal events, whereas normal activity contains more evolved olivine with Fo$_{70-73}$.

Additional constraints on the processes that control magmatic system, degassing behaviour and eruptive activity common to open-conduit volcanoes are provided by gas analysis of the degassing plume (Burton et al. 2007; Allard 2010; Aiuppa et al. 2010, 2011, 2021). While normal activity is generally explained by the periodic rise and explosion of “gas slugs” from within the volcanic edifice (e.g., maximum depths of ~2 km b.s.l.; Burton et al. 2007; Allard 2010), different models exist to explain the more violent events. The “magma blob” and the “gas blob” models are widely accepted, each focusing on a different “end-member” scenario (e.g., Métrich et al. 2001, 2005; Bertagnini et al. 2003, 2008; Francalanci et al. 2004; Allard 2010; Aiuppa et al. 2011; Vergniolle and Métrich 2022). These two models are not mutually exclusive but coexist, with the predominance of one or the other being the factor that controls eruption intensity and dynamics: in the “magma blob” scenario, the intrusion of hot, volatile-rich magma into the deep system acts as the dominant force, while in the “gas blob” model, CO$_2$-rich gas in the sub-volcano system leads to foam accumulation at the top of the LP reservoir and rapid rise of CO$_2$-rich gas to the surface (Allard 2010; Aiuppa et al. 2011, 2021). Gas plume observations indicate that CO$_2$ degassing escalates over SO$_2$ degassing prior to paroxysms, indicating a deep trigger and the

limited involvement of the shallow magmatic system in paroxysmal activity (Bertagnini et al. 1999; Allard 2010; Aiuppa et al. 2010, 2021). Detailed studies of crystal zoning indicate that paroxysms are preceded by a magmatic recharge ~60–100 days before the eruption, agreeing with the timescales observed in the degassing plume (Aiuppa et al. 2010, 2021; Métrich et al. 2021).

## Rationale

General consensus exists for the volcanological (Pioli et al. 2014) and geophysical (Ripepe et al. 2021) classification of paroxysms. In contrast, the activity spectrum lying in between normal Strombolian activity and paroxysmal events (i.e., major explosions) remains difficult to classify. One of the limiting factors in the study of major explosions is that, while these eruptions are relatively frequent (several events a year; Bevilacqua et al. 2020), the dispersal of their associated ejecta is mostly limited to the upper part of the volcano, often hampering sampling operations. As such, the source region of these eruptions, and their associated pre-event dynamics, are poorly constrained. Aiuppa et al. (2021) proposed that accumulation of a CO$_2$-rich foam may act as the driver of these events and identifies a possible relationship between the erupted tephra volume and the size of a foam layer that would accommodate the CO$_2$ prior to eruption. Accordingly, paroxysms would be associated with a larger foam layer that can accommodate CO$_2$ over longer timescales (several months), while major explosions would be characterised by lower incubation times (weeks) and lesser CO$_2$ excess. While these incubation timescales for paroxysms (e.g., ~4 months) are independently estimated by petrology (Métrich et al. 2021), no such independent constrain is available for major explosions, limiting interpretation of data streamed by the monitoring network (Ripepe et al. 2017, 2021; Aiuppa et al. 2011, 2021).

Previous studies on ejecta from recent major explosions have yielded ambiguous results in terms of their characteristics: some explosions were seen to only erupt HP magma fragments and lithic material (Jan 2005; Landi et al. 2008; Calvari et al. 2012; Pioli et al. 2014) as the 7 September 2008 event, while some others contain juvenile material intermediate in composition between HP and LP members (e.g., 3 May, 8 and 24 Nov 2009; Andronico and Pistolesi 2010; La Felice and Landi 2011; Pioli et al. 2014). Vesicular LP magma was described only in the most powerful major explosions (e.g., 23 Aug 1998, 3 May, 8 and 24 Nov 2009 and 25 and 30 Jun 2010; Métrich et al. 2010; Aiuppa et al. 2011; Pioli et al. 2014). In addition, a variety of
trigger mechanisms have been proposed for these events. Allard (2010) suggests that major explosions are the result of the rise of gas slugs that are larger in size than for normal activity; a similarly shallow source for major and normal explosions has been suggested by Calvari et al. (2012) using evidence from the major explosion in September 2008. However, intermediate melt compositions are common to many major explosions, indicating variable degrees of interaction/homogenisation between HP and LP magmas (Andronico and Pistolesi 2010; La Felice and Landi 2011; Pioli et al. 2014), and requiring the involvement of the deep magmatic system, at least to some extent (as inferred for the 1998 and the 2009 major explosions; Landi et al. 2008; Andronico and Pistolesi 2010; La Felice and Landi 2011; Pioli et al. 2014). While differences in the relative proportions of outgassing and deep magma contribution, crystal content, and vesicularity suggest a ponding depth different from paroxysms (Pioli et al. 2014; Métrich et al. 2021), Pioli et al. (2014) propose that major explosions could result from a continuous hybridisation process driven by newly arriving LP magma in the upper part of the deep system. For paroxysms, research on H₂O and CO₂ contents in melt inclusions has constrained a 250–300 MPa pressure range for the associated ponding region (Bertagnini et al. 2003; Métrich et al. 2001; 2005; 2010, 2021). However, no comparable information on H₂O and CO₂ volatile contents in melt inclusions exists for major explosions (Métrich et al. 2001; 2010; 2021).

A major explosion recently occurred on the 19 July 2020 and was characterised by the emission of fresh, deep-derived, vesicular products (Fig. 1c, d). Such class of intermediate-scale events ejecting volatile-rich magma is particularly valuable as it offers an opportunity to gain information on the plumbing system. In the present study, we focus on the deep-derived (LP) component, while a characterisation of the HP and intermediate component of this eruption can be found in Landi et al. (2022). We present an extensive dataset that combines (i) compositional profiles in olivine crystals, (ii) analysis of olivine-hosted melt inclusions, and (iii) geochemical characterisation of matrix glasses, to constrain the timescales of the processes that preceded this eruption. This data is complemented by dissolved S, H₂O, and CO₂ in melt inclusions to discuss the dynamics of the deep magma contribution and to add further pressure (i.e., depth) constraints on the source region feeding these frequent, violent events. The integration of the geochemical data with gas plume time series from the monitoring network yields insights into the degassing behaviour and into precursor timescales preceding the explosive events that can be successfully detected by the monitoring network.

# Methods and analytical procedures

## Sample preparation and imaging

Highly vesicular pumice clasts were collected the day after the 19 July 2020 eruption (Fig. 1c, d). They were gently crushed and sieved for 10 min into grain size classes of 1 mm, 500 µm, 250 µm, and < 250 µm. Pumice fragments for matrix glass analysis were selected from the > 500 µm fraction, while olivine crystals were hand-picked from the 500 µm and the 250 µm size fractions. While free pyroxene and plagioclase phenocrystals are present, this research focuses on olivine, as they are assumed to contain the least disturbed melt inclusions representing the earliest stages of magma evolution (cf. Métrich et al. 2001).

Pumice fragments and crystals were embedded in epoxy plugs, polished and carbon coated. Electron backscattering (BSE) images were obtained using a Zeiss EVO MA 10 scanning electron microscope at the Istituto Nazionale di Geofisica e Vulcanologia (INGV) in Pisa (Italy). Olivine crystals were positioned parallel to crystallographic faces, when identified and polished until the centre was exposed. Core-rim profiles were selected following visual orientation based on BSE images and in agreement with the suggestions described in Shea et al. (2015).

## Electron microprobe analyses

Compositional profiles for diffusion modelling were obtained on 34 olivine crystals through electron microprobe analysis. In addition, point analyses were obtained next to melt inclusions and embayments to characterise host olivine compositions for an additional 13 olivine crystals. A first dataset was obtained on both 500 µm and 250 µm olivine crystals performed using a WDS at JEOL 8200 Super Probe at the University of Milan (Italy). Analyses were performed at 15 kV with a focused beam current of 5 nA and steps that varied from 15 µm for apparently homogeneous olivine to 5 µm for complex (zoned) crystals. Profiles and spot analyses were also acquired using a CAMECA SX five and a CAMECA SX 100-electron microprobe in Paris (Camparis Platform, Sorbonne-University, France). Analytical conditions and standard deviations for individual elements are outlined in the supplementary material (Table S1). Elemental analytical standard deviations are generally low (< 5%) with Fe and Mg being max. 4% and 1%, respectively.

Major element and S and Cl contents were determined for melt inclusions and matrix glass. Glass analyses were performed at the CAMECA SX five and the SX 100-electron microprobe in Paris with a defocused beam with 10 µm diameter and a 15 kV accelerating potential. Beam

<figure><img src="image_7.png" /></figure>
currents vary between 10 and 30 nA, depending on the analysed element (Supplementary Table S1). Analytical conditions and individual standard deviations are detailed in the supplementary material, with oxide analytical standard deviations being below 8% and accuracy being well below 5%, relative to an ALV glass and a diopside standard (Table S1). In most cases, multiple spots were analysed within one melt inclusion and were averaged subsequently. While overall melt inclusions and embayments in 33 olivine crystals were characterised by electron microprobe for their major element and S composition, a subset of melt inclusions and embayments in 13 olivine crystals was additionally analysed by FTIR spectroscopy. The complete datasets are reported in Table S4.

Equilibrium and potential post-entrapment crystallisation (PEC) were calculated using Petrolog v.3 (Danyushevsky and Plechov 2011). A K_D value of 0.3 (cf. Toplis 2005) was used at NNO + 1, following the model of Borisov and Shapkin (1990). For density and viscosity estimates, the build-in models of Lange and Carmichael (1987) and Bottinga and Weill (1972) were used, respectively. Values range between 0.27 and 12.41% with a mode at 7%, indicating limited PEC. MI compositions were chosen to not be corrected to maintain coherency between melt inclusions in host olivine crystals of variable complexity. Post-entrapment corrections could have introduced a bias in the data, whereas overall trends, such as ratios between major elements (CaO/Al₂O₃ and K₂O/Na₂O) are not affected by PEC, which affects Fe-Mg-Si in olivine-hosted melt inclusions (Métrich et al. 2010; Rose-Koga et al. 2021). Similarly, H₂O, CO₂, and S are diluted to the same degree, with any variation in their contents being close or within respective error bars. As recommended in Rose-Koga et al. (2021), the raw data, not normalised to 100 wt%, is reported in the Supplementary Table S4, together with the restored melt inclusion compositions after PEC corrections.

## Diffusion timescale calculations

To calculate diffusion timescales of observed zonation, Fe-Mg profiles in 18 olivine crystals were selected. In 7 of these, control profiles with different orientations were acquired to test for diffusion anisotropy and to exclude section effects (Shea et al. 2015). Diffusion timescales were obtained performing step-size modelling following Girona and Costa (2013). Modeling element diffusion in zoned crystals to estimate timescales assumes that crystallisation occurs in equilibrium with the surrounding melt (Girona and Costa 2013; Shea et al. 2015; Costa et al. 2020). Accordingly, zoning patterns indicate a change in the surrounding magmatic environment (e.g., crystallisation in a more primitive melt in case of reverse zoning; Costa et al. 2020). Following this, element diffusion begins between the

different concentration levels when the zone is formed and ends when the crystal is erupted. If the parameters that influence diffusivity can be constrained (e.g., temperature, pressure, oxygen fugacity, orientation in 3D crystals, major element composition), this allows a reconstruction of the associated timescales (Girona and Costa 2013). For the model calculations of the 2020 eruption, temperature and pressure constraints derive from experimental studies on Stromboli LP melt and were set to 1150 °C and 250 MPa (Di Carlo et al. 2006; Pichavant et al. 2011). Data points that had an error > 30 days were excluded.

## Fourier-Transform Infrared Spectroscopy

To quantify the H₂O and CO₂ contents in melt inclusions (MIs), 13 olivine crystals from the 500 µm grain size were selected under a stereomicroscope, in preparation for analysis by Fourier-Transform Infrared Spectroscopy (FTIR). Subsequently, EPMA analysis was performed on the same MIs to constrain major element composition and S (Supplementary Table S4). For FTIR preparation, the crystals were embedded in orthodontic resin and polished in order to expose MIs from both sides. However, due to the small size of melt inclusions (~20 µm) and to avoid breakage of the crystal wafer, not all melt inclusions could be exposed from both sides. They were exposed from one side and corrected for the olivine spectra during data processing following the approach outlined in Nichols and Wysoczanski (2007). H₂O and CO₂ analyses were performed by transmission IR spectroscopy using a Nicolet-Thermo Fischer IN10 FTIR at the Dipartimento di Scienze della Terra of Pisa (Italy). Volatile concentrations were calculated following the Beer-Lambert law as

$$C = \frac{100 \times Abs \times M}{\varepsilon \times \rho \times t}$$

Olivine thickness (t) was measured using a calibrated microscope stage and 10 × and 50 × zoom with relative deviations of ± 11%. For large glassy areas and embayments, several spectra were collected in different areas to ensure homogeneous volatile distribution. Water absorption (Abs) was obtained by measuring the peak height at 3535 cm⁻¹, relative to a linear baseline. Carbonate peaks are measured at 1520 cm⁻¹ and 1430 cm⁻¹, and peak deconvolution was performed using PeakFit following the subtraction of the spectra of an oxidised degassed glass.

The absorption coefficient for ε3535 was chosen as 64.3 L/mol per cm (Métrich et al. 2001). An average value of ε1520 = 362 L/mol per cm was used as carbon absorption coefficient. This value is an average from absorption coefficients calculated for individual MI compositions and averaged subsequently following the method outlined in Dixon and Pan (1995). Density (ρ) was chosen as 2.69 ± 0.02 g/
cm³, based on the value determined for basaltic glass fragments with 2.5 to 3.5 wt% H₂O (Métrich et al. 2001). The detection limit for the CO₂ that is measurable in glassy MIs strongly depends on the thickness of the analysed crystal wafer and varies accordingly from sample to sample. In the following, a uniform detection limit of 50 ppm is estimated for all MIs, in agreement with values stated for comparable glass compositions and analysis methods (Luhr 2001; Métrich et al. 2001). For melt inclusions and embayments that were exposed from one side only, the olivine absorption was measured between 1600 and 2000 cm⁻¹ and was used to correct the melt inclusion/embayment thickness (Nichols and Wysoczanski 2007; von Aulock et al. 2014). The olivine absorption used for thickness correction is indicated in Table S8 and no trend between calculated dissolved volatile contents and melt inclusion/embayment thickness was evident. Thickness measurements constitute the main error source for H₂O and CO₂ concentrations. Based on the propagation of the thickness measurement deviation, we estimate a maximum error of 12% for H₂O and CO₂ concentrations, which is comparable to an error of 10% estimated in other, comparable studies (Luhr 2001; Métrich et al. 2001).

Volatile saturation pressures were calculated assuming a temperature of 1150 °C (Di Carlo et al. 2006; Pichavant et al. 2011; Métrich et al. 2021) and using the MagmaSat model of Ghiorso and Gualda (2015) within the VesiCal (v.1.0.1) web-application (Iacovino et al. 2021). To ensure data comparability, previously published saturation pressures for the 2007 eruption at Stromboli (Métrich et al. 2005, 2010, 2021), estimated with the model of Papale et al. (2006), were recalculated. A detailed discussion of the models used for saturation pressure calculations is reported in the supplementary material.

# Results

## The 19 July 2020 major explosion

The eruption occurred at 03:00:45 UTC on 19 July 2020 during a period of low activity. Based on the daily integrated evaluation of the parameters from the monitoring system delivered by LGS (Laboratorio di Geofisica Sperimentale – University of Florence) to the Italian Civil Protection, the overall activity of the volcano during 2020 was fluctuating between low and medium levels, with sporadic, weeks-long periods of higher intensity between January and April.

Visual and thermal observations of the 19 July explosion by the LGS monitoring system showed three eruption pulses, reaching a column height of ~1 km (Fig. 1b). Eruptive products were dispersed around the crater area down to an elevation of ~500 m a.s.l., setting fire to the vegetation. Spatter bombs up to 1 m in diameter fell as a scatter shower

around the crater area and were characterised by mingling of a deep, LP-type and a scoriaceous, HP-type component (Fig. 1c, d).

Geophysical parameters associated to the event showed a ground deformation of 3.5 µrad that started 4 min prior to the event (Fig. 2b). Furthermore, a seismic velocity 4.0×10⁻⁴ m/s and a seismic amplitude of 1×10⁻⁴ m were registered, while the infrasonic network recorded a maximum pressure of 250 Pa at 450 m distance and localized the event in the central/SW sectors. Based on geophysical parameters used to describe the variability of the activity at Stromboli, the eruption was classified as a “major explosion”; particularly, ground deformation recorded along the volcano slope by tiltmeters, placed this eruption at the boundary between “typical” major explosions (0.8 µrad on average, onset of the inflation ramp 300 s prior to eruption) and “typical” paroxysms (10 µrad, 600 s prior to eruption; Fig. 2b; Ripepe et al. 2021).

On a longer timescale, the 19 July 2020 eruption was preceded by a large increase of the volcanic gas plume CO₂/SO₂ ratio (median at 20 to 30 molar ratio) and of the volcanic CO₂ flux (median at ~2000 tons/day) (data from Aiuppa et al. 2021). This burst of deeply sourced (Aiuppa et al. 2021) CO₂-rich gas started abruptly on the 29 June 2020, i.e., ~20 days before the blast, persisted until the major explosion and then rapidly vanished (see Fig. 3).

## Olivine chemical profiles and zoning

A total of 47 olivine crystals were characterised in terms of composition and zoning patterns, with detailed compositional profiles being acquired on 34 crystals. No correlation between composition and grain size is evident. Olivine compositions range between Fo₆₆ and Fo₈₇, predominantly showing values of Fo₈₀-₈₇ (Fig. 4a). Zoning patterns are predominantly reverse (53%), followed by homogeneous compositions (26%) and subordinate normally zoned crystals (21%). Multiple population modes can be distinguished when considering compositional frequencies, defining five groups (Fig. 4a, b): in the first group, the prevailing compositions have a mode of Fo₈₅—measured either in homogeneous crystals or in cores of zoned crystals—to a maximum of Fo₈₇. This group represents the least evolved olivine compositions of this eruption. Within the first group, rims tend to be slightly more evolved compared to core compositions. Associated chemical profiles show a distinct core plateau and an exponential increase or decrease towards the rim for reverse and normal zoning, respectively (Fig. 4c). Some of them are chemically homogeneous (Fig. 4d).

A second composition group comprises crystals Fo₈₂-₈₄ that occur both as cores and rims (Fig. 4b). Where they form the core, they show systematically
**Fig. 3** Volcanic gas plume time-series throughout year 2020 (redrawn from data in Aiuppa et al. 2021). Large star indicates the 19 July 2020 major explosion; small stars indicate three smaller major explosions in November 2020. Results are presented as centred 7-day averaged normalised distributions of a CO₂/SO₂ molar ratios in the plume and b CO₂ fluxes. Dark red colour tones illustrate temporal variations of the median of the frequency distributions, while light blue colour tones correspond to frequency bins in which measurements are 30 to 60% less frequent than the median. An obvious acceleration in CO₂-rich gas release is evident since late June 2020

<figure><img src="image_8.png" /></figure>

reverse zoning patterns and occur surrounded by a group 1 rim (Fo₈₄-₈₇; Fig. 4f). In the third group, compositions range from Fo₇₉-₈₂ and exist exclusively as crystal cores, whereas a fourth group with Fo₇₆-₇₉ is found exclusively as rim compositions. Finally, the fifth population (Fo₆₆-₇₆) is subordinate and corresponds to inherited cores.

Olivine-hosted MIs are mostly associated with group 1 compositions (Fo₈₅, yellow line Fig. 4a). Less frequently, MIs are hosted by crystals with Fo₈₃ and Fo₈₁ (groups 2 and 3; Fig. 4a). Among the 33 analysed crystals hosting MIs and embayments, 48.5% show homogenous compositions, whereas 51.5% are reversely (33.3%) or normally (18.2%) zoned, respectively.

As a whole, the compositional spectrum of olivine in the July 2020 pumice clasts agree well with that outlined for previous major explosions (Pichavant et al. 2011, 2022; Métrich et al. 2021). The least evolved compositions (Fo₈₅-₈₇) overlap with the range that was previously described for small-scale paroxysms (e.g., the 2003, 2007 events; Métrich et al. 2021). No olivine with compositions Fo>₈₇, as is common in large paroxysms (Métrich et al. 2010, 2021), are identified.

## Matrix glass and melt inclusion compositions

### Major elements

The 19 July 2020 matrix glasses plot in the domain of the crystal-poor basalts (LP) that were reported in the literature for previous eruptions (e.g., 1998 ME, 2003 SSP, 2007 SSP, and 1456 LSP; Bertagnini et al. 2003; Métrich et al. 2005, 2010), having CaO/Al₂O₃ = 0.48–0.72 ± 4% and K₂O/Na₂O = 0.67–1.34 ± 11% (Fig. 5a). Interestingly, although the July 2020 matrix glass compositions cluster in the range typical of the LP magma, they also extend compositionally towards the more evolved HP member, with some analyses having K₂O/Na₂O ratios close to 1 (Fig. 5a, Table S4). Overall, matrix glass compositions highlight a wide variability, comprising intermediate compositions similar to those that have been reported for other major explosions (i.e., 3 May, 8 and 24 Nov 2009, La Felice and Landi 2011). Contrastingly, matrix glass compositions observed in previous small-scale paroxysms of 2003 and 2007 as well as in the 1998 major explosion show a limited variability both for matrix glass as well as olivine compositions (Fig. 5, Fig. S3).
<figure><img src="image_9.png" /><figcaption>**Fig. 4** **a** Frequency diagram (bin size Fo mol% = 1) showing compositions for analysed olivine crystals marked by a solid line. Fo contents of MI hosting olivine are shown in yellow, while olivine compositions next to embayments are marked by a dashed line. Olivine group compositions are marked by numbers with Gr. 1- Fo<sub>84-87</sub>, Gr. 2- Fo<sub>82-84</sub>, Gr. 3- Fo<sub>79-82</sub>, Gr. 4- Fo<sub>76-79</sub>, and Gr. 5- Fo<sub>&lt;76</sub>. Olivine compositional ranges are reported for normal Strombolian activity (SA), major explosions (ME), and paroxysms (P). **b** Compositional zoning pattern within different olivine crystals. Arrow direction indicates compositional change from core to rim. Grey arrows mark zon-</figcaption></figure>

ing profiles that have an exponential decrease or increase in the Fo contents towards the rim for normally and inversely zoned crystals, respectively. Red arrows mark complex zoning, where more multiple compositions were observed, while diamonds highlight homogeneous composition. Examples of **c** a normally zoned olivine with an exponential pattern; **d** a compositionally homogeneous crystal; **e** an inversely zoned olivine with core and rim characterised by two distinct compositional plateaus; **f** multiple compositions. Corresponding backscattered electron images of olivine crystals are also shown (white scale bar corresponds to 100 µm)

Glassy melt inclusions display equivalent or less evolved compositions than matrix glass, with ratios of CaO/Al₂O₃ = 0.6–0.96 and K₂O/Na₂O = 0.54–0.81 (Figs. 5a, 6a, b, Table S5). These compositions agree with previously obtained results for the 1998 major explosion and for the 2003 and 2007 small paroxysms but extend towards (and beyond) the 1456 LSP dataset. Embayments show compositions that overlap with matrix glasses but are generally more evolved than MIs, with CaO/Al₂O₃ = 0.57–0.68, K₂O/Na₂O = 0.61–0.86 (Fig. 6b). Melt inclusion CaO/Al₂O₃ ratios are variable, possibly due to minor clinopyroxene crystallisation in the melt that leads to decreasing CaO/

Al₂O₃ ratios, while K₂O/Na₂O ratios remain relatively constant (e.g., La Felice and Landi 2011). Host olivine compositions remain constant at ~Fo₈₅, except for those MIs that are associated with host compositions Fo₇₉-₈₄ and cluster in the embayment domain (Fig. 6a, b).

## Dissolved volatiles

Sulphur and chlorine in matrix glass show values of S = 0.008–0.117 wt% (± 26% average analytical error, where higher S contents have an error of ± 7%, while values lower than 0.018 have a maximum error of ± 65%) and
<figure><img src="image_10.png" /></figure>

**Fig. 5** Glass compositions for matrix glass and MIs. The compositional field for matrix glass for high porphyricity melts (HP) is reported in grey, following data on scoriae from ¹Métrich et al. (2001, 2010) and Landi et al. (2008), while the compositional field for low porphyricity (LP) matrix glass drawn for the 1998 major explosion, the 2003 and 2007 small paroxysms and the 1456 CE large paroxysm is from ²Métrich et al. (2005, 2010). Glass and MI compositions for the large-scale paroxysm 1456 CE (LSP) are taken from Bertagnini et al. (2003) and Métrich et al. (2021). **a** K₂O/Na₂O and CaO/Al₂O₃ ratios for matrix glass, melt inclusions, and embayments. Matrix glasses are marked by filled symbols, while MIs are represented by empty symbols. New data from this study is marked by black diamonds for the 2020 major explosion matrix glass and red crosses for MIs and embayments from the 2020 eruption. Average errors are marked by an error bar with K₂O/Na₂O=11% and CaO/Al₂O₃=4%. **b** Ratios of S/Cl and CaO/Al₂O₃ in melt inclusions, embayments, and matrix glass. Symbol coding is the same as for **a** and average errors are reported as S/Cl=9% and CaO/Al₂O₃=4%. For matrix glass compositions, the average error for S/Cl is ±28% with S and Cl values being close to the detection limit. The respective error in the error bar is marked by a dashed line

Cl=0.075–0.169 wt% (±9%). While S/Cl ratios show a wide range (0.066–0.733 ± 28%, with ± 10% for high S values and maximum error of ± 66% for low S values), the majority of data points cluster between S/Cl=0.066 and S/Cl=0.277 (Fig. 5b). Overall, these compositional ranges agree well with previously obtained compositions for matrix glass volatiles obtained in LP melts from the 1998 major

explosion and the 2003, 2007, and 1456 paroxysms (Bertagnini et al. 2003; Métrich et al. 2005, 2010, 2021).

Glassy MIs are characterised by variable sulphur contents of S=0.082–0.227 wt% (± 6%) whereas chlorine contents show limited variability with Cl=0.157–0.208 wt% (± 7%). Compared to matrix glass, S/Cl ratios in MIs show a larger variability, extending to higher values up to a range of 0.52–1.14 (± 9%) (Fig. 5b). While no clear relationship between volatile composition of MIs and composition of the surrounding host olivine can be recognised, MIs associated with groups 2 and 3 hosts (Fo₇₀₋₈₄) tend to exhibit slightly lower S contents (Fig. 6c). It is noteworthy that the S contents observed in the 2020 major explosion MIs compare well with those obtained for the 2007 and the 1456 paroxysms (Fig. 6c) as well as with the 1998 and the 2003 events (Fig. S4). At almost constant sulphur compositions, 2020 MIs span a wide range in H₂O contents (Fig. 6c), showing slightly lower values than melt inclusions from both 2007 and the 1456 events. Sulphur contents in 2020 embayments show lower values between 0.039 and 0.153 (± 6%) at almost constant H₂O contents (Fig. 6c) and constant CaO/Al₂O₃ ratios of ~0.64 (Fig. S4). In comparison, 2007 MIs and embayments highlight a simultaneous decrease in S and H₂O with embayments, showing consistently lower S and H₂O values (Fig. 6c).

Dissolved H₂O and CO₂ contents show no systematic relationship with melt inclusion size and span a wide range with H₂O=1.53–3.64 wt% and CO₂=580–1866 ppm, with primary mode compositions of H₂O=2.8 wt% and CO₂=1500 ppm (Fig. 7a). Glass embayments show lower volatile contents with a primary mode at H₂O=2.4 wt% and CO₂=750 ppm and display a more limited compositional range (H₂O=1.72–2.4 wt%; CO₂=450–1639 ppm, Fig. 7a) than MI. Water and carbon dioxide contents previously described for the 2003 and 2007 small paroxysms show a comparable range (primary modes at 2.8 wt% H₂O and 1750 and 1250 ppm CO₂ for 2007 MIs, 2.6 wt% H₂O and 750 ppm CO₂ for the 2007 embayments, and 2.2 wt% H₂O and 1750 ppm CO₂ for 2003 MIs; Métrich et al. 2005, 2010). It is interesting to note the bimodality between melt inclusions and embayments in the 2007 data, with two clear distinct fields showing lower volatile compositions (H₂O=0.77–2.95 wt% and CO₂=343–1622 ppm) for the embayments, while melt inclusions outline higher volatile amounts (H₂O=2.64–3.78 wt% and CO₂=1126–1704 ppm, Fig. 7a).

The calculated volatile pressures for the July 2020 encompass a wide pressure range, both for MIs with a range of 95–333 MPa, a main peak at 250–300 MPa, and a secondary peak at 150–200 MPa, as well as for the embayments that show pressures between 121 and 296 MPa, with a prominent peak around 175 MPa. Assuming an average
**Fig. 6** Melt inclusion compositions for the 2020 major explosion in comparison with literature data. **a** Relationship between olivine host composition (Fo composition in mol%) and MI and embayment major elements (CaO/Al₂O₃). Error bar is Fo=2% and CaO/Al₂O₃=4%. For comparison, MI data is also shown for the large-scale 1456 CE paroxysm (LSP, solid black outline, Bertagnini et al. 2003; Métrich et al. 2021) and for small-scale paroxysms eruptions 2007, 2003 (SSP, dashed) and the 1998 major explosion (ME, dotted outline, Métrich et al. 2005, 2010, 2021). **b** Major element ratios K₂O/Na₂O vs CaO/Al₂O₃ in MIs and embayments. MI compositions are shown in red: Fo₈₄.₈₇, orange Fo₈₂.₈₄, yellow Fo₇₉.₈₂, and blue Fo₇₆.₇₉. Error bar is 11% and 4% for K₂O/Na₂O and CaO/Al₂O₃, respectively. **c** MI and embayment sulphur and H₂O compositions. Symbols are as in insets **a** and **b**, while MI data is reported from literature for the 1456 CE paroxysm (LSP, dark circles) and for the 2007 paroxysm (SSP, grey squares). Embayment data is available only for the 2007 paroxysm (empty grey squares). The error bar is 6% for S and a maximum of 12% for H₂O

crustal density of 2700 kg/m³ pressures of 150, 175, and 275 MPa correspond to depths of 4.8, 5.7, and 9.5 km b.s.l., respectively.

## Timescales from diffusion modelling

Diffusion timescales from Fe–Mg compositional profiles were calculated for different zoning patterns in 18 olivine crystals (total of 25 Fe–Mg profiles, 7 of which were used as control profiles; Fig. 8a) using the DIPRA model of Girona and Costa (2013). Timescales range between a minimum of 3 days and a maximum of 2 months, generally being shorter than 40 days (Fig. 8a), with two frequency modes at 5–10 and 20–25 days (dotted blue line in Fig. 8b). Timescales show no clear relationship with the different zoning pattern but vary slightly for the different zoning types, with exponential trends showing timescales of 3 days to 5 weeks, while olivine crystals with two compositional plateaus indicate 1–4 weeks. For olivine with complex, multiple zoning, the outermost zoning indicates timescales of less than 3 weeks. Timescales calculated from Fe–Mg diffusion profiles in olivine crystals for 5 different paroxysms (July and August 2019 SSP, 1456 LSP, sixteenth century LSP, 1930 LSP; Métrich et al. 2021) range between 5–100 days for SSP and 2–60 days for LSP. Interestingly, considering the frequency distribution of paroxysm timescales (dashed red line in Fig. 8b, with only the shorter interval ≤40 days reported here, while the timeframe 60–100 days is not shown; cf. Métrich et al. 2021), these highlight a single mode at <5 days.

## Discussion

### Source depth of the 19 July 2020 eruption

The comprehensive presented geochemical characterisation of the tephra ejected during the 19 July 2020 eruption

<figure><img src="image_11.png" /></figure>

includes olivine, MI, and matrix glass analyses. Olivine compositions and zoning patterns show large variability; however, the prevailing mode at Fo₈₅.₈₇, which represents the least evolved olivine compositions of the 2020 eruption, distinguishes it from large-scale paroxysms that contain more magnesian (Fo>87; Fig. 4, Fig. S3) olivine crystals (Métrich et al. 2010, 2021). Matrix glasses of the 19 July 2020 eruptive products also cover a large

<figure><img src="image_12.png" /></figure>
**Fig. 7** **a** Dissolved H₂O and CO₂ contents in MIs and embayments. ► Literature data for the small-scale 2003 and 2007 paroxysms is taken from Métrich et al. (2005, 2010). Data from this study is shown as red crosses for MIs, and as red empty diamonds for embayments. Melt inclusions and embayments that are distributed along the x-axis have CO₂ contents that are below the detection limit and are reported to have CO₂=0 for representative purposes only. A maximum error is estimated at 12%. **b** Shows the same data for the 2020 eruption as in **a** and also shows the evolution of H₂O and CO₂ contents for different scenarios with a starting composition of H₂O=2.9 wt% and CO₂=2500 ppm. Solid black line=open-system degassing; black dotted line and dashed line=evolution where the melt is in equilibrium with 1 wt% and 3 wt% of an exsolved gas phase. The blue arrow highlights the evolution of melt composition assuming CO₂ buffering, while the red arrow indicates a trend with melt subjected to mixing with melts with lower volatile contents. **c** Frequency diagram (bin size 50 MPa) with saturation pressures calculated based on H₂O and CO₂ compositions, assuming a temperature of 1150 °C and using MagmaSat (Ghiorso and Gualda 2015) and VesiCal (Iacovino et al. 2021). Saturation pressures for the small-scale paroxysm 2007 were recalculated with the same models for consistency, using the volatile data from Métrich et al. (2010). Depth is indicated in km below sea level (b.s.l.) and is calculated from lithostatic depths assuming an averaged density for surrounding rocks as 2700 kg/m³

compositional spectrum including both typical LP com-
positions, as observed in previous eruptions, as well as
"intermediate" compositions between the HP and LP fields
(Fig. 5). MIs and embayments record similarly variable
compositions, which generally compare well with those
observed in past major explosions (e.g., 1998) and small-
scale paroxysms (2003, 2007), but also suggest that erup-
tion source and ascent dynamics may reflect a more com-
plex history (see below).

Further constraints to allow better discussion of the
ponding depth of major explosions are gained from
total volatile (CO₂ + H₂O) pressures. It is well known
that absolute values for pressure and depth estimates
are strongly dependent on the used solubility model and
therefore need to be treated with caution (see the supple-
mentary material for a more detailed discussion). Previ-
ously published pressures for the 2007 paroxysm from
Métrich et al. (2010) were recalculated with the same
solubility model used for the 2020 major explosion data
(Fig. S2, Table S3), to ensure the comparability of the
different datasets. H₂O and CO₂ contents identify two
primary pressure modes at 275 MPa and 175 MPa for the
2020 eruption that correspond to the pressures recorded
by MIs and embayments, respectively (Fig. 7c). For an
average crustal density of 2700 kg/m³, this corresponds
to melt inclusion entrapment depths at 9.5 km b.s.l.. MIs
from the 2007 paroxysm, on the other hand, record higher
total volatile contents (Fig. 7a), yielding relatively deeper
(recalculated) entrapment depths of 11.4 km b.s.l.. The
wider range in 2020 MI volatile contents points towards
a secondary peak in MI pressure distributions at lower
pressures of 175 MPa (Fig. 7a, c). These pressures agree

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>
<figure><img src="image_16.png" /></figure>

**Fig. 8** **a** Timescales calculated from Fe–Mg compositional profiles using diffusion modelling in DIPRA (Girona and Costa 2013). A number of 7 control profiles was analysed on selected crystals: their data is reported by grey crosses. The uncertainty associated with diffusion timescale modelling for the 2020 major explosion is shown by error bars. **b** Timescales (days) with arrow direction marking the zoning direction from core to rim and arrow length showing the associated olivine compositions. The associated frequency distribution for the 2020 major explosion is shown by a dotted blue line and highlights two modes at 5–10 and 20–25 days, respectively. The frequency distribution for the timescales following the Fe–Mg diffusion modelling in Métrich et al. (2021) is reported by a red dashed line

with a peak at 150–175 MPa as indicated by 2020 and 2007 embayment contents, and point to a second, shallower pressure region, where olivine crystallisation and glass entrapment took place during magma ascent and evolution.

Different from the narrow pressure distribution outlined by MI volatile contents from the 2007 paroxysm, which suggest a single entrapment event at 325 MPa (Fig. 7c), melt inclusions of the 2020 eruption highlight that glass entrapment took place at depths indicated by pressures ranging from 300 to 120 MPa. This leads to two possible interpretations: (1) the melt rises from a deeper (10 km b.s.l.) zone and stagnates at 5–6 km b.s.l., prior to the final ascent; (2) the ascending melt from the 10 km depth encounters a stationary, shallower melt at 5–6 km b.s.l., where self-mixing occurs, enhancing volatile exsolution from the volatile-rich deep magma, eventually triggering the eruption. We argue that the presence of a magma ponding zone at 5–6 km b.s.l. (Fig. 9a), where olivine Fo<sub>80-84</sub> crystallisation and self-mixing occur, is justified by the compositional range of matrix glass (Fig. 5a), and the continuous range in MIs and embayments chemistry that are hosted mostly in reversely zoned

and highlights a mode at <5 days. The analysed paroxysms include Jul and Aug 2019 SSP, and the 1456 CE, the undated 16th and the 1930 LSP. Colour coding and nomenclature for zoning pattern is the same as in Fig. 4 with zoning patterns in grey representing an exponential increase for inversely zoned or decrease in Fo-composition towards the rim for normally zoned crystals. Zoning patterns showing two distinct compositional plateaus for core and rim are reported in black and zoning pattern with multiple compositions are shown in red. For the latter, only compositions that have been used for diffusion modelling are shown (i.e., the outermost two concentration levels). Triangles mark olivine crystals with compositional variation that is restricted within one compositional group

or homogeneous crystals. These observations lead us to propose a continuous influx and hybridisation of this shallow zone by a melt that is in equilibrium with Fo<sub>84-87</sub>, pointing towards the second scenario.

## Insights into the ascent dynamics of the 19 July 2020 eruption

Under open-system conditions, CO₂ is early exsolved from the melt while H₂O remains dissolved until lower pressures, while closed-system conditions are consistent with a scenario in which the melt is in equilibrium with variable amounts of an exsolved gas phase. The exsolved gas phase is composed mostly of CO₂, decreasing the CO₂ melt composition, while H₂O exsolution becomes increasingly dominant at lower pressures. Sulphur exsolution tends to be coupled with H₂O exsolution and occurs generally at pressures < 150 MPa (Métrich et al. 2010, 2021; Lesne et al. 2011a,b; Aiuppa et al. 2021; Pichavant et al. 2013; 2022, Laiolo et al. 2022). Based on the trends observed in MI and embayment volatile contents, we highlight several observations: first, the majority of MIs and embayments follow closed-system degassing in equilibrium with ~2

<figure><img src="image_17.png" /></figure>
<figure><img src="image_18.png" /><figcaption>**Fig. 9** **a** Schematic overview of the 19 July 2020 major explosion, summarising the main characteristics such as input (1), processes in the reservoir (2) and during ascent (3) and changes that were observed by the monitoring network prior to eruption (4). Insets **b**, **c**,</figcaption></figure>

and **d** represent a simplified conceptual model of Stromboli's plumbing system, where different source regions are activated for **b** low-intensity, **c** high-intensity major explosions, and **d** paroxysms and are highlighted in colour

wt% of a gas phase (cf. with dotted and dashed lines in Fig. 7b, which mark 1 wt% and 3 wt% exsolved gas phase, respectively), but are inconsistent with open-system degassing (solid line in Fig. 7b). Comparable conditions were previously discussed for the 2007 paroxysm by Métrich et al. (2010) and Aiuppa et al. (2010), who described the melt to be in equilibrium with 2.4 wt% CO₂-rich gas phase. Second, the clustering of few MIs and embayments along a mixing trend (red line “mixing” in Fig. 7b) corroborates the previously discussed hypothesis of the interaction between a more volatile-rich and a more degassed melt at a shallower (5–6 km b.s.l.) depth. Third, contrary to the expected coupled degassing behaviour of S and H₂O, MIs show a trend of relatively constant S, at markedly decreasing H₂O levels (Fig. 6c). A similar trend of near constant, but slightly increasing CO₂ contents at decreasing water contents can be observed in several MIs (trend “CO₂ buffering” in inset Fig. 7b), suggesting that the analysed 2020 melt inclusions testify partial H₂O loss in the melt. MI dehydration has been related both to post-entrapment hydrogen diffusion through olivine (e.g., Lloyd et al. 2013) as well as to magma buffering by a CO₂-rich phase of deeper origin (Blundy et al. 2010; Métrich et al. 2010). The trend of water decreasing at nearly constant sulphur compositions in MIs diverges from that recorded by embayments, which highlights decreasing sulphur (Fig. 6c). Sulphur compositions of ~1500 ppm have been demonstrated to mark the sulphur saturation limit in hydrous and oxidised magmas

at Stromboli, both experimentally (Lesne et al. 2011b) and from MI analyses (Métrich et al. 2010). We interpret these trends to indicate that most melt inclusions might have not reached S exsolution pressures (Lesne et al. 2011a, b). This would imply that processes leading to magma dehydration occurred at pressures > 150 MPa, which agrees with a lack of relationship between H₂O and CaO/Al₂O₃ ratios or with K₂O and embayments showing lower sulphur contents (Fig. 6c). The high CO₂ values in 2020 melt inclusions (Fig. 7a, b) concur with previously described high CO₂ values of 1689 ppm, observed in melt inclusions containing carbonate crystals within bubbles (Métrich et al. 2001). Also, embayment compositions with high CO₂, in combination with low CaO/Al₂O₃ ratios, confirm a continuous flux in the system, in agreement with the 2 wt% gas phase in equilibrium with the melt as previously described.

In summary, MI and embayment volatile data indicate multiple processes that occurred in the storage region, involving the interaction between melts with different volatile contents at depths between 5–6 km b.s.l., but also highlighting a buffering process of the melt by a CO₂-rich phase within the deep storage region at 9.5 km b.s.l.. It is worth highlighting that while there are important differences in the eruption source depth and the ascent and degassing behaviour of the 19 July 2020 major explosion compared to paroxysmal eruptions, the underlying factors are strikingly similar, testifying the steady-state conditions that control the volcanic activity at Stromboli.
## Role of the CO₂-phase in major explosions and paroxysms: towards a comprehensive scheme

Prior to the 19 July 2020 major explosion, degassing occurs at a pressure range of ~150–175 MPa, as indicated by embayment H₂O and CO₂ compositions (Fig. 7a,c). This pressure range corresponds to a depth of 5–6 km b.s.l. and is shallower than the top of the LP deep reservoir (Fig. 9; Métrich et al. 2021). The CO₂ gas-buffered signature witnessed by MIs suggests that the melt coexists with a large fraction of exsolved gas in this pressure range (dashed and dotted curves in Fig. 7b). Therefore, MIs of the 2020 eruption are fully consistent with mechanisms (Allard 2010; Aiuppa et al. 2010, 2021) in which Stromboli's more violent explosions are triggered by accumulation of gas bubbles to form a gas-rich “foam” that ultimately collapses once a volumetric threshold is reached. It is interesting to note that this gas foam collapse process is similar for both major explosions as well as paroxysms, although occurring at different depths and involving different volumetric thresholds for the foam layer (Aiuppa et al. 2021). In fact, while for the paroxysms, an accumulation zone is proposed at relatively greater depths, we show that a shallower (5–6 km b.s.l.) zone is activated during the 2020 major explosion, extending the potential source region for intermediate and LP-type magma to shallower depths than previously known. Although the 2020 explosion represents an extreme case for major explosions, almost at the lower boundary of paroxysms (Fig. 2a, b), it also provides the opportunity to discuss the similarities between eruptions of different intensity and magnitude and to extend it to major explosions. This adds crucial insights into Stromboli's plumbing system but also raises the question about the underlying mechanism, and whether at Stromboli eruption intensity and magnitude might follow a continuous scale, potentially being related to the activation of different regions (depths) of the LP system.

Similar to a distinction between small-scale and large-scale paroxysms (Andronico et al. 2008; Métrich et al. 2010, 2021), a distinction between “ideal” end-members of variable eruption intensity and magnitude can be approached by discussing the dominant processes (Fig. 9).

Previous studies have suggested that eruption intensity and magnitude are controlled by a delicate balance between outgassing and the deep system contribution, with only the more violent events being influenced by deeper processes (e.g., Pioli et al. 2014). Accordingly, low-intensity, major explosions that are just above normal activity represent the lower end-member of intermediate-scale eruptions. They are characterised by ejection of high proportions of HP magma and lithics, with scarce to no LP magma emission (e.g., 9 Jan 2005, Andronico et al. 2008, Landi et al. 2008; 7 Sept 2008, Calvari et al. 2012), and are dominated by more transient and shallow outgassing processes (Fig. 9b). In line with

the processes discussed above, the melt is degassed, with melt inclusion compositions being similar to those of the surrounding matrix glass (Landi et al. 2009; Métrich et al. 2010). Intermediate matrix glass compositions can be identified, but the LP component is negligible, and the influx of deeply sourced gas less significant, as indicated by more subtle CO₂ flux variations (compare the CO₂ flux changes prior to the 19 July 2020 event with the very minor change preceding the November 2020 series of smaller events; Fig. 3 and Aiuppa et al. 2021). With increasing intensity, but still within the major explosions field, the influence from the deep system increases (Fig. 9c), as witnessed by those major explosions that contain a larger amount of LP magma (e.g., 23 August 1998, Bertagnini et al. 1999; 8 and 24 November 2009, La Felice and Landi 2011; Pioli et al. 2014) and are accompanied by a larger change in the CO₂ flux (Fig. 3). These higher-intensity major explosions, such as the one studied here, share geophysical characteristics that are similar to small-scale paroxysms such as those occurred in 2003, 2007, and 2019 (Fig. 2). From a geochemical point of view, olivine compositions of Fo₈₅-₈₇ evidence comparable equilibrium crystallisation to those of small-scale paroxysms. However, several characteristics differ significantly from the latter, including a wider glass matrix compositional range, a larger interval of MI compositions, and a complex entrapment pressure distribution, suggesting the activation of at least two different magmatic ponding regions. This is also consistent with different timescales of magma ascent, with paroxysms showing an earlier first perturbation of the deep magmatic system up to 4 months prior to eruption (Aiuppa et al. 2021; Métrich et al. 2021). For the 2020 eruption, timescales are significantly shorter, with Fe–Mg diffusion profiles in olivine testifying a major system perturbation ~20–25 days prior to eruption onset (Fig. 8), agreeing with a change in the CO₂/SO₂ ratio in the degassing plume (Fig. 3; Aiuppa et al. 2021). We ascribe this major system perturbation to the arrival of a CO₂-rich gas phase in the deep storage zone and the self-mixing between volatile-rich and volatile-poorer magma in the upper part of the LP reservoir and the shallow ponding region (Fig. 9a). This, in turn, leads to increasing gas accumulation and a concomitant increase in CO₂/SO₂ gas flux, as observed in the degassing plume (Fig. 3). A continuous deep influx and hybridisation is consistent with a wide range in melt compositions, resulting in olivine crystallisation and MI and embayment entrapment at a wide pressure range between 120 and 300 MPa, while contributing to gas foam accumulation (Fig. 9a). The shorter timescales associated with the first system perturbation compared to those described for paroxysms, are consistent with the model discussed in Aiuppa et al. (2021) that relates precursory timescales with foam layer size. Accordingly, for a major explosion, a combination of system geometry and relatively higher gas/melt ratios lead to shorter incubation

<figure><img src="image_19.png" /></figure>
times compared to paroxysms that involve larger magma and gas volumes. Moreover, for paroxysms, an earlier first system perturbation, several months before eruption onset (Métrich et al. 2021; Aiuppa et al. 2021) is coherent with a direct ascent of magma from greater depths, as testified by the Mg-rich olivine (Fo<sub>89-90</sub>) in the largest paroxysms (Fig. 9d; Métrich et al. 2021). This deepest component is inferred to be sourced from the crust-mantle interface, where low V<sub>p</sub> and V<sub>s</sub> are interpreted to represent a zone of partial melt (e.g., at 16 km b.s.l., ~450 MPa; Panza et al. 2007; Pichavant et al. 2011) and is absent in major explosions, as evidenced by the lack of Mg-rich-olivine (Fo<sub>89-90</sub>).

Shorter timescales of 2–10 days were observed both in this work as well as for paroxysms (Fig. 8b) and are interpreted to represent the very final magma ascent in the conduit (Métrich et al. 2021). The frequency mode of <5 days highlighted by the paroxysm data is slightly shorter than the 5–10-day mode outlined by the 2020 major explosion (compare dashed red line with dotted blue line in Fig. 8b), while the ascent occurs from greater depths for paroxysms compared to major explosions, as confirmed by embayment pressures (Fig. 7b). In other words, the final ascent for paroxysms occurs over shorter time of <5 days despite starting from greater depth possibly suggesting higher ascent velocities, while major explosions are characterised by a more gradual, ~5–10 days-long slower ascent from shallower depths.

## Implications for monitoring and eruption forecasting

Accurate eruption modelling and forecasting needs to start from correct interpretation of magma source and dynamics. While the final journey of magma ascent through the conduit for paroxysms at Stromboli has been successfully modelled and is used for early-warning systems by exploiting the mechanical response of the edifice (i.e., deformation of the upper part of the cone; Giudicepietro et al. 2020; Ripepe et al. 2021), the modelling of major explosions still requires important advancements. We have shown here that the source for the most violent major explosions lies in the deep reservoir but involves slightly shallower depths compared to that of paroxysms. While paroxysms, both of small- and large-scale, rise directly from the deep LP reservoir, with absolute, recalculated depths of 11.4 km b.s.l., major explosions have more complex magma dynamics: glass and olivine chemistry, and the volatile contents of MIs and embayments suggest a potential activation of a shallower zone placed at 5–6 km b.s.l.. This shallower magmatic source could be either the result of a fast evolution of the deeper magma during ascent, or a specific magmatic region with those characteristics, activated by the fast rise of the deeper magma. Whatever the case, this region should be

used as a source for modelling of this category of explosive events, considering that magma/gas volumes involved are to be considered smaller if compared with paroxysms of any scale. In line with this, potential precursory signals that can be revealed by the monitoring system can be expected to be of similarly smaller magnitude. Given the large range of characteristics of major explosions (intensity, erupted material, areal extent), petrological studies cannot be applied to all this variability, particularly to those of smaller scale, characterised by the absence of suitable juvenile material. For classification purposes, quantitative parametrization by means of geophysical parameters, such as the amplitude and timing of ground deformation (Fig. 2b), may allow a discrimination of most of the major events from the normal activity. Intermediate compositions described for the 2020 major explosion fit well within the framework of a previously proposed relationship between olivine compositions and eruption intensity/magnitude (Métrich et al. 2021) and provide a promising outlook for future, more detailed petrological classification schemes. In this perspective, the geochemical signal of the CO₂/SO₂ ratio and the CO₂ flux from the degassing plume (Aiuppa et al. 2021) and its correlation with magmatic processes in the deep system from petrological studies (e.g., Métrich et al. 2021) revealed to be an effective tool to document perturbations of the deep magmatic system over timescales of days to months. For the majority of cases, this appeared to add crucial information to the long-term forecasting of major explosions occurrence. In a more general scheme, the presented comprehensive dataset on a violent major explosion adds further crucial insights in terms of magma storage depths and pre-eruptive timescales that will be integrated within the database that provides the framework for short- and long-term eruption forecasting.

## Conclusions

The 19 July 2020 major explosion provides a unique opportunity to improve the existing understanding of Stromboli's plumbing system and its pre-eruptive processes, as it allows integration of geochemically-derived petrological data on the deep, volatile-rich magmatic component with time series for the degassing plume. Using our new data, we can summarise four main outcomes and implications. First, total dissolved volatile contents (H₂O + CO₂) in MIs suggest a magma ponding depth of 9.5 km b.s.l., which is slightly shallower than the recalculated paroxysm source depth of 11.4 km b.s.l.. Overall, this also implies that processes that control intermediate and paroxysmal activity within Stromboli's deep plumbing system take place at greater depths than previously envisioned. Second, melt inclusions and embayments witness a difference in ascent behaviour, involving multiple depths (where olivine crystallisation is
favoured) and continuous ascent, as evidenced by the wide pressure range covered by melt inclusions and embayments. The final ascent of the volatile-rich magma begins from a depth of 5–6 km b.s.l. (or a source pressure of 150–175 MPa). Third, the analysis of matrix glass, olivine compositions as well as melt inclusions and embayments outline compositions intermediate between typical LP and HP melts, which is witnessing processes of hybridisation with magma in equilibrium with olivine Fo<sub>84-87</sub>. This is also consistent with self-mixing between melts with different volatile contents within the upper part of the deep reservoir. Contrary to paroxysmal eruptions, no involvement of deep-seated magma is recorded in the mineralogy of the basaltic juvenile clasts. Fourth, timescales calculated from diffusion profiles in olivine crystals indicate a major system perturbation ~20–25 days before eruption onset, represented by the arrival of a CO<sub>2</sub>-rich phase in the deep magma ponding zone and the self-mixing of magmas with slightly different volatile contents at shallower depths. These timescales are consistent with a sudden increase observed in the CO<sub>2</sub>/SO<sub>2</sub> gas flux, as measured in the surface degassing plume by the monitoring network. Moreover, they are significantly shorter than the 60–100 days observed for paroxysms, consistent with a first system perturbation due to magma ascent from greater depths (Métrich et al. 2021). In line with the observed differences in ascent behaviour, the timescales for final ascent prior to eruption onset involve 5–10 days, compared to <5 days for paroxysms, suggesting a faster magma ascent for the latter.

The 2020 eruption is fully consistent with a bottom-up process that includes magma self-mixing and the tapping of a shallow zone at 5–6 km b.s.l., where CO₂ accumulation leads to foam formation. Despite the differences for the geophysical signals and the geochemical data between the 19 July 2020 major explosion and paroxysms, it appears that the main controlling processes are similar, such as buffering by a CO₂ phase, magma hybridisation, and gas foam accumulation. The shorter incubation times observed for major explosions are consistent with the previously suggested scaling relationship between foam layer size and system perturbation time (Aiuppa et al. 2021). The question whether eruptive activity at Stromboli follows a continuous magnitude and intensity range needs to be addressed in future multi-parametric studies, integrating geophysical monitoring data, field data, gas time-series, and petrological analysis of the deep-derived component to investigate a comprehensive scaling scheme. Finally, the observation of the reactivation of a shallower area within the intermediate-deep reservoir has important implications for eruption forecasting. The final stagnation (and degassing) depth of 5–6 km needs to be taken into account when identifying precursory signals and modelling early-warning timescales.

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1007/s00445-023-01647-y.

**Acknowledgements** The authors thank N. Rividi and M. Fialin at the Sorbonne-University, France (Camparis) and A. Risplendente at Università degli Studi di Milano for help with microprobe analyses. M. Bitetto is thanked for the photo in Fig. 1. This work was supported by the Italian Civil Protection in the framework of the agreements between the Universities of Florence and Pisa and between INGV and the University of Pisa. We thank J. Hammer and A. Harris for editorial handling and helpful suggestions on the manuscript as well as M. Myers, R. Spallanzani, M. J. Rutherford, and an anonymous reviewer for their detailed comments and suggestions. Finally, we would like to thank F. Van Wyk de Vries for help with handling the paper through the EdMan system.

**Funding** Open access funding provided by Università di Pisa within the CRUI-CARE Agreement.

**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://creativecommons.org/licenses/by/4.0/.

## References

Aiuppa A, Bertagnini A, Métrich N, Moretti R, Di Muro A, Liuzzo M, Tamburello G (2010) A model of degassing for Stromboli volcano. Earth Planet Sci Lett 295(1–2):195–204. https://doi.org/10.1016/j.epsl.2010.03.040

Aiuppa A, Burton M, Allard P et al (2011) First observational evidence for the CO2-driven origin of Stromboli’s major explosions. Solid Earth 2:135–142. https://doi.org/10.5194/se-2-135-2011

Aiuppa A, Bitetto M, Delle Donne D et al (2021) Volcanic CO2 tracks the incubation period of basaltic paroxysms. Sci Adv. https://doi.org/10.1016/j.earscirev.2006.08.004

Allard P (2010) A CO2-rich gas trigger of explosive paroxysms at Stromboli basaltic volcano, Italy. J Volcanol Geoth Res 189:363–374. https://doi.org/10.1016/j.jvolgeores.2009.11.018

Andronico D, Pistolesi M (2010) The November 2009 paroxysmal explosions at Stromboli. J Volcanol Geoth Res 196:120–125. https://doi.org/10.1016/j.jvolgeores.2010.06.005

Andronico D, Corsaro RA, Cristaldi A, Polacci M (2008) Characterizing high energy explosive eruptions at Stromboli volcano using multidisciplinary data: an example from the 9 January 2005 explosion. J Volcanol Geoth Res 176(4):541–550. https://doi.org/10.1016/j.jvolgeores.2008.05.011

Barberi F, Rosi M, Sodi A (1993) Volcanic hazard assessment at Stromboli based on review of historical data. Acta Vulcanol 3:173–187

Barberi F, Civetta L, Rosi M, Scandone R (2009) Chronology of the 2007 eruption of Stromboli and the activity of the Scientific Synthesis Group. J Volcanol Geoth Res 182(3–4):123–130. https://doi.org/10.1016/j.jvolgeores.2008.09.019
Bertagnini A, Coltelli M, Landi P et al (1999) Violent explosions yield new insights into dynamics of Stromboli volcano. Eos (Washington DC) 80:2–8. https://doi.org/10.1029/99EO00415

Bertagnini A, di Roberto A, Pompilio M (2011) Paroxysmal activity at Stromboli: lessons from the past. Bull Volcanol 73:1229–1243. https://doi.org/10.1007/s00445-011-0470-3

Bertagnini A, Metrich N, Landi P, Rosi M (2003) Stromboli volcano (Aeolian Archipelago, Italy): an open window on the deep-feeding system of a steady state basaltic volcano. Journal of Geophysical Research 108. https://doi.org/10.1029/2002JB002146

Bertagnini A, Métrich N, Francalanci L, et al (2008) Volcanology and Magma Geochemistry of the Present-Day Activity: Constraints on the Feeding System. In: Calvari S, Inguaggiato S, Ripepe M, Rosi M (eds) The Stromboli Volcano, An integrated study of the 2002–2003 eruption. AGU Geophysical Monograph. Series, Washington DC, pp 19–38. https://doi.org/10.1029/182GM04

Bertagnini A, Métrich N, Francalanci L, et al (2013) Volcanology and magma geochemistry of the present-day activity: constraints on the feeding system. 19–37. https://doi.org/10.1029/182gm04

Bevilacqua A, Bertagnini A, Pompilio M et al (2020) Major explosions and paroxysms at Stromboli (Italy): a new historical catalog and temporal models of occurrence with uncertainty quantification. Sci Rep 10:1–18. https://doi.org/10.1038/s41598-020-74301-8

Blundy J, Cashman KV, Rust A, Witham F (2010) A case for CO2-rich arc magmas. Earth Planet Sci Lett 290(3–4):289–301. https://doi.org/10.1016/j.epsl.2009.12.013

Borisov AA, Shapkin AI (1990) A new empirical equation rating Fe3+/Fe2+ in magmas to their composition, oxygen fugacity, and temperature. Geochem Int 27:111–116

Bosman A, Chiocci FL, Romagnoli C (2009) Morpho-structural setting of Stromboli volcano revealed by high-resolution bathymetry and backscatter data of its submarine portions. Bull Volcanol 71(9):1007–1019. https://doi.org/10.1007/s00445-009-0279-5

Bottinga Y, Weill DF (1972) The viscosity of magmatic silicate liquids: a model for calculation. Am J Sci 272:438–475. https://doi.org/10.2475/ajs.272.5.438

Burton M, Allard P, Muré F (1979) la Spina A (2007) Magmatic gas composition reveals the source depth of slug-driven Strombolian explosive activity. Science 317:227–230. https://doi.org/10.1126/science.1141900

Bombrun M, Harris A, Gurioli L, Battaglia J, Barra V (2015) Anatomy of a Strombolian eruption: Inferences from particle data recorded with thermal video. Journal of Geophysical Research: Solid Earth 120(4) 2367–2387 https://doi.org/10.1002/2014JB011556

Calvari S, Büttner R, Cristaldi A et al (2012) The 7 September 2008 Vulcanian explosion at Stromboli volcano: multiparametric characterization of the event and quantification of the ejecta. J Geophys Res: Solid Earth 117:1–17. https://doi.org/10.1029/2011JB009048

Coppola D, Laiolo M, Cigolini C et al (2020) Thermal remote sensing for global volcano monitoring: experiences from the MIROVA system. Front Earth Sci 7(January):1–21. https://doi.org/10.3389/feart.2019.00362

Costa F, Shea T, Ubide T (2020) Diffusion chronometry and the timescales of magmatic processes. Nat Rev Earth Environ 1(4):201–214. https://doi.org/10.1038/s43017-020-0038-x

Danyushevsky LV, Plechov P (2011) Petrolog3: integrated software for modeling crystallization processes. Geochem Geophys Geosyst 12:Q07021. https://doi.org/10.1029/2011GC003516

Delle Donne D, Lacanna G, Marchetti E, Ripepe M, & Ulivieri G (2006). Monitoring explosive volcanic activity using thermal images, Stromboli volcano, Italy. AGU Fall Meeting Abstracts, V43B-1795.

Di Carlo I, Pichavant M, Rotolo SG, Scaillet B (2006) Experimental crystallization of a high-K arc basalt: the golden pumice, Stromboli volcano (Italy). J Petrol 47(7):1317–1343. https://doi.org/10.1093/petrology/egl011

Di Traglia F, Nolesini T, Intrieri E, Mugnai F, Leva D, Rosi M, Casagli N (2014) Review of ten years of volcano deformations recorded by the ground-based INSAR monitoring system at Stromboli volcano: a tool to mitigate volcano flank dynamics and intense volcanic activity. Earth Sci Rev 139:317–335. https://doi.org/10.1016/j.earscirev.2014.09.011

Dixon JE, Pan V (1995) Determination of the molar absorptivity of dissolved carbonate in basanitic glass. Am Miner 80:1339–1342. https://doi.org/10.2138/am-1995-11-1224

Edmonds M, Liu EJ, Cashman KV (2022) Open-vent volcanoes fuelled by depth-integrated magma degassing. Bull Volcanol 84(3). https://doi.org/10.1007/s00445-021-01522-8

Francalanci L, Tommasini S, Conticelli S (2004) The volcanic activity of Stromboli in the 1906–1998 AD period: Mineralogical, geochemical and isotope data relevant to the understanding of the plumbing system. J Volcanol Geoth Res 131(1–2):179–211. https://doi.org/10.1016/S0377-0273(03)00362-7

Francalanci L, Davies GR, Lustenhouwer W, Tommasini S, Mason PRD, Conticelli S (2005) Intra-grain Sr isotope evidence for crystal recycling and multiple magma reservoirs in the recent activity of Stromboli volcano, Southern Italy. J Petrol 46(10):1997–2021. https://doi.org/10.1093/petrology/egi045

Francalanci L, Bertagnini A, Métrich N, Renzulli A, Vannucci R, Landi P, Del Moro S, Menna M, Petrone C M., & Nardini I (2008). Mineralogical, geochemical, and isotopic characteristics of the ejecta from the 5 April 2003 paroxysm at Stromboli, Italy: Inferences on the Preeruptive Magma Dynamics BT - The Stromboli Volcano: An Integrated Study of the 2002–2003 Erup tion. The Stromboli Volcano: An Integrated Study of the 2002–2003 Erup- tion 182:331–345. http://www.earth-prints.org/handle/2122/4635%0Apapers3://publication/uuid/7CD14DA7-D40F-4A25-B2CE-FD4275AD4284

Genco R, Ripepe M (2010). Inflation-deflation cycles revealed by tilt and seismic records at Stromboli volcano. Geophysical Research Letters 37:(12). https://doi.org/10.1029/2010GL042925

Ghiorso MS, Gualda GAR (2015) An H2O–CO2 mixed fluid saturation model compatible with rhyolite-MELTS. Contrib Miner Petrol 169:1–30. https://doi.org/10.1007/s00410-015-1141-8

Giordano G, De Astis G (2020) The summer 2019 basaltic Vulcanian eruptions (paroxysms) of Stromboli. Bull Volcanol 83:1. https://doi.org/10.1007/s00445-020-01423-2

Girona T, Costa F (2013) DIPRA: a user-friendly program to model multi-element diffusion in olivine with applications to timescales of magmatic processes. Geochem Geophys Geosyst 14:422–431. https://doi.org/10.1029/2012GC004427

Giudicepietro F, López C, Macedonio G et al (2020) Geophysical precursors of the July–August 2019 paroxysmal eruptive phase and their implications for Stromboli volcano (Italy) monitoring. Sci Rep 10(1):1–16. https://doi.org/10.1038/s41598-020-67220-1

Giudicepietro F, Calvari S, Alparone S, et al (2019). Integration of ground-based remote-sensing and in situ multidisciplinary monitoring data to analyze the eruptive activity of stromboli volcano in 2017–2018. Remote Sensing, 11(15). https://doi.org/10.3390/rs11151813

Gurioli L, Harris AJL, Colò L, Bernard J, Favalli M, Ripepe M, Andronico D (2013) Classification, landing distribution, and associated flight parameters for a bomb field emplaced during a single major explosion at Stromboli, Italy. Geology, 41(5), 559–562. https://doi.org/10.1130/G33967.1

Harris A, Delle Donne D, Dehn J, Ripepe M, Worden, AK (2013) Volcanic plume and bomb field masses from thermal infrared camera
imagery. Earth and Planetary Science Letters, 365, 77–85. https://doi.org/10.1016/j.epsl.2013.01.004

Harris A, Ripepe M (2007) Temperature and dynamics of degassing at Stromboli. J Geophys Res: Solid Earth 112(3):1–18. https://doi.org/10.1029/2006JB004393

Iacovino K, Matthews S, Wieser PE et al (2021) VESIcal Part I: an open-source thermodynamic model engine for mixed volatile (H2O-CO2) solubility in silicate melts. Earth and Space Sci 8:1–55. https://doi.org/10.1029/2020EA001584

La Felice S, Landi P (2011) The 2009 paroxysmal explosions at Stromboli (Italy): magma mixing and eruption dynamics. Bull Volcanol 73:1147–1154. https://doi.org/10.1007/s00445-011-0502-z

Laiolo M, Delle Donne D, Coppola D et al (2022) Shallow magma dynamics at open-vent volcanoes tracked by coupled thermal and SO2 observations. Earth and Planetary Science Letters 594:117726. https://doi.org/10.1016/j.epsl.2022.117726

Landi P, D'Oriano C, Petrelli M, Nazzari M, Andronico D (2022). Inferences on the magmatic plumbing system at Stromboli volcano (Italy) from trace element geochemistry of matrix glasses and minerals in different types of explosive eruptions. Contrib Mineral Petrol 1. https://doi.org/10.1007/s00410-022-01962-1

Landi P, Metrich N, Bertagnini A, Rosi M (2004) Dynamics of magma mixing and degassing recorded in plagioclase at Stromboli (Aeolian Archipelago, Italy). Contrib Miner Petrol 147:213–227. https://doi.org/10.1007/s00410-004-0578-y

Landi P, Métrich N, Bertagnini A, Rosi M (2008) Recycling and “re-hydration” of degassed magma inducing transient dissolution/crystallization events at Stromboli (Italy). J Volcanol Geoth Res 174:325–336. https://doi.org/10.1016/j.jvolgeores.2008.02.013

Landi P, Corsaro RA, Francalanci L, Civetta L, Miraglia L, Pompilio M, Tesoro R (2009) Magma dynamics during the 2007 Stromboli eruption (Aeolian Islands, Italy): mineralogical, geochemical and isotopic data. J Volcanol Geoth Res 182(3–4):255–268. https://doi.org/10.1016/j.jvolgeores.2008.11.010

Lange RA, Carmichael ISE (1987) Densities of Na2O-K2O-CaO-MgO-FeO-Fe2O3-Al2O3-TiO2-SiO2 liquids: New measurements and derived partial molar properties. Geochim Cosmochim Acta 51:2931–2946. https://doi.org/10.1016/0016-7037(87)90368-1

Lesne P, Kohn SC, Blundy J, Witham F, Botcharnikov RE, Behrens H (2011b) Experimental simulation of closed-system degassing in the system basalt-H2O-CO2-S-Cl. J Petrol 52(9):1737–1762. https://doi.org/10.1093/petrology/egr027

Lesne P, Scaillet B, & Pichavant M (2011a). The H2O solubility of alkali basaltic melts: an experimental study. 133–151. https://doi.org/10.1007/s00410-010-0588-x

Lloyd AS, Plank T, Ruprecht P, Hauri EH, Rose W (2013) Volatile loss from melt inclusions in pyroclasts of differing sizes. Contrib Miner Petrol 165:129–153. https://doi.org/10.1007/s00410-012-0800-2

Luhr JF (2001) Glass inclusions and melt volatile contents at Parícutin volcano, Mexico. Contrib Miner Petrol 142:261–283. https://doi.org/10.1007/s004100100293

Marani M, Gamberi F, Bortoluzzi G, Carrara G, Ligi M, Penitenti D (2004) Tyrrhenian sea bathymetry. In: Marani MP, Gamberi F, Bonatti E (eds) From seafloor to deep mantle: architecture of the Tyrrhenian backarc basin. APAT, Mem Descr Carta Geologica d’Italia 44 Plate 1–2

Métrich N, Bertagnini A, Landi P, Rosi M (2001) Crystallization driven by decompression and water loss at Stromboli volcano (Aeolian Islands, Italy). J Petrol 42:1471–1490. https://doi.org/10.1093/petrology/42.8.1471

Métrich N, Bertagnini A, di Muro A (2010) Conditions of magma storage, degassing and ascent at Stromboli: New insights into the volcano plumbing system with inferences on the eruptive dynamics. J Petrol 51:603–626. https://doi.org/10.1093/petrology/egp083

Métrich N, Bertagnini A, Pistolesi M (2021) Paroxysms at Stromboli volcano (Italy): source, genesis and dynamics. Front Earth Sci 9:1–17. https://doi.org/10.3389/feart.2021.593339

Métrich N, Bertagnini A, Landi P, et al (2005) Triggering mechanism at the origin of paroxysms at Stromboli (Aeolian Archipelago, Italy): the 5 April 2003 eruption. Geophys Res Lett 32. https://doi.org/10.1029/2004GL022257

Nichols ARL, Wysoczanski RJ (2007) Using micro-FTIR spectroscopy to measure volatile contents in small and unexposed inclusions hosted in olivine crystals. Chemical Geology 242(3–4):371–384. https://doi.org/10.1016/j.chemgeo.2007.04.007

Panza GF, Peccerillo A, Aoudia A, Farina B (2007) Geophysical and petrological modelling of the structure and composition of the crust and upper mantle in complex geodynamic settings: the Tyrrhenian Sea and surroundings. Earth Sci Rev 80(1–2):1–46. https://doi.org/10.1016/j.earscirev.2006.08.004

Papale P, Moretti R, Barbato D (2006) The compositional dependence of the saturation surface of H2O+CO2 fluids in silicate melts. Chem Geol 229:78–95. https://doi.org/10.1016/j.chemgeo.2006.01.013

Patrick MR, Harris AJL, Ripepe M et al (2007) Strombolian explosive styles and source conditions: insights from thermal (FLIR) video. Bull Volcanol 69:769–784. https://doi.org/10.1007/s00445-006-0107-0

Petrone CM, Braschi E, Francalanci L, Casalini M, Tommasini S (2018) Rapid mixing and short storage timescale in the magma dynamics of a steady-state volcano. Earth Planet Sci Lett 492:206–221. https://doi.org/10.1016/j.epsl.2018.03.055

Pichavant M, Pompilio M, D'Oriano C, Di Carlo I (2011) Petrography, mineralogy and geochemistry of a primitive pumice from Stromboli: implications for the deep feeding system. Eur J Mineral 23:499–517. https://doi.org/10.1127/0935-1221/2011/0023-2109

Pichavant M, Di Carlo I, Pompilio M, le Gall N (2022) Timescales and mechanisms of paroxysm initiation at Stromboli volcano, Aeolian Islands. Italy Bullet Volcanol. https://doi.org/10.1007/s00445-022-01545-9

Pichavant M, Di Carlo I, Rotolo SG, Scaillet B, Burgisser A, Le Gall N, Martel C. (2013). Generation of CO2-rich melts during basalt magma ascent and degassing. 545–561. https://doi.org/10.1007/s00410-013-0890-5

Pioli L, Pistolesi M, Rosi M (2014) Transient explosions at open-vent volcanoes: the case of Stromboli (Italy). Geology 3–6. https://doi.org/10.1130/G35844.1

Pistolesi M, Delle Donne D, Pioli L et al (2011) The 15 March 2007 explosive crisis at Stromboli volcano, Italy: assessing physical parameters through a multidisciplinary approach. J Geophys Res: Solid Earth 116:1–18. https://doi.org/10.1029/2011JB008527

Ponte G (1919) La catastrofica esplosione dello Stromboli. R Accad Naz Lincei 28:89–94

Ripepe M, Harris AJL (2008) Dynamics of the 5 April 2003 explosive paroxysm observed at Stromboli by a near-vent thermal seismic and infrasonic array. Geophysical Research Letters 35(7). https://doi.org/10.1029/2007GL032533

Ripepe M, Harris AJL, Carniel R (2002) Thermal, seismic and infrasonic evidences of variable degassing rates at Stromboli volcano. J Volcanol Geoth Res 118:285–297. https://doi.org/10.1016/S0377-0273(02)00298-6

Ripepe M, Delle Donne D, Harris A, Marchetti E, Ulivieri G (2008) Dynamics of Strombolian activity. Washington DC Am Geophys Union Geophys Monograph Series 182:39–48

Ripepe M, Delle Donne D, Lacanna G et al (2009) The onset of the 2007 Stromboli effusive eruption recorded by an integrated geophysical network. J Volcanol Geoth Res 182:131–136. https://doi.org/10.1016/j.jvolgeores.2009.02.011

Ripepe M, Pistolesi M, Di C et al (2017) Forecasting effusive dynamics and decompression rates by magmatic model at
open-vent volcanoes. Sci Rep 7:1–9. https://doi.org/10.1038/
s41598-017-03833-3

Ripepe M, Lacanna G, Pistolesi M, et al (2021) Ground deformation reveals the scale-invariant conduit dynamics driving explosive basaltic eruptions. Nature Communications 12:. https://doi.org/
10.1038/s41467-021-21722-2

Rittmann A (1931) Der Ausbruch des Stromboli am 11 September
1930. Zeitschrift Für Vulkanologie 14:47–77

Rose-Koga EF, Bouvier AS, Gaetani GA, Wallace PJ, Allison CM, Andrys
JA, Angeles de la Torre CA, Barth A, Bodnar RJ, Bracco Gartner
AJJ, Butters D, Castillejo A, Chilson-Parks B, Choudhary BR, Cluzel
N, Cole M, Cottrell E, Daly A, Danyushevsky LV, ... Zhou T (2021).
Silicate melt inclusions in the new millennium: a review of recom-
mended practices for preparation, analysis, and data presentation.
Chem Geol 570. https://doi.org/10.1016/j.chemgeo.2021.120145

Rosi M, Bertagnini A, Landi P (2000) Onset of the persistent activity at Stromboli volcano (Italy). Bull Volcanol 62:294–300. https://doi.org/10.1007/s004450000098

Rosi M, Bertagnini A, Harris AJL et al (2006) A case history of paroxysmal explosion at Stromboli: timing and dynamics of the April 5, 2003 event. Earth Planet Sci Lett 243:594–606. https://doi.org/10.1016/j.epsl.2006.01.035

Rosi M, Pistolesi M, Bertagnini A et al (2013) Chapter 14 Stromboli volcano, Aeolian Islands (Italy): present eruptive activity and hazards. Geol Soc, London, Memoirs 37:473–490. https://doi.org/10.1144/M37.14

Shea T, Costa F, Krimer D, Hammer JE (2015) Accuracy of time-scales retrieved from diffusion modeling in olivine: a 3D perspective. Am Miner 100(10):2026–2042. https://doi.org/10.2138/am-2015-5163

Toplis MJ (2005) The thermodynamics of iron and magnesium partitioning between olivine and liquid: criteria for assessing and predicting equilibrium in natural and experimental systems. Contrib Miner Petrol 149(1):22–39. https://doi.org/10.1007/s00410-004-0629-4

Vergniolle S, Métrich N (2022) An interpretative view of open-vent volcanoes. Bull Volcanol 84(9):1–45. https://doi.org/10.1007/s00445-022-01581-5

von Aulock FW, Kennedy BM, Schipper CI, Castro JME, Martin D, Oze C, Watkins JM, Wallace PJ, Puskar L, Bégué F, Nichols ARL., & Tuffen H (2014) Advances in Fourier transform infrared spectroscopy of natural glasses: from sample preparation to data analysis. Lithos, 206–207(Supplement C), 52–64. https://doi.org/10.1016/j.lithos.2014.07.017