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

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 CO2 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 CO2/SO2 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.


Introduction
In spite of intensive study, the controls on explosion intensity, magma source and ascent conditions that determine unexpected changes towards more violent events at openconduit systems, remain incompletely understood and the subject of debate (cf. Pioli et al. 2014;Aiuppa et al. 2021, Vergniolle & Métrich 2022Edmonds 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 intermediate 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 easternmost 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 volcanoes 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;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(Ripepe et al. , 2017(Ripepe et al. , 2021Giudicepietro et al. , 2020. Data streamed by the monitoring network have been employed for the classification 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).
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. 2006Rosi et al. , 2013Pistolesi 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. 2006Rosi et al. , 2013Pistolesi et al. 2011;Giordano and De Astis 2020). Tephra volumes range from 10 4 -10 5 m 3 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(e.g., CE, 1456 CE, sixteenth century) have tephra 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-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 volumes at least one order of magnitude larger ("large-scale paroxysms (LSP)"; Rittmann 1931;Métrich et al. 2010Métrich et al. , 2021Bertagnini 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 syneruptive processes Landi et al. 2004;Francalanci et al. 2004Francalanci et al. , 2005Petrone et al. 2018;Métrich et al. 2010Métrich et al. , 2021Pichavant 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 2 O < 0.6 wt%; CO 2 < 100 ppm; S < 1300 ppm; Cl < 2700 ppm), highpotassium, shoshonitic melt (high SiO 2 = 48.5-51.5 wt%, low CaO/Al 2 O 3 of ⁓0.48, and high K 2 O/Na 2 O 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. 2001Métrich et al. , 2010Landi et al. 2004Landi et al. , 2008Francalanci et al. 2008;Bertagnini et al. 2011Bertagnini et al. , 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; 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, 2003eruptions (Ponte 1919Rosi et al. 2006;Bertagnini et al. 2011;Pistolesi et al. 2011;Pioli et al. 2014;Giordano andDe 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 andGenco andRipepe 2010). The ground deformation for the 19 July 2020 eruption is shown in pink and is from this study 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(Métrich et al. , 2005Francalanci et al. 2004;Landi et al. 2004Landi et al. , 2009Bertagnini 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 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. 2010Métrich et al. , 2021Pichavant 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. 2010Aiuppa et al. , 2011Aiuppa et al. , 2021. While normal activity is generally explained by the periodical 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. 2001Métrich et al. , 2005Bertagnini et al. 2003Bertagnini et al. , 2008Francalanci 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. 2011Aiuppa et al. , 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. 2010Aiuppa et al. , 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(Aiuppa et al. , 2021Mé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(Ripepe et al. , 2021Aiuppa et al. 2011Aiuppa et al. , 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 2009and 25 and 30 Jun 2010Mé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 2 O and CO 2 contents in melt inclusions has constrained a 250-300 MPa pressure range for the associated ponding region Métrich et al. 2001;2005;2010. However, no comparable information on H 2 O and CO 2 volatile contents in melt inclusions exists for major explosions (Métrich et al. 2001;2010;. 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 2 O, and CO 2 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.

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 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 2 O 3 and K 2 O/Na 2 O) are not affected by PEC, which affects Fe-Mg-Si in olivinehosted melt inclusions Rose-Koga et al. 2021). Similarly, H 2 O, CO 2 , 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 stepsize modelling following Girona and Costa (2013). Modelling 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 2 O and CO 2 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 2 O and CO 2 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 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 −1 , relative to a linear baseline. Carbonate peaks are measured at 1520 cm −1 and 1430 cm −1 , 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 3 , based on the value determined for basaltic glass fragments with 2.5 to 3.5 wt% H 2 O (Métrich et al. 2001). The detection limit for the CO 2 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 −1 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 2 O and CO 2 concentrations. Based on the propagation of the thickness measurement deviation, we estimate a maximum error of 12% for H 2 O and CO 2 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(Métrich et al. , 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.

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 -4 m/s and a seismic amplitude of 1 × 10 -4 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 2 /SO 2 ratio (median at 20 to 30 molar ratio) and of the volcanic CO 2 flux (median at ~ 2000 tons/day) (data from Aiuppa et al. 2021). This burst of deeply sourced (Aiuppa et al. 2021) CO 2 -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 66 and Fo 87 , predominantly showing values of Fo 80-87 (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 85measured either in homogeneous crystals or in cores of zoned crystals-to a maximum of Fo 87 . 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 82-84 that occur both as cores and rims (Fig. 4b).
Where they form the core, they show systematically reverse zoning patterns and occur surrounded by a group 1 rim (Fo 84-87 ; Fig. 4f). In the third group, compositions range from Fo 79-82 and exist exclusively as crystal cores, whereas a fourth group with Fo 76-79 is found exclusively as rim compositions. Finally, the fifth population (Fo 66-76 ) is subordinate and corresponds to inherited cores.
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(Pichavant et al. , 2022Métrich et al. 2021). The least evolved compositions (Fo 85-87 ) 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 >87 , as is common in large paroxysms (Métrich et al. , 2021, are identified.

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., 1998ME, 2003SSP, 2007SSP, and 1456Bertagnini et al. 2003;Métrich et al. 2005Métrich et al. , 2010, having CaO/Al 2 O 3 = 0.48-0.72 ± 4% and K 2 O/ Na 2 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 2 O/Na 2 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 and24 Nov 2009, La Felice andLandi 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).  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 2 /SO 2 molar ratios in the plume and b CO 2 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 2 -rich gas release is evident since late June 2020 Glassy melt inclusions display equivalent or less evolved compositions than matrix glass, with ratios of CaO/ Al 2 O 3 = 0.6-0.96 and K 2 O/Na 2 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 2 O 3 = 0.57-0.68, K 2 O/ Na 2 O = 0.61-0.86 (Fig. 6b). Melt inclusion CaO/Al 2 O 3 ratios are variable, possibly due to minor clinopyroxene crystallisation in the melt that leads to decreasing CaO/ Al 2 O 3 ratios, while K 2 O/Na 2 O ratios remain relatively constant (e.g., La Felice and Landi 2011). Host olivine compositions remain constant at ⁓Fo 85 , except for those MIs that are associated with host compositions Fo 79-84 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 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-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) Page 11 of 21 34 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, and 1456paroxysms (Bertagnini et al. 2003Métrich et al. 2005Métrich et al. , 2010Métrich et al. , 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 70-84 ) 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 2 O contents (Fig. 6c) (Fig. 6c).
Dissolved H 2 O and CO 2 contents show no systematic relationship with melt inclusion size and span a wide range with H 2 O = 1.53-3.64 wt% and CO 2 = 580-1866 ppm, with primary mode compositions of H 2 O = 2.8 wt% and CO 2 = 1500 ppm (Fig. 7a). Glass embayments show lower volatile contents with a primary mode at H 2 O = 2.4 wt% and CO 2 = 750 ppm and display a more limited compositional range (H 2 O = 1.72-2.4 wt%; CO 2 = 450-1639 ppm, Fig. 7a 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  Bertagnini et al. (2003) and Métrich et al. (2021). a K 2 O/Na 2 O and CaO/Al 2 O 3 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 2 O/Na 2 O = 11% and CaO/Al 2 O 3 = 4%. b Ratios of S/Cl and CaO/Al 2 O 3 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 2 O 3 = 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 crustal density of 2700 kg/m 3 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.

Source depth of the 19 July 2020 eruption
The comprehensive presented geochemical characterisation of the tephra ejected during the 19 July 2020 eruption includes olivine, MI, and matrix glass analyses. Olivine compositions and zoning patterns show large variability; however, the prevailing mode at Fo 85-87 , 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. , 2021. Matrix glasses of the 19 July 2020 eruptive products also cover a large 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 2 O 3 ). Error bar is Fo = 2% and CaO/Al 2 O 3 = 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 , 2003 and the 1998 major explosion (ME, dotted outline, Métrich et al. 2005Métrich et al. , 2010Métrich et al. , 2021 (Fig. 5). MIs and embayments record similarly variable compositions, which generally compare well with those observed in past major explosions (e.g., 1998) and smallscale paroxysms (2003,2007), but also suggest that eruption source and ascent dynamics may reflect a more complex history (see below).
Further constraints to allow better discussion of the ponding depth of major explosions are gained from total volatile (CO 2 + H 2 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 supplementary material for a more detailed discussion). Previously 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 2 O and CO 2 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 3 , 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  Métrich et al. (2005Métrich et al. ( , 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 2 contents that are below the detection limit and are reported to have CO 2 = 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 2 O and CO 2 contents for different scenarios with a starting composition of H 2 O = 2.9 wt% and CO 2 = 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 2 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 2 O and CO 2 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 3 ▸ 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 80-84 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 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 84-87 , pointing towards the second scenario.

Insights into the ascent dynamics of the 19 July 2020 eruption
Under open-system conditions, CO 2 is early exsolved from the melt while H 2 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 2 , decreasing the CO 2 melt composition, while H 2 O exsolution becomes increasingly dominant at lower pressures. Sulphur exsolution tends to be coupled with H 2 O exsolution and occurs generally at pressures < 150 MPa (Métrich et al. , 2021Lesne et al. 2011a,b;Aiuppa et al. 2021;Pichavant et al. 2013;, 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 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 and highlights a mode at < 5 days. The analysed paroxysms include Jul and Aug 2019 SSP, and the 1456 CE, the undated 16 th 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 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 2 -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 2 O, MIs show a trend of relatively constant S, at markedly decreasing H 2 O levels (Fig. 6c). A similar trend of near constant, but slightly increasing CO 2 contents at decreasing water contents can be observed in several MIs (trend "CO 2 buffering" in inset Fig. 7b), suggesting that the analysed 2020 melt inclusions testify partial H 2 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 2 -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 . 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 2 O and CaO/Al 2 O 3 ratios or with K 2 O and embayments showing lower sulphur contents (Fig. 6c). The high CO 2 values in 2020 melt inclusions (Fig. 7a, b) concur with previously described high CO 2 values of 1689 ppm, observed in melt inclusions containing carbonate crystals within bubbles (Métrich et al. 2001). Also, embayment compositions with high CO 2 , in combination with low CaO/Al 2 O 3 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 2 -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.  (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, and d represent a simplified conceptual model of Stromboli's plumbing system, where different source regions are activated for b lowintensity, c high-intensity major explosions, and d paroxysms and are highlighted in colour Role of the CO 2 -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 2 O and CO 2 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 2 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. 2010Aiuppa et al. , 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 largescale paroxysms (Andronico et al. 2008;Métrich et al. 2010Métrich et al. , 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. 20087 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 2 flux variations (compare the CO 2 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 and24 November 2009, La Felice andLandi 2011;Pioli et al. 2014) and are accompanied by a larger change in the CO 2 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 85-87 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 2 /SO 2 ratio in the degassing plume ( Fig. 3; Aiuppa et al. 2021). We ascribe this major system perturbation to the arrival of a CO 2 -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 2 /SO 2 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 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 89-90 ) 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 p and V s 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 89-90 ).
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 2 /SO 2 ratio and the CO 2 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 2 O + CO 2 ) 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 84-87 . 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 deepseated 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 2 -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 2 / SO 2 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 2 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 2 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.