Magma recharge and mush rejuvenation drive paroxysmal activity at Stromboli volcano (original) (raw)

Introduction

Basaltic volcanoes can remain active for tens to thousands of years with the continual presence of magma1, requiring storage and transport conditions that can sustain persistently eruptible melt. Magma storage conditions beneath these volcanoes may significantly change with time, leading to sudden and dramatic changes in explosivity2,3. Determining the rates and causes of these changes and how they modulate eruptive style over societally relevant timescales is of paramount importance for evaluating potential hazards. The three large-scale explosions (one major explosion and two paroxysms) that occurred at Stromboli volcano (Italy) over an uncommonly short time interval between 25 June and 28 August 2019 offer a unique opportunity to study the short-term variations in a basaltic plumbing system that can lead to paroxysmal events.

Stromboli is an active open conduit basaltic volcano (Fig. 1a) well-known for its persistent mild (normal) Strombolian activity occasionally interrupted by sudden, short-lived events ranging in size and intensity from major (violent Strombolian) to paroxysmal explosions4,5. Four paroxysms occurred at Stromboli over the last 70 years of recorded eruptive history prior to 2019. Each was preceded by week- to month-long periods of effusive activity6,7. Notably, the 2003 and 2007 paroxysms followed eruption of a similar cumulative volume of lava7, suggesting that lava volume may herald future paroxysms7 or that slow decompression rates may control the transition from effusive to violent (paroxysmal) explosions8. The prolonged effusive phase in May–October 2014, however, did not culminate in a violent explosion8,9.

Fig. 1: Map of Stromboli volcano and mineralogical map of 2019 tephra.

a Map of Stromboli volcano with its location in respect to Southern Italy shown in the inset. Mineralogical maps with colour-coded mineral phases of 2019 hp (b) and lp (c) tephra from the 3 July (samples P30-09 hp and P01-5 lp) paroxysmal events. Typical hp black scoriae with high-porphyritic seriate textures (~45–50 vol%), low vesicularity (~10–20 vol.%) and microlite-rich glassy groundmass have erupted during normal Strombolian activity, major explosions, and paroxysmal events. Light lp pumices characterised by weakly porphyritic seriate textures (<10 vol.% phenocrysts and microphenocrysts), high vesicularity (~30–40 vol.%) and a glassy groundmass are erupted during paroxysmal events and less frequently during major explosions. Phenocrysts and microphenocrysts of plagioclase (~20–30 vol.% of the total crystal content), clinopyroxene (~10–15 vol.%) and minor olivine (≤5 vol.%) are present in both hp and lp tephra. However, the larger plagioclase phenocrysts visible in the lp pumice in c (cyan crystals), are derived by exchange from _hp-_magma.

Full size image

The two paroxysms of 3 July and 28 August 2019 occurred without appreciable lava effusion, but Strombolian activity intensified abruptly on 9 June 201910, prior to the major explosion on 25 June. Frequent and unusual sustained lava fountaining persisted until 2 July11. The first paroxysm on 3 July was preceded by a small lava flow and edifice inflation ~45 and 10 minutes before the event11,12,13, respectively. The explosion formed an 8.4 km high eruption column11. Two pyroclastic currents descended the NW flank of the volcano (Sciara del Fuoco) and caused a ~1 m high tsunami10,11,14. After the 3 July paroxysm, explosive activity increased significantly for several weeks10 and discontinuous discharge of lava overtopped the southern crater area10,11,14,[15](/articles/s41467-022-35405-z#ref-CR15 "Stromboli Daily Bulletin. http://lgs.geo.unifi.it/index.php/reports/stromboli-daily

            ."). A second paroxysm on 28 August produced a 6.4 km high eruption column, a pyroclastic current, a smaller tsunami[10](/articles/s41467-022-35405-z#ref-CR10 "Andronico, D. et al. Uncovering the eruptive patterns of the 2019 double paroxysm eruption crisis of Stromboli volcano. Nat. Commun. 12, 4213 (2021)."),[11](/articles/s41467-022-35405-z#ref-CR11 "Giordano, G. & De Astis, G. The summer 2019 basaltic Vulcanian eruptions (paroxysms) of Stromboli. Bull. Volcanol. 83, 1 (2021)."), inflation of the edifice some 7 minutes before the event[16](/articles/s41467-022-35405-z#ref-CR16 "Giudicepietro, F. et al. Geophysical precursors of the July-August 2019 paroxysmal eruptive phase and their implications for Stromboli volcano (Italy) monitoring. Sci. Rep. 10, 1–16 (2020).") and a lava flow lasting a few hours. On 29 August, two major explosions occurred followed by lava effusion that lasted \~24 hours[14](/articles/s41467-022-35405-z#ref-CR14 "Plank, S. et al. The July/August 2019 Lava Flows at the Sciara del Fuoco, stromboli-analysis from multi-sensor infrared satellite imagery. Remote Sens. 11, 2879 (2019)."). Volcanic activity returned to a normal level by 20 September 2019[10](/articles/s41467-022-35405-z#ref-CR10 "Andronico, D. et al. Uncovering the eruptive patterns of the 2019 double paroxysm eruption crisis of Stromboli volcano. Nat. Commun. 12, 4213 (2021).").

Previous work on Stromboli’s 2019 paroxysms combining petrological and geophysical data17, implicated shallow conduit processes as triggers, while changing CO2 flux18 was interpreted to reflect an increasing amount of volatiles exsolved in the plumbing system since 2006, with an escalation in CO2 and SO2 degassing in the weeks to months before the paroxysms19. The focus of these earlier studies on the shallow plumbing system and concentration on either eruption precursors10,12, or possible triggering mechanisms17,18,19 have left key aspects of the role of the deeper plumbing system largely unexplored.

Here, we compare textural information and the compositional variability of clinopyroxene, glassy groundmass and bulk rocks from the major explosion on 25 June and the paroxysms on 3 July and 28 August with the products of normal Strombolian activity of 11 May 2019 and similar data of previous activity, particularly for 2003–20173. By decoding complex zonation in clinopyroxene, and integrating thermobarometry and diffusion modelling results, we reconstruct the temporal paths of clinopyroxene populations to illuminate the role of magma recharge and mush remobilisation in modulating eruption dynamics. Our approach allows the interpretation of magmatic processes and their surficial expression, as recorded by monitoring data, showing that magma-mush dynamics exercise a key control on eruption style and magnitude at persistently active basaltic volcanoes.

Results and discussion

Petrography and geochemistry of the 2019 major explosion and paroxysms

The persistent mild Strombolian activity, effusive eruptions and major explosions, all involve a degassed, highly porphyritic (hp) magma from a shallow reservoir, seismically identified at 1–3 km depth20. Deep-seated (10–12 km depth21) more mafic and, volatile-rich low-porphyritic (lp) magma is erupted, alongside _hp-_magma, during paroxysms21,22, and in smaller quantities during some of the major explosions23,24.

Two products were erupted during both the 3 July and 28 August paroxysms: (1) typical hp black scoriae with high-porphyritic seriate texture, low vesicularity and microlite-rich glassy groundmass. This lithotype also characterises normal Strombolian activity and was the only one erupted during the major explosion on 25 June (Fig. 1b); (2) light lp pumice characterised by a weakly porphyritic seriate texture, high vesicularity and a glassy groundmass (Fig. 1c). Both endmembers have paragenesis of plagioclase, clinopyroxene and minor olivine (Fig. 1b, c), typical of products of the Present-day activity of Stromboli3,22,25,26 (i.e., <1.2 ka27). The two endmembers were syn-eruptively mingled, resulting in some homogenisation and crystal exchange between the two magmas (e.g., plagioclase in lp ejecta was largely derived by exchange from _hp-_magma, see Fig. 1c).

Bulk compositions of summer 2019 tephra are shoshonitic basalts typical of the Present-day activity and chemically similar to the _lp-_magma groundmass glass (Supplementary Fig. 1 and Supplementary Tables 1, 2). However, the glassy groundmass is compositionally bimodal (Supplementary Fig. 1).

Clinopyroxene textures and compositions

Clinopyroxene (Supplementary Tab. 3) displays a broad compositional range from diopside (Di69–80Hd9–16En8–15Fs1–3, Mg#80–90) to augite (Di50-67Hd18-25En10-20Fs2-9, Mg#69-79) (see Methods) (Fig. 2b), as previously documented for the Present-day activity and interpreted as caused by the continuous lp-hp magma mixing and antecryst recycling from crystal mush3,28,29,30. Augitic compositions are in equilibrium with the more evolved resident _hp_-magma. Diopsidic compositions are in equilibrium with the more mafic _lp_-magma and their presence record inheritance from a mush region or recharging events of the more primitive _lp_-magma. Clinopyroxene can be texturally and chemically categorised into two main groups and several subgroups (Fig. 2).

Fig. 2: Textural classification, chemical composition, and frequency of 2019 clinopyroxenes.

Clinopyroxenes from the major explosion of 25 June and the two paroxysms of 3 July and 28 August are euhedral with well-defined edges and tabular forms, with both phenocrysts (Lmax > 0.3 mm, on the basis of the longest size dimension) and microphenocrysts (Lmax = 0.3–0.1 mm; cf. Mollo et al.67) present (a) BSE (back-scattered electron) images for the two groups and respective subgroups of clinopyroxenes populations identified in the 25 June major explosions and the 3 July and 28 August paroxysms. The white numbers indicate the Mg# (Mg/(Mg+Fe2+tot) at.) for the different portions of clinopyroxene and correspond to augitic (Mg#69–79) and diopsidic (Mg#80–90) compositions. Subgroup 1A (zoned augite): Px3 sample P30-09-2-hp, 3 July paroxysm; subgroup 1B (diopsidic overgrowth): Px4-5 sample P24-3-lp, 3 July paroxysm; subgroup 1C (diopsidic band): Px35 sample 25062019, 25 June major explosion; subgroup 2A (zoned diopside): Px2 sample P01-05-lp, 3 July paroxysm and Px 7 (homogenous diopside), sample P44-1-lp, 28 August paroxysm; subgroup 2B (augitic overgrowth): Px1 sample P21-1-lp, 3 July paroxysm; b Mg#-overgrowth vs Mg#-core diagram for the 2019 clinopyroxenes. The chemical composition of each clinopyroxene core is plotted against the corresponding composition of the overgrowth for 25 June (black), 3 July hp (red) and lp (orange), 28 August hp (blue) and lp (cyan). The two main groups with the relative subgroups are indicated. Only the subgroup 1C is not plotted due to the low abundance of this crystal population. Mg# error bars are reported in the lower right quadrant. c Abundance (%) of clinopyroxene in each subgroup according to eruption (major, paroxysm) and magma (hp, lp, mingled) type. The total for all 2019 clinopyroxene is also given. N: number of clinopyroxenes for each subgroup.

Full size image

Group 1 consists of clinopyroxene phenocrysts and microphenocrysts with an augitic core. Crystals are nearly homogeneous or slightly zoned (i.e., zoned augite subgroup 1A) or, alternatively, crystals have augitic cores with diopsidic rims (i.e., diopsidic overgrowth; subgroup 1B) reaching hundreds of microns (Fig. 2). Rare crystals from Group 1 show diopsidic recharge bands, appearing as sharp to diffused boundaries around euhedral to resorbed augitic cores and euhedral augitic rims (subgroup 1C). Group 1, particularly 1A is modally dominant (85% of crystals; Fig. 2c).

Group 2 clinopyroxene phenocrysts and microphenocrysts (15% of crystals; Fig. 2c) have diopsidic cores (Fig. 2a, b). These crystals are homogeneous or slightly zoned (subgroup 2A zoned diopside), or they have diopsidic cores with dissolution-reaction features (i.e., cores are resorbed and likely antecrystic) and augitic overgrowths (subgroup 2B; Fig. 2). It is worth noting that not all subgroups are represented in each tephra. The main differences in clinopyroxene cargoes are summarised in Fig. 2c.

Chemical mapping of augite crystals of subgroup 1A and augitic overgrowths of subgroup 2B crystals (Fig. 3, Supplementary Figs. 2–3) shows a finely banded, oscillatory zoning pattern in Cr, Al, and Ti. In subgroup 1A crystals, this oscillatory zoning is superimposed on the main low-amplitude augitic core-rim zoning where ΔMg# (=Mg#overgrowth – Mg#core) ranges between + 5 and −5. The increase in Mg# in the zoned augite crystals of subgroup 1A may or may not correspond to a modest increase in Cr content (Fig. 3, Supplementary Figs. 2–3). By contrast, augitic cores of subgroup 1B crystals are homogenous (Fig. 3, Supplementary Fig. 2). Group 2 and diopsidic rims of subgroup 1B crystals typically have elevated concentrations of Mg, Al, Ni and Cr (>2500 ppm) (Fig. 3, Supplementary Figs. 2–3).

Fig. 3: Major and trace elements chemical mapping for 2019 clinopyroxenes.

BSE images (top) and quantitative chemical mapping of Mg, Al, Ti (wt%), Cr and Ni (ppm) for the main four subgroups of clinopyroxenes. Subgroup 1A: Px 11 sample P27-1-hp, 3 July paroxysm; subgroup 1B: Px6 sample P30-07-lp, 3 July paroxysm; subgroup 2B: Px1 sample P-21-1-lp, 3 July paroxysm; subgroup 2A: Px13 sample P30-07-lp, 3 July paroxysm. The scale bar for Cr is on the left for subgroup 1A, 1B and 2B and on the right for subgroup 2A.

Full size image

2003-2019 clinopyroxene geochemical evolution and _P_-_T_-H2O crystallisation conditions

Composition of diopsidic and augitic clinopyroxenes from the major explosion of 25 June and the paroxysms of 3 July and 28 August 2019 notably differ from those erupted during the 2003–2017 activity3 (Fig. 4a, b). In particular, crystal populations from 2019 eruptions are more enriched in TAl + M1Ti + M1Fe3+ and depleted in TSi + M1Mg + M1Fe2+ compared with those from 2003–2017 eruptions (Fig. 4a) with Tschermak components preferentially incorporated in 2019 crystals at the expense of diopside (Di; Fig. 4b). Nevertheless, both 2019 and 2003–2017 clinopyroxenes align along a single trend of Tschermak component (CaTs + CaTiTs + CaFeTs) substitution (Fig. 4a), suggesting a common geochemical evolution, which is consistent with the relative magma uniformity across 2003-2019 and the entire Present-day activity (Supplementary Fig. 1).

Fig. 4: Chemical composition of 2019 clinopyroxenes and P-T crystallisation conditions.

a Al+Ti+Fe3+ (apfu) vs. Si+Mg+Fe2+ (apfu) and (b) CaTs+CaFeTs+CaTiTs (mol.%) vs. Diopside (Di, mol%) diagrams showing the chemical composition and range of variability of augitic and diopsidic portions of the 2019 clinopyroxenes (black triangle: 25 June major explosion; orange and red triangle: 3 July paroxysm; and cyan and blue triangle: 28 August paroxysm) compared with the augitic and diopsidic composition of 2003–2017 clinopyroxenes (grey field and square) from Di Stefano et al.3. At the scale of the figure, error bars are smaller than the symbols. CaTs+CaFeTs+CaTiTs (mol.%) show the Tschermak-style cation substitution [TSi, M1Mg, M1Fe2+] → [TAl, M1Ti, M1Fe3+]. Probability density function of temperature (c) and pressure (d) for 2019 augite (black line: 25 June, based on 28 datapoints; red line: 3 July, 82 datapoints; and blue line: 28 august, 56 datapoints) and diopside (orange line: 3 July, 37 datapoints; cyan: 28 August, 44 datapoints) compositions compared with the P-T range of 2003-2007 eruptive period (grey band) from Di Stefano et al.3. Clinopyroxene-based pressure estimates do not vary significantly among the different diopsidic-augitic crystal portions due to the lack of substantial changes in the _P_-sensitive jadeite-melt exchange reaction at relative low-pressure conditions68,69. Estimated H2O contents for 2019 (this study) and 2003–20173 are also given. Augite equilibrates at a slightly higher melt-H2O content relative to diopside due to the higher portion of solid fraction during the late stage of clinopyroxene growth, before abundant _hp_-magma degassing and plagioclase crystallisation during magma ascent towards the surface3. Indeed, the amount of H2O dissolved in Stromboli magma reaches ~2.7–3.5 wt.% prior to outgassing21, in agreement with the maximum melt-H2O content recorded by 2019 clinopyroxene. There is also a positive correlation between the ratio of iron to magnesium in clinopyroxene and the amount of H2O dissolved in the melt67. In particular, the M1Fe2+/M1Mg ratio of 2019 clinopyroxene is systematically higher than that of 2003–2017 crystals, which is consistent with their higher melt-H2O concentrations.

Full size image

Enhanced Tschermak substitution may reflect higher temperature crystallisation31,32 or more hydrous conditions33. Probability density functions3 for T show that 2019 diopsidic and augitic compositions peak at ~1175 ± 16 and ~1155 ± 16 °C, respectively (Fig. 4c, see Methods for T-P-H 2 O calculation schemes). Comparing 2019 diopsidic clinopyroxene with those erupted from the Present-day activity3,34 indicates preferential equilibration of diopside within the deep-seated _lp_-magma at T > 1170 °C, with near-liquidus mineral phase appearance at ~1200 °C. Conversely, 2019 augitic clinopyroxene saturates the _hp_-magma at temperatures higher than those estimated for the Present-day activity (i.e., T < 1145 °C3,34), indicating a hotter reservoir in 2019.

Probability density functions for P show that both 2019 diopsidic and augitic compositions peak at ~150 ± 100 MPa (Fig. 4d), consistent with a near-invariant pressure range obtained for the 2003–2017 crystals (~180 ± 100 MPa3). The stability of the magmatic reservoir, with magma ponding at the same level (also observed in olivine5), confirms that Stromboli’s plumbing system is characterised by a shallow reservoir located at ~50–150 MPa (i.e., ~1.5–6 km depth assuming a crustal density of 2700 kg/m3;5,21,25). In this reservoir, the resident _hp_-magma crystallises augitic clinopyroxene, olivine and abundant plagioclase22,29,30,35,36. An additional deeper reservoir, inferred from melt inclusion and water solubility data21, is located at ~190–260 MPa (~7–10 km depth21), where the feeding _lp_-magma undergoes limited crystallisation of diopsidic clinopyroxene and olivine3,21,30. Thus, the steady-state dynamics established since at least ~1.7 ka ago, and sustained by continuous refilling of a shallow _hp-_reservoir by volatile-rich _lp_-magma30, continue to control the system’s evolution. Syn-eruptive mingling of compositionally similar lp-hp basaltic magmas leads to the observed range in crystallinity3,22,28,29,30,36.

Melt-H2O concentrations in equilibrium with 2019 clinopyroxene range from 0.5–3.1 wt.% (diopside) and 0.5–3.6 wt.% (augite), slightly higher than those of the 2003–2017 period (0.5–2.4 and 0.5–2.9 wt.% for diopside and augite, respectively) and consistent with the Fe:Mg ratio in clinopyroxene (Fig. 4).

Mafic recharge as eruption trigger of the 2019 paroxysms

Diopsidic portions of subgroups 1B, 1C, 2A and 2B exhibit high Cr concentrations, elevated concentrations of Ni, Mg and Al and lower concentrations of Ti (Fig. 3, Supplementary Figs. 2–3). The elevated Cr reflects crystallisation from mafic, Cr-rich _lp_-magma. Diopsidic cores of Group 2 crystals likely equilibrated at depth, whereas diopsidic rims (and diopsidic bands) on augitic clinopyroxene formed because of mafic recharge due to convective mass transfer37 from depth of _lp-_magma into _hp-_magma. Thus, subgroup 2A crystals reflect crystallisation entirely within the deep lp magmatic domain that is compositionally distinct from the shallow hp reservoir, where augitic clinopyroxene crystallises3,22,29,30,36. By contrast, resorbed diopsidic and augitic cores of subgroup 2B and some 2C crystals, respectively, are interpreted as disrupted portions of a vertically-extended mush column situated at the base of the _hp_- reservoir28,36,38 or encompassing both the shallow and deeper reservoirs3,39 and, thus, as antecrysts in line with the previous studies3,28,30,36,38,39. The predominance of subgroup 1B crystals in lp products and their abundance in hp and mingled products erupted during the paroxysms suggest interaction between hp and lp magmas prior to eruption. This further suggest that new injected mafic _lp-_magma(s) induced mixing and crystallisation (i.e., the diopsidic overgrowth of subgroup 1B) and triggered the 2019 paroxysms. At the same time, the absence of subgroup 1B and group 2 (particularly subgroup 2B) crystals, alongside the scarcity of subgroup 1C in the major explosion of 25 June and the presence of subgroups 1B and group 2 exclusively in paroxysmal products of 3 July and 28 August (Fig. 2c) clearly indicates mush remobilisation only occurred during these powerful events.

Fe-Mg diffusive timescales of recharge and magma mixing events were determined on 65 clinopyroxenes from the major explosion of 25 June and the two paroxysms of 3 July and 28 August 2019 through the non-isothermal approach (NIDIS)40 (see Methods). Timescales for all three eruptive events range from a few days to 6 years (Table 1, Fig. 5 and Supplementary Figs. 6–64), with 74% of crystals recording timescales \(\le{1}\) month pre-eruption (Fig. 5a, b). Clinopyroxenes from lp tephra erupted during the two paroxysms return short timescales from a few days to a couple of months (Table 1), with mingled and hp products from the two paroxysms indicating timescales in the order of days to ~4 years. The major explosion on 25 June shows longer timescales between 6 ± 2 months and 6 ± 2 years, with only one crystal recording a timescale of a few days (Fig. 5b, Table 1). In stark contrast are the longer timescales (~1–182 years) reported for the 2003–2017 period3, although these timescales were determined for clinopyroxenes characterised by diopsidic resorbed cores and diopsidic bands (similar to our subgroup 2B and 1C). In fact, according to Di Stefano et al.3 the 2003–2017 period is characterised by a lack of subgroup 1B crystals and a predominance of diopsidic resorbed cores, but clinopyroxenes with diopsidic rims occur in both 2003 _lp_-pumices and 2007 _hp_-lava flows39,41,[42](/articles/s41467-022-35405-z#ref-CR42 "Francalanci, L. et al. Mineralogical, geochemical, and isotopic characteristics of the ejecta from the 5 april 2003 paroxysm at stromboli, italy: inferences on the preeruptive magma dynamics. In Geophysical Monograph Series, 331–345 (2008). https://doi.org/10.1029/182GM27

            ."), although no timescales information is available for these crystals.

Fig. 5: Fe-Mg timescales for the 2019 clinopyroxenes.

a Histogram showing the number of clinopyroxene vs timescales (\(\triangle\)time injection – eruption, in days and years) from mafic injection to eruption for 25 June (black) major explosion, 3 July (red) and 28 August (blue) paroxysms. Timescales for the 2003–2017 eruptive period from Di Stefano et al.3 are also shown (grey), alongside the timescales calculated from olivine crystals by Métrich et al.5 for the two 2019 paroxysms (green). n = number of timescales estimate. b Timescales are presented in order of increasing \(\triangle\)time and are plotted according to the different type of erupted magmas (filled circle) – black: 25 June major explosion _hp_-magma; orange: 3 July and 28 August paroxysms _lp_-magma; purple: 3 July and 28 August paroxysms _hp_-magma; green: 3 July and 28 August paroxysms _mingled_-magma. Timescales for the 2003–2017 eruptive period from Di Stefano et al.3 are also shown (grey triangle). c Timescales (in days) for the 2019 major explosion and paroxysms are presented according to different pyroxenes subgroups. The majority of timescales obtained from the 3 July and 28 August (lp and mingled composition) paroxysms were calculated from subgroup 1B clinopyroxenes, which therefore represent the time between growth of the diopsidic rim and eruption caused by injection of mafic magma into the shallow reservoir prior to eruption. These short timescales correspond to the eruption-triggering event as discussed in the text.

Full size image

Table 1 Fe-Mg timescales calculated for clinopyroxenes from 2019 eruptions

Full size table

The very short timescales from the diopsidic overgrowth of subgroup 1B clinopyroxenes (a few days- to a couple of months, with two outliers of 128 and 224 days, Fig. 5c) indicate the magmatic system responded rapidly to injection of mafic magma(s) prior to the two paroxysms. In contrast to other proposed models5,17,19, interaction between the lp- and _hp-_magma in the Stromboli’s shallow reservoir extends beyond syn-eruptive mingling with lp mafic recharge arriving in the shallow reservoir months before the 2019 paroxysmal events inducing limited and incomplete magma mixing associated with diopsidic overgrowth induced crystallisation. Accordingly, 1B crystals confirm rapid overgrowth of diopsidic portions around augitic cores following injection of _lp_-magma into the shallow _hp_-reservoir (Fig. 3, Supplementary Figs. 2–3). The importance of lp mafic recharge actively interacting with the _hp-_magma in the months prior the 2019 paroxysms is emphasised by the prevalence of subgroup 1B crystals in lp products and their abundance in hp and mingled samples of the paroxysms (Fig. 2). Moreover, the abundance in the lp samples of zoned diopsidic 2A clinopyroxene, which formed in the deeper more mafic lp recharge magma and lack augitic overgrowths (Figs. 2–3, Supplementary Figs. 2–3 and 59–62), suggests incomplete mixing and limited residence time of the recharge magma(s) in the shallow reservoir. Short crystal residence times (5–110 days for 3 July and 10-115 days for 28 August, albeit without uncertainties) were also derived from olivine by Métrich et al.5 (Fig. 5a). Our recorded recharges from month(s) to days before the 3 July and 28 August 2019 events (Fig. 5b, c) substantiate the role of the deeper magmatic system in triggering both paroxysms, well beyond volatile flushing13,17,18,19 by feeding fresh _lp_-magma continuously into the shallow _hp_-reservoir from ~1–2 months prior to each paroxysmal event.

Regarding the major explosion of 25 June, the absence of subgroup 1B clinopyroxenes, and the longer timescales (days to 6 years) measured in subgroup 1A and 1C crystals (see below) demonstrates that _lp_-magma did not directly triggered this eruption.

Convective mixing regime controlling magma dynamics

Both diopsidic and augitic clinopyroxenes from the major explosions and paroxysms in 2019 overlap compositionally with crystals from 2003–2017 (Fig. 6, Supplementary Tab. 4). However, notable differences exist in the uptake of elements in the lattice sites (Fig. 6, Supplementary Fig. 4). Dynamic stirring experiments37 indicate that fine low-amplitude concentric zoning may result from kinetic crystallisation under convective stirring. Therefore, Cr-rich concentric zones may not necessarily correlate with other geochemical indices of recharge37. Augitic portions and overgrowth of subgroup 1A, 1C, and 2B crystals, and some diopsidic portions of subgroup 2A crystals are characterised by fine oscillatory zoning, sometimes associated with Cr-rich concentric zones (Figs. 2 and 3, Supplementary Figs. 2–3) with limited correlation with other elements. This contrasts with the incorporation of Cr in clinopyroxene resulting from mafic recharge. This fine oscillatory zoning in the augitic portions of 1A, 1C and 2B (and some 2A) crystals may result from convective stirring, hinting that during the 2019 major and paroxysmal activity, the shallow _hp_-reservoir, and possibly also the deeper _lp_-reservoir, underwent continuous convection, promoting crystallisation and minimising kinetic effects on clinopyroxene zoning37.

Fig. 6: Trace elements of 2019 clinopyroxenes.

a Yb; b Sc; c Zr, and d Cr vs La (ppm) diagrams showing the trace elements range of variability of augitic (Aug) and diopsidic (Di) portions of the 2019 clinopyroxenes (black triangle: 25 June major explosion; orange and red triangle: 3 July paroxysm; and cyan and blue triangle: 28 August paroxysm) compared with the augitic and diopsidic composition of 2003–2017 clinopyroxenes (grey field) from Di Stefano et al.3. At the scale of the figure, error bars are smaller than the symbols.

Full size image

Superimposed on the fine oscillatory zoning, zoned augites and diopsides (subgroups 1A and 2A) show slightly larger but still relatively low-amplitude compositional variations (Figs. 2 and 3). Several subgroup 1A crystals and some subgroups 2A (e.g., Fig. 2, Supplementary Figs. 2–3, 6–14, and 60–61) crystals have rounded cores, indicating dissolution and reaction with the crystallising magma (hp or lp respectively) but with some minor compositional variation resulting from local stronger convective regime whereby growth kinetics and diffusion boundary layers were not reduced by the homogenising effect of stirring37. Modification of these compositional boundaries by diffusion, where the slow diffusing elements Al and Ca show a sharp change that mimics the Mg profile (Supplementary Figs. 6–14 and 59–62) can therefore be used to obtain time constraints of the convective regime.

Fe-Mg diffusion timescales of the larger compositional zoning for zoned augite of subgroups 1A range from 543 ± 228 to a few days, for the major explosions of 25 June and the two paroxysms, whereas those for zoned diopside of subgroup 2A crystals are on average shorter, ranging from days to <2 months (Table 1; Fig. 5c, Supplementary Figs. 6–14 and 59–62). These timescales overlap with those of the mafic triggering events (see above), suggesting magnification of convection by mafic recharge event(s). The presence of zoned augites, alongside diopsidic and augitic overgrowth, indicates that mafic recharge destabilised the shallower parts of the system, thereby establishing a chaotic mixing regime similar to that postulated for Post-Pizzo period30 (1.7–1.5 ka43), favouring the formation of subgroup 1A (i.e., crystallisation of the oscillatory zoned augitic overgrowth). Similarly, the presence of zoned diopside subgroup 2A crystals with shorter timescales suggest that during June-August 2019, dynamic crystallisation also characterised the deeper _lp_-reservoir, leading to diopsidic overgrowth of 2 A crystals around pre-existing diopsidic cores, shortly before mafic recharges rapidly triggered the paroxysms.

The lifecycle of basaltic crystal mush

Despite the steady-state and similar configuration of Stromboli’s plumbing system, the marked textural and compositional differences between the 2019 clinopyroxenes and those from the 2003–2017 period3,38,39,[42](/articles/s41467-022-35405-z#ref-CR42 "Francalanci, L. et al. Mineralogical, geochemical, and isotopic characteristics of the ejecta from the 5 april 2003 paroxysm at stromboli, italy: inferences on the preeruptive magma dynamics. In Geophysical Monograph Series, 331–345 (2008). https://doi.org/10.1029/182GM27

            .") (Figs. [2](/articles/s41467-022-35405-z#Fig2) and [4](/articles/s41467-022-35405-z#Fig4)), but also from previous activity (i.e., 1984-1996[28](/articles/s41467-022-35405-z#ref-CR28 "Francalanci, L. et al. Intra-grain Sr isotope evidence for crystal recycling and multiple magma reservoirs in the recent activity of Stromboli volcano, Southern Italy. J. Pet. 46, 1997–2021 (2005).")), require significantly different magma dynamics (Fig. [7](/articles/s41467-022-35405-z#Fig7)). A striking difference is the paucity of crystals with resorbed (antecrystic) diopsidic and augitic cores (subgroup 2B and some subgroup 1C crystals, Fig. [2](/articles/s41467-022-35405-z#Fig2), Supplementary Figs. [3](/articles/s41467-022-35405-z#MOESM1), [56](/articles/s41467-022-35405-z#MOESM1)–[57](/articles/s41467-022-35405-z#MOESM1) and [63](/articles/s41467-022-35405-z#MOESM1)–[64](/articles/s41467-022-35405-z#MOESM1)), which are rare in the tephra from the 3 July and 28 August paroxysms (Fig. [2c](/articles/s41467-022-35405-z#Fig2)) and absent (subgroup 2B) from the 25 June ejecta. Abundant resorbed diopsidic and augitic cores texturally and compositionally similar to our subgroup 2B and 1C crystals were reported in 2003–2017 products[3](/articles/s41467-022-35405-z#ref-CR3 "Di Stefano, F. et al. Mush cannibalism and disruption recorded by clinopyroxene phenocrysts at Stromboli volcano: new insights from recent 2003–2017 activity. Lithos 360–361, 105440 (2020)."),[28](/articles/s41467-022-35405-z#ref-CR28 "Francalanci, L. et al. Intra-grain Sr isotope evidence for crystal recycling and multiple magma reservoirs in the recent activity of Stromboli volcano, Southern Italy. J. Pet. 46, 1997–2021 (2005)."),[38](/articles/s41467-022-35405-z#ref-CR38 "Francalanci, L. et al. Crystal recycling in the steady-state system of the active Stromboli volcano: a 2.5-ka story inferred from in situ Sr-isotope and trace element data. Contrib. Mineral. Petrol. 163, 109–131 (2012)."),[39](/articles/s41467-022-35405-z#ref-CR39 "Landi, P. et al. Magma dynamics during the 2007 Stromboli eruption (Aeolian Islands, Italy): mineralogical, geochemical and isotopic data. J. Volcanol. Geotherm. Res. 182, 255–268 (2009)."),[41](/articles/s41467-022-35405-z#ref-CR41 "Métrich, N., Bertagnini, A., Landi, P., Rosi, M. & Belhadj, O. Triggering mechanism at the origin of paroxysms at Stromboli (Aeolian Archipelago, Italy): the 5 April 2003 eruption. Geophys. Res. Lett. 32, 1–4 (2005)."),[42](/articles/s41467-022-35405-z#ref-CR42 "Francalanci, L. et al. Mineralogical, geochemical, and isotopic characteristics of the ejecta from the 5 april 2003 paroxysm at stromboli, italy: inferences on the preeruptive magma dynamics. In Geophysical Monograph Series, 331–345 (2008). 
              https://doi.org/10.1029/182GM27
              
            ."), and in ejecta erupted since at least 1900[5](/articles/s41467-022-35405-z#ref-CR5 "Métrich, N., Bertagnini, A. & Pistolesi, M. Paroxysms at Stromboli Volcano (Italy): source, genesis and dynamics. Front. Earth Sci. 9, 1–17 (2021)."),[22](/articles/s41467-022-35405-z#ref-CR22 "Francalanci, L., Tommasini, S. & Conticelli, S. 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. Geotherm. Res. 131, 179–211 (2004)."),[26](/articles/s41467-022-35405-z#ref-CR26 "Pompilio, M., Bertagnini, A. & Métrich, N. Geochemical heterogeneities and dynamics of magmas within the plumbing system of a persistently active volcano: evidence from Stromboli. Bull. Volcanol. 74, 881–894 (2012)."),[28](/articles/s41467-022-35405-z#ref-CR28 "Francalanci, L. et al. Intra-grain Sr isotope evidence for crystal recycling and multiple magma reservoirs in the recent activity of Stromboli volcano, Southern Italy. J. Pet. 46, 1997–2021 (2005)."). They were interpreted as portion of a high 87Sr/86Sr crystal mush formed since 2.5 ka ago located at the base of the _hp_\-reservoir[36](/articles/s41467-022-35405-z#ref-CR36 "Bragagni, A., Avanzinelli, R., Freymuth, H. & Francalanci, L. Recycling of crystal mush-derived melts and short magma residence times revealed by U-series disequilibria at Stromboli volcano. Earth Planet. Sci. Lett. 404, 206–219 (2014)."),[38](/articles/s41467-022-35405-z#ref-CR38 "Francalanci, L. et al. Crystal recycling in the steady-state system of the active Stromboli volcano: a 2.5-ka story inferred from in situ Sr-isotope and trace element data. Contrib. Mineral. Petrol. 163, 109–131 (2012)."), (i.e., >3 km depth) or encompassing both reservoirs[21](/articles/s41467-022-35405-z#ref-CR21 "Metrich, N., Bertagnini, A. & Muro, A. D. 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 (2010)."),[44](/articles/s41467-022-35405-z#ref-CR44 "Aiuppa, A. et al. Unusually large magmatic CO2 gas emissions prior to a basaltic paroxysm. Geophys. Res. Lett. 37, 1–5 (2010)."). Clinopyroxene and plagioclase antecrysts are periodically remobilised from the bottom by the ascent _lp_\-magma throughout this period, and recycled from the top by sinking of the dense degassed _hp-_magma left in the conduit during the normal Strombolian activity[28](/articles/s41467-022-35405-z#ref-CR28 "Francalanci, L. et al. Intra-grain Sr isotope evidence for crystal recycling and multiple magma reservoirs in the recent activity of Stromboli volcano, Southern Italy. J. Pet. 46, 1997–2021 (2005)."),[30](/articles/s41467-022-35405-z#ref-CR30 "Petrone, C. M., Braschi, E., Francalanci, L., Casalini, M. & Tommasini, S. Rapid mixing and short storage timescale in the magma dynamics of a steady-state volcano. Earth Planet. Sci. Lett. 492, 206–221 (2018)."),[36](/articles/s41467-022-35405-z#ref-CR36 "Bragagni, A., Avanzinelli, R., Freymuth, H. & Francalanci, L. Recycling of crystal mush-derived melts and short magma residence times revealed by U-series disequilibria at Stromboli volcano. Earth Planet. Sci. Lett. 404, 206–219 (2014)."),[38](/articles/s41467-022-35405-z#ref-CR38 "Francalanci, L. et al. Crystal recycling in the steady-state system of the active Stromboli volcano: a 2.5-ka story inferred from in situ Sr-isotope and trace element data. Contrib. Mineral. Petrol. 163, 109–131 (2012)."). Accordingly, we interpret the resorbed diopsidic and augitic cores of subgroups 2B and 1C, respectively, as antecrysts remobilised from this long-lived[28](/articles/s41467-022-35405-z#ref-CR28 "Francalanci, L. et al. Intra-grain Sr isotope evidence for crystal recycling and multiple magma reservoirs in the recent activity of Stromboli volcano, Southern Italy. J. Pet. 46, 1997–2021 (2005)."),[38](/articles/s41467-022-35405-z#ref-CR38 "Francalanci, L. et al. Crystal recycling in the steady-state system of the active Stromboli volcano: a 2.5-ka story inferred from in situ Sr-isotope and trace element data. Contrib. Mineral. Petrol. 163, 109–131 (2012).") extant crystal mush.

Fig. 7: Schematic model of the plumbing system at Stromboli volcano.

a 2003–2017 eruptive period adapted from Di Stefano et al.;3 b June–August 2019 (this study). During the 2003–2017 eruptive period extensive mush remobilisation and cannibalisation operated by the uprising deeper undegassed mafic recharge _lp_-magma (in red) characterises magma dynamics. The _lp_-magma is fed by the deeper reservoir as inferred by Métrich et al.21. Abundant recycled diopsidic antecrysts are remobilised from the crystal mush and transported by the _lp_-magmas permeating the mush, into the shallower _hp_-reservoir (yellow-orange). The continuous and efficient remobilisation of the mush between 2003 and 2017 led to a rejuvenated plumbing system configuration in 2019, characterised by the _lp_-magma extensively permeating the mush with limited mush remobilisation. A few recycled diopsidic and augitic antecrysts are remobilised and transported by the _lp_-magma into the shallower _hp_-reservoir where the augitic overgrowth forms around resorbed diopsidic core (subgroup 2B) or resorbed augitic core surrounded by diopsidic band (subgroups 2C). Zoned or homogenous diopside (subgroup 2A) and augite (subgroup 1A) formed within the _lp_- and _hp_-reservoirs, respectively. The continuous arrival of _lp-_magma determines a convective mixing (black rounded arrows) in both reservoirs. Such highly dynamic environments originate both the zoned crystals and the low-amplitude zonation of the overgrowth portion. The abundance of diopsidic overgrowth (subgroup 1B) and their short timescales indicate that the 2019 paroxysmal events were triggered by injection of new _lp_-magma from a few months to a few days before the explosion onset. The schematic model has been modified and adapted from Di Stefano et al.3.

Full size image

In the 2019 products, about half of the handful of crystals of subgroup 1C (2% of the total, Fig. 2c, Supplementary Figs. 56–57) have resorbed augitic cores with diopsidic overgrowth and augitic rims. Resorbed diopsidic core (i.e., subgroup 2B) are slightly more abundant (4%), bringing at ~5% the total proportion of clinopyroxene with resorbed cores (Fig. 2c). The proportion of diopsidic clinopyroxene antecrysts was estimated at ~33% for the 2003–2017 products3, which also include normal and major activity, demonstrating that the proportion of clinopyroxene antecrysts in 2019 ejecta was significantly lower than in the 2003–2017 period. As our focus is on the clinopyroxene population, only a less direct comparison can be drawn from the proportion of clinopyroxene and plagioclase antecrysts in the products of the 5 April 2003 paroxysm estimated at ~17% with an increase of ~8% in a lava erupted after 5 April 2003 paroxysm38. However, considering that the total crystallinity of both _hp_-scoriae and _lp_-pumice has remained relatively constant since at least A.D. 776-135027,28 and that resorbed antecryst augitic and diopsidic cores are reported as a common feature of 1984-1996 ejecta28, our data suggest that the proportion of clinopyroxene antecrysts in the eruptive products of summer 2019 is also significantly lower than stretching back to at least 1984. Notably, resorbed augitic cores are not reported among the large clinopyroxene antecryst population examined for the 2003–2017 period3, although their presence is noted in the pumice of 5 April 2003 paroxysm and 2007 lavas39,41,[42](/articles/s41467-022-35405-z#ref-CR42 "Francalanci, L. et al. Mineralogical, geochemical, and isotopic characteristics of the ejecta from the 5 april 2003 paroxysm at stromboli, italy: inferences on the preeruptive magma dynamics. In Geophysical Monograph Series, 331–345 (2008). https://doi.org/10.1029/182GM27

            ."). Resorbed augitic cores are also less abundant than the diopsidic core in the scarce antecryst population in 2019 products suggesting a greater contribution from the deeper diopsidic-dominated portion of the crystal mush. Thus, at least since 2003 (unfortunately similar quantitative data are unavailable for older periods of activity) clinopyroxene has recorded a higher degree of remobilisation from a diopside-dominated mush which, we suggest, forms the lower portion of the mushy system (Fig. [7](/articles/s41467-022-35405-z#Fig7)).

The overall scarcity of crystals with resorbed cores indicates that the dynamics of Stromboli during June-August 2019 involved limited interaction of ascending magma with extant mush, contrasting with the 2003–2017 period of activity (Fig. 7) and, possibly also, for activity since at least 1984 if not for the entire 20th century, as testified by the abundance of antecrysts5,22,26,28. Efficient mush disruption and cannibalism has been proposed for the 2003–2017 period when a large proportions of diopsidic antecrysts were entrained by the _lp_-magma3 (Fig. 7a). We postulate that the continuous remobilisation of the mush led to a rejuvenated plumbing system in 2019, at least in the central conduit. The paucity of antecrysts in the 2019 paroxysms and their absence during major explosive activity supports the hypothesis of a mush rejuvenation process, with limited interaction between the ascending _lp-_recharge magma and the mush (Fig. 7b).

Furthermore, such a model is supported by timescales from subgroup 2B crystals and those with resorbed cores from subgroup 2C (Fig. 5c, Supplementary Figs. 56–57 and 63–64, Table 1), which indicate that the resorbed antecrystic (diopsidic and augitic) cores were remobilized between 5 ± 2 months and 6 ± 2 years before eruption, with the majority of crystals recording timescales between 2 and 6 years. These timescales correspond to approximately mid 2013 to mid-2017, with only one corresponding to early-mid 2019. While acknowledging we have only four estimates, we believe our calculated timescales are significant given the low proportion of antecrysts encountered in the 2019 products. Additionally, these timescales strongly contrast with those calculated for antecryst remobilisation during the 2003–2017 period3, which range from 4 to 182 years, including the paroxysmal events of 2003 and 2007 (Fig. 5a, b). Our data indicate that the majority of the 2019 clinopyroxene antecrysts were remobilised from the crystal mush by “older” mafic magma recharge during 2013-2017. Accordingly, there was very limited interaction between ascending _lp_-magma and crystal mush, or at least with the portion of the mush producing resorbed diopsidic and augitic antecrystic cores. This contrast sharply with the documented significant remobilisation of diopsidic antecryst cores in 2003–20173, suggesting that cannibalisation of the deeper diopsidic-dominated portion of the mushy system, generated a rejuvenated plumbing system where _lp-_magma is efficiently permeating the mush without significant clinopyroxene antecryst remobilisation (Fig. 7).

As discussed above, zoned 2 A diopsidic crystals indicate timescales from a few months to a few days (Fig. 5c, Supplementary Figs. 59–62, Table 1) and their cores may have higher Cr and Ni concentrations (Fig. 3), resorption textures (Supplementary Figs. 3 and 60–61) or normal compositional zoning (Supplementary Figs. 59–60). Rarely observed homogenous or slightly zoned diopside crystals (subgroup 2A) at Stromboli[42](/articles/s41467-022-35405-z#ref-CR42 "Francalanci, L. et al. Mineralogical, geochemical, and isotopic characteristics of the ejecta from the 5 april 2003 paroxysm at stromboli, italy: inferences on the preeruptive magma dynamics. In Geophysical Monograph Series, 331–345 (2008). https://doi.org/10.1029/182GM27

            ."),[45](/articles/s41467-022-35405-z#ref-CR45 "Bertagnini, A. et al. The Stromboli volcano: an integrated study of the. Erupt. Geophys. Monogr. Ser. 182, (2008)."), are interpreted to crystallise in the _lp_ magmatic environment, highlighting active crystallisation in a deeper reservoir where new crystal mush may build (Fig. [7](/articles/s41467-022-35405-z#Fig7)).

Implications for forecasting changes in eruptive style at basaltic volcanoes

Our data indicate an active role for the deeper mushy magmatic system in triggering the paroxysms at Stromboli in July–August 2019, not only through volatile supply19, but also through mafic magma ascent and recharge(s). Remarkable agreement exists between our calculated timescales for lp mafic magma injections into the shallow system and precursory monitoring signals before the two 2019 paroxysms.

The kernel density function for the 3 July paroxysm (red line and symbols in Fig. 8a) illustrates that mafic recharge started ~7 months pre-eruption with escalating activity from June 2019. Strikingly no corresponding timescales are recorded in the 28 August 2019 samples (blue line and symbols in Fig. 8a). Instead, clinopyroxenes predominantly record mafic recharge event(s) from ~4 July up to a few days before the 28 August paroxysm. This demonstrates complete evacuation of the shallow system, at least with respect to the _lp_-magma, which triggered the first paroxysm on 3 July 2019, before being refilled and re-destabilised by newly injected recharge _lp_-magma(s) soon after the first paroxysm and up to a few days before the 28 August event, eventually triggering the second paroxysm. This new model differs profoundly from previous models, that linked the 2019 paroxysms to gas flux blockage in the shallow conduit occurring a few minutes pre-eruption, with undegassed _lp_-recharge magma merely representing passive upward migration of a gas slug restoring the conduit gas flux17. From gas data, a mafic recharge trigger mechanism (bottom-up) was suggested for the 3 July paroxysm while a top-down mechanism (i.e., depressurisation due to the preceding lava flow and passive uprising of _lp_-magma) was proposed for the 28 August paroxysm19. By contrast, we provide strong evidence for the arrival of mafic magma(s) up to a few days before both paroxysms. In addition, the lack of overlapping short timescales between the 3 July and 28 August paroxysms suggests the arrival of distinct batches of gas-rich _lp_-magmas before the two events pointing to bottom-up mechanism also for the 28 August paroxysm.

Fig. 8: Relationship between timescales of mafic recharge events and monitoring data for the 2018- 2019 period.

a Kernel density estimates of timescales of magma recharge (subgroups 1B, 1C and Group 2; filled circles) and convective mixing in the shallow _hp-_reservoir (subgroup 1A; open circles) during the January 2018—August 2019 period for the two paroxysmal events (3 July, red curve and circles and 28 August blue curve and circles) are shown alongside the timescales determined for the major explosion on 25 June 2019 (black circles). Timescales of convective mixing for the 25 June major explosion are shown as single data (black open circle) due to the low density of data. b July 2018—September 2019 variation of excess of SO2 and CO2 flux (green and purple line respectively) as calculated by Aiuppa et al.19. January—September 2019 variation of hourly explosion frequency (orange line) data from Andronico et al.10. November 2018-September 2019 VLP size (amplitude in counts) time series (dashed cyan line) data from Giudicepietro et al.16. Dashed vertical light grey lines marked with light grey asterisks: major explosions between January 2018 and September 2019 (data from INGV[46](/articles/s41467-022-35405-z#ref-CR46 "Bollettini Multidisciplinari. https://www.ct.ingv.it/index.php/monitoraggio-e-sorveglianza/prodotti-del-monitoraggio/bollettini-settimanali-multidisciplinari

            .")); continuous vertical dark grey lines marked with red and grey starts: 3 July and 28 August 2019 paroxysmal events.

Full size image

The sequence of events highlighted here coincides well with monitoring data10,19 (Fig. 8b). Increased excess SO2 and CO2 flux19 (green and purple lines in Fig. 8b) correspond to an increasing density in the timescales of mafic injection recorded by clinopyroxenes from the 3 July paroxysm. The first peak of the kernel distribution for the 3 July event represents a timescale of 224 ± 98 days and subsequent peaks correspond to timescales of 165 ± 65 and 128 ± 52 days before eruption (Table 1). These latter timescales agree well with the increasing CO2 and SO2 fluxes ~180 days and ~120 days pre-eruption19. The emerging view is that, contrary to the previous model19, magma and gas remained coupled during several months in the build-up to the 3 July paroxysm until a few days before eruption. Chaotic mixing in the shallow reservoir, recorded by subgroup 1A crystals, contemporaneously with the mafic recharge, shows the perturbation of the shallow _hp-_reservoir by the arrival of deeper undegassed _lp_-magma. A similar scenario is also envisaged for the 28 August paroxysm (Fig. 8b).

The hourly explosion frequency10 (orange line in Fig. 8b), reveals that the sudden increased explosivity from February-March 2019 coincides with an increasing density of recharge timescales. Furthermore, the highest peaks in the kernel distributions recorded by clinopyroxenes match perfectly with the recorded periods of high to very high levels of normal activity both in terms of explosive frequency and intensity observed from 9 June 2019, ~1 month before the 3 July paroxysm, and lasting until 20 September 201910. The size and amplitude of VLP (very-long period) seismic signals (dashed blue line in Fig. 8b) significantly increased between November 2018 and January 2019, with a strong change in the VLP activity also occurring one month before the 3 July paroxysm16. The rapid response of the shallow system to the injections of _lp_-magma is interpreted to result from a more direct link between the two reservoirs, allowing the _lp-_magma permeating more efficiently the lower portion of the mush and resulting in increased explosive activity and changes in the VLP activity. Our model also agrees well with CO2 soil degassing18, which significantly increased since 2016 and culminated with the 2019 paroxysms, and with our observations[46](/articles/s41467-022-35405-z#ref-CR46 "Bollettini Multidisciplinari. https://www.ct.ingv.it/index.php/monitoraggio-e-sorveglianza/prodotti-del-monitoraggio/bollettini-settimanali-multidisciplinari

            .") of more variable eruptive activity since September 2019, where at least eight major explosions, intense spattering episodes and lava overflows have occurred.

In conclusion, our investigation of Stromboli’s 2019 paroxysms reveals a key role for batches of undegassed _lp_-magma recharge arriving in the shallow reservoir up to a few days before these events. Our data indicate a rejuvenated plumbing system where _lp-_magma permeates more efficiently the lower diopside-dominated portion of the mush allowing a direct linkage between deeper (lp) and shallow (hp) reservoirs and limited interaction between _lp-_recharge magma and extant crystal mush (Fig. 7). This sustains the current variability of eruptive styles with near immediate eruptive response to mafic magma recharge. The remarkable agreement between our calculated recharge timescales and the observed variation in time of various monitoring signals strongly supports such a model (Fig. 8).

The processes observed at Stromboli are likely common to other persistently active basaltic volcanoes (e.g., Etna, Italy;47 Volcan de Fuego, Guatemala;[48](/articles/s41467-022-35405-z#ref-CR48 "Liu, E. J. et al. Petrologic monitoring at Volcán de Fuego, Guatemala. (2020) https://doi.org/10.1016/j.jvolgeores.2020.107044

            .") Kilauea, Hawaii[49](/articles/s41467-022-35405-z#ref-CR49 "Patrick, M. R. et al. The cascading origin of the 2018 Kīlauea eruption and implications for future forecasting. Nat. Commun. 2020 111 11, 1–13 (2020).")) that show similar sudden shifts from mild Strombolian explosions to violent paroxysms. Stromboli’s post-2019 activity, characterised by an efficient connection between deeper and shallower reservoirs, has important implications for forecasting future paroxysmal eruptions and interpreting volcano monitoring data. We suggest that for well-monitored volcanoes, where detailed knowledge of the tempo and processes governing the plumbing system exists, devising a high-frequency petrological monitoring protocol can help discriminate gas-driven and magma-driven eruption triggering. Such information can be integrated with volcanological, geophysical, and geochemical data for a fast response during a volcanic unrest phase or crisis.

Methods

Analytical techniques

Bulk rock analyses for major and trace elements were conducted at Actlabs (Activation Laboratories Ltd., Ancaster, Ontario, Canada). Results are reported in Supplementary Table 1. Major elements were analysed by lithium metaborate/tetraborate fusion—ICP-OES (inductively coupled plasma optical emission spectrometry) and run on a Thermo Jarrell Ash ENVIRO II ICP. FeO was determined through titration, using cold acid digestion of ammonium metavanadate, and hydrofluoric acid in an open system. Ferrous ammonium sulphate was added after digestion and potassium dichromate was the titrating agent. Trace elements were measured by lithium metaborate/tetraborate fusion—ICP-MS (inductively coupled plasma mass spectrometry). Fused samples were diluted and analysed using Perkin Elmer Sciex ELAN 6000, 6100, or 9000 instruments. A bulk rock quality control table measured against certified standards used for bulk rock analyses at the Activation Laboratories Ltd. (Actlabs) is provided in the supplementary material (Supplementary Table 5).

Automated SEM-EDS maps were undertaken at Imaging and Analysis Centre (IAC) at the Natural History Museum in London (NHM) on a Tescan Integrated Mineral Analyser (TIMA) with four EDAX EDS detectors using a field emission source at 25 kV and 19.75 beam intensity. Analyses were carried out on a 5 μm pixel spacing with 1000 X-ray counts collected per pixel amounting to around 12 million analysis points and 12 billion X-ray counts per sample. Mineral classifications were manually assigned using only spectra acquired from this sample set.

Back Scattered electron (BSE) images were acquired with the Jeol IT500 FEG-SEM and the FEI QUANTA 650 FEG-SEM, at IAC-NHM, operating at 15–20 and 15 KeV electron current, and 10–11 and 15 mm WD (working distance), respectively. The spatial resolution of the SEM at our analytical condition was tested imaging the Richardson Test Slide (model 80302, serial n. 10461) for the 53 and 50 nm resolution. This test sample is routinely used to test the spatial resolution of SEM. We performed the test at different voltages (20, 15, 10 and 8 kV), working distances (WD) of 8 and 15 mm and different beam spot sizes. The images were all taken at the FEI Quanta 650 FEG SEM, where the spot size is a numerical value and we typically operate between 3 and 5. The results of the test are reported in Supplementary Figures 65–67, where it is visible that the resolution remains unchanged under all the different tested conditions. Therefore, at our operation conditions of 20–15 kV and 15–10 mm WD we have a resolution better than 50 nm, which is four times smaller than a normal pixel size of 0.2 μm and on average 20–40 times smaller than our smallest diffusion profiles (~1–2 μm).

To obtain the best results for diffusion modelling and minimise the effect of not sectioning perpendicular to the core-rim boundary, careful screening of the crystals in terms of reducing the orientation effect was carried out in the analysed thin sections. Diffusion profiles were performed along the direction perpendicular to the (100) or (010) crystallographic plane, according to ref. 40.

Electron microprobe Chemical composition of clinopyroxenes and matrix glasses alongside high-resolution core-rim compositional profile, ~3–5 µm step size, of major and selected trace elements of clinopyroxene were obtained using a Cameca SX 100 electron microprobe at IAC-NHM. The Cameca SX100 is equipped with five WDS (wavelength dispersive) spectrometers and one EDS (energy dispersive) spectrometer. It was operated at 20KeV and 20 nA with a focused beam for clinopyroxenes and a defocused 10 μm spot size for matrix glasses. Sodium was measured as first elements and counted for 10 s, all the other elements (Si, Mg, Al, Ca, Ti, Cr, Mn, Fe, Ca, and Ni) were measured for 20 s. Matrix effects were corrected using the Cameca PAP built-in protocol[50](/articles/s41467-022-35405-z#ref-CR50 "Pouchou, J.-L. & Pichoir, F. Quantitative analysis of homogeneous or stratified microvolumes applying the model “PAP”. In Electron Probe Quantitation pp 31–75 (1991). https://doi.org/10.1007/978-1-4899-2617-3_4

            ."). Detection limits for most elements were <0.2 wt.% and standard deviation of the analyses was typically <0.5 wt.%. A quality control table of measured against certified standard used for electron microprobe analyses of clinopyroxene at NHM-IAC laboratories is provided in the Supplementary Material (Supplementary Table [6](/articles/s41467-022-35405-z#MOESM8)). From each profile, representative analysis of the core-mantle-rim portion of the crystal were extracted and used to discuss the compositional variability and determine the _P-T-_H2O contents of clinopyroxenes. Data are reported in Supplementary Tables [2](/articles/s41467-022-35405-z#MOESM1) and [3](/articles/s41467-022-35405-z#MOESM1).

Trace element analysis Concentrations of trace elements in pyroxenes and glasses were determined by laser ablation inductively coupled plasma mass spectrometry (LA-ICPMS) spot analysis at the IAC-NHM and at the Institute of Geochemistry and Petrology, ETH Zürich. The analyses at NHM were performed using a 193 nm Iridia laser ablation system (Teledyne Photon Machines, Bozeman, MT, USA) coupled to an Agilent 8900 ICP-QQQ-MS. The Iridia is equipped with a 193 nm MLase short pulse (<7 ns) laser and the low-dispersion Cobalt (2-volume) ablation chamber51. The laser was operated using an 8 Hz repetition rate, a fluence of ~3.5 J cm−2 and a spot diameter of 35 µm. Samples were ablated in a He atmosphere (c. 0.5 l min−1) and the ablated aerosol was mixed with Ar (c. 0.8 l min−1) before being introduced to the interface of the ICP-MS. The ICP-MS was calibrated for maximum sensitivity, while minimising oxide production (232Th16O+/232Th+ <0.15 %) and laser induced fractionation (232Th+/238U+ > 90%) using a line scan of the NIST SRM 612 at the start of each analytical session. Analyses at ETH were conducted using a 193 nm Resonetics ArF excimer laser coupled to a Thermo Element XR ICPMS within. Analytical spot sizes were 43 µm with the energy of laser output maintained at ~3.5 J/cm2.

For all trace element determinations, a set of two NIST SRM 612 reference glass was measured every 25 unknowns to calibrate and monitor for instrumental drift. A further set of two GSD-1g reference glass was measured every 25 unknowns as part of a long-term lab reproducibility. Data reduction was performed using Iolite v. 3.6152 (NHM) and MATLAB-based programme SILLS53 (ETH) using the Ca concentration determined by electron microprobe as the internal standard. For each analysis a baseline of 25 s was measured prior to 40 s of ablation. The following isotopes were measured: 24Mg, 27Al, 43Ca, 45Sc, 47Ti, 51V, 52Cr, 55Mn, 60Ni, 66Zn, 89Y, 90Zr, 139La, 140Ce, 141Pr, 146Nd, 147Sm, 153Eu, 157Gd, 159Tb, 163Dy, 165Ho, 166Er, 169Tm, 172Yb, and 175Lu. Dwell times for each element was 10 ms, except for 24Mg, 27Al (5 ms); 52Cr, 60Ni, 66Zn, 153Eu, 159Tb, 165Ho, 169Tm, 175Lu (20 ms). All trace elemental determinations are considered to have precision better than 5% relative, based on homogeneous standard reproducibility where abundances are significantly above (typically 3x) the limits of detection. Elements of particular interest to this study (e.g. La, Yb, Zr, Sc, Cr, Fig. 6) all occur at abundances an order of magnitude or greater above typical limits of detection for these analytical conditions. Data are reported in Supplementary Table 4.

LA-ICP-MS trace element mapping of clinopyroxene grains was conducted over six analytical sessions at IAC-NHM (see above for instrument details). Samples were ablated in a He atmosphere (c. 0.8 l min−1) and the ablated aerosol was mixed with Ar (c. 0.8 l min−1) and introduced to the interface of the ICP-MS using the Aerosol Rapid Introduction System (ARIS; Teledyne Photon Machines, Bozeman, MT, USA; Van Acker et al., 2016). For all mapping experiments aerosol dispersion (washout time) was minimised by optimising the full width of a single laser pulse at 10% of the maximum peak intensity (FW0.1 M) to <5 ms using a live signal monitor (LA Monitor) in the HDIP software[54](/articles/s41467-022-35405-z#ref-CR54 "Teledyne CETAC Technologies. HDIP LA-ICP-MS Image Processing Software [WWW Document]. https://www.teledynecetac.com/products/laser-ablation/hdip

             (2020).") to improve signal sensitivity and reduce imaging time.

Laser ablation compositional maps were acquired through a series of adjacent line scans, obtained by moving the stage at constant speed under a fixed ablation site. The effect of analytical artefacts (aliasing) caused by the interaction of the pulsed laser signal and the total sweep cycle time55, was minimised by adjusting the total sweep cycle of the ICP-MS and laser parameters (frequency, scan speed, spot size) using the LA Optimise tool in HDIP (based on equations of Van Marderen et al.51,55 and Šala et al.56).

All mapping experiments were carried out using a 3 × 3 µm square beam, a fluence of 3.5 J cm−2, a repetition rate of 375 Hz and a scan speed of 93.75 µm s−1, corresponding to 12 laser pulses per 3 µm pixel. Six isotopes were measured in no-gas mode, with a total ICP-MS sweep cycle of 32 ms. A line scan of a representative pyroxene was used to optimise dwell times for each measured isotope to ensure sufficient signal/noise was achieved on the low abundance isotopes at the small spot size. Isotopes measured in order of increasing dwell time are: 24Mg, 27Al (0.5 ms); 47Ti (3 ms); 52Cr (5.2 ms); 43Ca (5.4 ms); 60Ni (7.4 ms). Each analytical session comprised the acquisition of pyroxene compositional maps during a single ICP-MS acquisition lasting between c. 2 and 9 hours. Individual compositional maps took between c. 30 minutes and c. 3.5 hours to gather. A gas blank was gathered between each successive line scan. NIST SRM 61257 was measured as line scans every ~30 minutes for calibration and drift correction. Larger images were separated into multiple blocks by reference materials to better constrain instrumental drift.

Image processing and quantification was carried out using HDIP (v. 1.6). After synchronising the timestamps of the laser with the ICP-MS signal, the background of each isotope was calculated during periods the laser was not firing and subtracted from each isotope measured. Drift correction was applied to the measured counts of each isotope using the NIST SRM 612. Pyroxene grains were isolated from each created image based on measured counts of 43Ca using the labelled segments tool in HDIP. Glass reference materials BCR-2g and GSD-1g were used as secondary reference material. Trace element concentrations in pyroxenes were calculated on a pixel-by-pixel basis using the Sum Normalisation tool in HDIP using the Ca concentration obtained from EMPA as an internal standard. For crystals zoned in Ca concentration each zone was isolated based on the BSE image using the labelled segments tool in HDIP and the appropriate Ca concentration was used as the internal standard for each zone.

Clinopyroxene calculation scheme and modelling

Following the calculation scheme of Putirka et al.58, the total number of oxygens in the formula of clinopyroxene is fixed to six. Then, it is assumed that X Si + X TAl = 2.00 in T-site, where X i is the cation fraction of the element of interest i. Jd (jadeite) = X M1Al if X M1Al < X_ Na or, alternatively, Jd = X Na. CaTs (Ca-Tschermark) = X M1Al – Jd, CaFeTs = (X M1Fe3+ + X Cr)/2, and CaTiTs = (X TAl – CaTs)/2 (if X TAl > CaTs). Di + Hd (diopside + hedenbergite) = X Ca – CaTs – CaFeTs – CaTiTs and En + Fs (enstatite + ferrosilite) = [(X TOTMg + X TOTFe) – (Di + Hd)]/2. A charge balance equation is applied to all the analyses for determining X Fe3+ [i.e., (X M1Al + X M1Fe3+ + X Cr + 2_X Ti) = (X TAl + X Na)]. X M1Mg and X M1Fe2+ are calculated as: X M1Mg = X TOTMg/(X TOTMg + X TOTFe) × (1 – X M1Al – X Cr – X Ti) and X M1Fe2+ = X TOTFe – X M2Fe2+, where X M2Fe2+ = 1 – (X Ca + X M2Mg) and X M2Mg = X TOTMg – X M1Mg, respectively.

The error associated with estimates of clinopyroxene-based thermometer (_P_-H2O-independent equation with error of ±6 °C[59](/articles/s41467-022-35405-z#ref-CR59 "Scarlato, P., Mollo, S., Petrone, C. M., Ubide, T. & Di Stefano, F. Interpreting Magma Dynamics Through a Statistically Refined Thermometer. In Crustal Magmatic System Evolution: Anatomy, Architecture, and Physico‐Chemical Processes (eds. Masotta, M., Beier, C. & Mollo, S.) 195–212 (AGU, Geophysical Mongraph Series, 2021). https://doi.org/10.1002/9781119564485.ch9

            ")), barometer (_T_\- H2O-independent equation with error of ±136 MPa[60](/articles/s41467-022-35405-z#ref-CR60 "Putirka, K., Johnson, M., Kinzler, R., Longhi, J. & Walker, D. Thermobarometry of mafic igneous rocks based on clinopyroxene-liquid equilibria, 0-30 kbar. Contrib. Miner. Pet. 123, 92–108 (1996).")), and hygrometer (_P-T_\-dependent equation with error of 0.5 wt.%[61](/articles/s41467-022-35405-z#ref-CR61 "Perinelli, C. et al. An improved clinopyroxene-based hygrometer for Etnean magmas and implications for eruption triggering mechanisms. Am. Mineral. 101, 2774–2777 (2016).")) is minimised by testing for equilibrium between clinopyroxene and melt (Supplementary Fig. [5](/articles/s41467-022-35405-z#MOESM1)). Poorly differentiated matrix glasses of _lp_\-products and their bulk rock analyses that are compositionally similar to these matrix glasses are equilibrated with diopside. Conversely, more evolved matrix glass compositions of _hp_\-products are equilibrated with augite. The use of bulk rock compositions is restricted to _lp_\-products, in agreement with the textural observation that their crystal content is only \~5 vol.% and implying that bulk rock and matrix glass compositions are almost identical. The most likely _P_\-_T_\-H2O conditions used for probability density functions are determined by a mathematical procedure based on minimisation of crystal-melt equilibrium parameters (Supplementary Fig. [5](/articles/s41467-022-35405-z#MOESM1)). Thus, barometric, thermometric, and hygrometric estimates are adjusted within the calibration errors of the predictive equations in order to minimise (1) the difference between measured vs. predicted diopside + hedenbergite components (_Δ_DiHd from Mollo et al.[62](/articles/s41467-022-35405-z#ref-CR62 "Mollo, S., Putirka, K., Misiti, V., Soligo, M. & Scarlato, P. A new test for equilibrium based on clinopyroxene–melt pairs: clues on the solidification temperatures of Etnean alkaline melts at post-eruptive conditions. Chem. Geol. 352, 92–100 (2013).")) and (2) the difference between measured vs. predicted clinopyroxene-melt partition coefficients for Na and Ti (_Δ_DNa and _Δ_DTi from Blundy et al.[63](/articles/s41467-022-35405-z#ref-CR63 "Blundy, J. D., Falloon, T. J., Wood, B. J. & Dalton, J. A. Sodium partitioning between clinopyroxene and silicate melts. J. Geophys. Res. Solid Earth 100, 15501–15515 (1995).") and Mollo et al.[33](/articles/s41467-022-35405-z#ref-CR33 "Mollo, S. et al. An integrated P-T-H2O-lattice strain model to quantify the role of clinopyroxene fractionation on REE+Y and HFSE patterns of mafic alkaline magmas: Application to eruptions at Mt. Etna. Earth-Sci. Rev. 185, 32–56 (2018)."), respectively). By integrating these thermodynamic-based models (Supplementary Fig. [5](/articles/s41467-022-35405-z#MOESM1)), the near-equilibrium state of clinopyroxene is verified when values of _Δ_DiHd, _Δ_DNa, and _Δ_DTi approaches the ideal equilibrium condition of _Δ_ \= 0 (cf. Mollo et al.[33](/articles/s41467-022-35405-z#ref-CR33 "Mollo, S. et al. An integrated P-T-H2O-lattice strain model to quantify the role of clinopyroxene fractionation on REE+Y and HFSE patterns of mafic alkaline magmas: Application to eruptions at Mt. Etna. Earth-Sci. Rev. 185, 32–56 (2018).")).

Fe-Mg elemental diffusion modelling

Timescales of mafic recharges, convective mixing and antecrysts remobilisation from the crystal mush were determined using NIDIS (Non-Isothermal Diffusion Incremental Step) model of Petrone et al.40. The NIDIS model considers how temperature changes during the growth of a zoned crystal modulate the rate of diffusion, allowing to calculate the storage time in different magmatic environments characterised by different temperatures and chemical compositions. In the case of clinopyroxenes, the model uses grey values extracted from BSE images as proxy for Mg# (i.e., Fe-Mg elemental diffusion) and the diffusive profile for each compositional boundary is modelled at the temperature of the liquid in equilibrium with that specific composition40. For crystals like those of subgroup 1C (Fig. 2, Supplementary Figs. 57–58), the NIDIS model calculates the total residence time (Δt = Δt 1 + Δt 2, labelled as diopsidic band in Fig. 5c) by estimating the temporal interval between the formation of the diopsidic band (Δt 1) and the eruption (Δt 2). The partial timescale Δt 1 corresponds to the partial diffusive time after the formation of the internal compositional boundary (i.e., augitic core—diopsidic band) at the higher T of diopsidic band. Δt 2 corresponds to the partial diffusive time after the formation of the external compositional boundary (i.e., diopsidic band —augitic overgrowth) until eruption, and it is calculated at the T of augitic overgrowth. According to NIDIS model, Δt is calculated using a backward approach starting from the last formed compositional boundary (Δt 2) and proceeding inward towards the more internal boundary (Δt 1, see Petrone et al.30,40 for further details). Residence time in the last magmatic environment before eruption was calculated for subgroups 1B and 2B crystals (Figs. 2 and 5, Supplementary 15–55 and 63–64, Table 1), adapting the NIDIS model30,40. The overgrowth of crystals records the last magmatic event prior to eruption and in the case of subgroup 1B crystals the diopsidic overgrowth suggests a mafic recharge prior to eruption. It has been argued30 that if the mafic recharge event occurred shortly (i.e., from a few days to a few months) prior to eruption, it represents the immediate eruption trigger, pointing to a bottom-up mechanism (i.e., an active role of the deeper system as opposed to a passive role driven by the shallow conduit dynamics19). The diopsidic overgrowths of subgroup 1B clinopyroxenes record the timescales of mafic recharge event(s), as confirmed by the core-rim compositional profiles (Supplementary Figs. 15–55) where the compositional changes of Ca and Al, the slow diffusive elements, preserve sharp variations suggestive of diffusion profiles due to magma mixing events as opposed to chemical variation due to crystal fractionation (see § Mafic recharge as eruption trigger of 2019 paroxysms). The augitic overgrowth of 2B events provided time constraints on remobilisation of resorbed diopsidic antecryst cores from the crystal mush and their subsequent storage in the shallow hp magmatic environment (see § The lifecycle of basaltic crystal mush). Time constraint on antecryst remobilisation form the crystal mush have been also obtained from resorbed augitic cores surrounded by diopsidic band of some subgroup 1C crystals (Supplementary Figs. 56–57). As discussed in the text (see § Convective mixing regime controlling magma dynamics), we have calculated pre-eruptive timescales for subgroups 1A and 2A crystals (Figs. 2 and 5c, Supplementary Figs. 6–14, 59–62). These diffusive timescales set constraints on the convective mixing processes in both the shallow hp and the deeper mafic lp magmatic environments.

The greyscales values of zoned boundaries of clinopyroxenes were extracted by BSE-photomicrographs using the greyvalues script40, calibrated against the Fe-Mg (Mg#) chemical profile acquired at the electron microprobe and fitted with the cratefit script40 to calculate the timescale. Electron beam convolution effect can be important particularly for shorth diffusion profile (i.e., short timescales). To assess if the greyscale profiles were at least in part modified by convolution effect, nine clinopyroxenes with timescales ranging from a few days to years were re-imaged at the FEG-SEM collecting BSE images at 10 kV and 8 kV with the working distance set at 8 mm. These data were compared with data extracted from BSE images taken at higher voltages (15–20 kV) and working distances (10–15 mm) (i.e., the analytical conditions used in this study). Each crystal was imaged under the new conditions (3–6 photos and timescales), but we carefully selected the same area from which the original profile was extracted, reproducing a similar greyscale profile to avoid any further bias. For simplicity we report the calculated timescales at the different conditions (the goodness of the fit of the greyscale profiles measured by the r2 were all between 1.00 and 0.99, except for 3 profiles that returned 0.98–0.97). As shown in Supplementary figures 68–69, all timescales calculated for each crystal under different analytical conditions agree within the analytical error of the calculated timescale. Therefore, we can exclude any smoothing of the profile due to convolution effect.

Grey values were acquired along transect normal to the {1 0 0} and {0 1 0} sectors avoiding the c axis. According to thermobarometric results and considering 2-sigma standard error of ±16 °C on the estimated temperatures (Fig. 4) Fe-Mg diffusion coefficients were calculated at 1155 and 1175 °C for augitic and diopsidic composition respectively. The values of the pre-exponential factor and activation energy are from Dimanov & Sautter64. As the diffusion temperature decreases, the timescale error increases due to variations in activation energy and intra-crystal cation diffusivity40,[65](/articles/s41467-022-35405-z#ref-CR65 "Petrone, C. M. & Mangler, M. F. Elemental Diffusion Chronostratigraphy. In: Crustal Magmatic System Evolution: Anatomy, Architecture, and Physico‐Chemical Processes (eds. Masotta, M., Beir, C. & Mollo, S.) 177–193 (AGU, Geophysical Monograph Series, 2021). https://doi.org/10.1002/9781119564485.ch8

            .") (and references therein). In light of these considerations, timescales from this study represent conservative crystal residence times for the magma dynamics at Stromboli. As discussed in Petrone et al.[30](/articles/s41467-022-35405-z#ref-CR30 "Petrone, C. M., Braschi, E., Francalanci, L., Casalini, M. & Tommasini, S. Rapid mixing and short storage timescale in the magma dynamics of a steady-state volcano. Earth Planet. Sci. Lett. 492, 206–221 (2018)."), under the used analytical condition and considering the uncertainty on the fit parameter \\(\\sqrt{4{Dt}}\\) (always better than 5% and in most cases better than 2%), the pixel resolution of the BSE-SEM greyscale (better than 0.3 μm) and the values of the diffusion coefficient D, the NIDIS method allows to start detecting change of the diffusion profile from the initial condition after only 1 day (see Supplementary Fig. [6](/articles/s41467-022-35405-z#MOESM1) in Petrone et al.[30](/articles/s41467-022-35405-z#ref-CR30 "Petrone, C. M., Braschi, E., Francalanci, L., Casalini, M. & Tommasini, S. Rapid mixing and short storage timescale in the magma dynamics of a steady-state volcano. Earth Planet. Sci. Lett. 492, 206–221 (2018).")).

Data availability

The entire set of mineral chemistry data is available at NERC EDS Dataset Petrone CM 2022[66](/articles/s41467-022-35405-z#ref-CR66 "Petrone, C. M. Chemical and Mineralogical analysis, SEM images and elemental diffusion timescales. NERC EDS Natl Geosci. Data Cent. (Dataset). https://doi.org/10.5285/93f8a6f4-4f99-41b5-a1bd-107a438b2cd7

             (2022)."). All other data are provided in the paper and/or the [Supplementary Materials](/articles/s41467-022-35405-z#MOESM1). All clinopyroxenes BSE images, core-rim chemical profiles and calculated timescales are provided in the [Supplementary Materials](/articles/s41467-022-35405-z#MOESM1).

References

1. Di Renzo, V., Corsaro, R. A., Miraglia, L., Pompilio, M. & Civetta, L. Long and short-term magma differentiation at Mt. Etna as revealed by Sr-Nd isotopes and geochemical data. Earth Sci. Rev. 190, 112–130 (2019).
   Article ADS Google Scholar
2. Maclennan, J. Mafic tiers and transient mushes: evidence from Iceland. Philos. Trans. R. Soc. A Math. Phys. Eng. Sci. 377, 20180021 (2019).
3. Di Stefano, F. et al. Mush cannibalism and disruption recorded by clinopyroxene phenocrysts at Stromboli volcano: new insights from recent 2003–2017 activity. Lithos 360–361, 105440 (2020).
   Article Google Scholar
4. Bevilacqua, A. et al. Major explosions and paroxysms at Stromboli (Italy): a new historical catalog and temporal models of occurrence with uncertainty quantification. Sci. Rep. 10, 1–18 (2020).
   Article Google Scholar
5. Métrich, N., Bertagnini, A. & Pistolesi, M. Paroxysms at Stromboli Volcano (Italy): source, genesis and dynamics. Front. Earth Sci. 9, 1–17 (2021).
   Article Google Scholar
6. Barberi, F. Volcanic hazard assessment at Stromboli based on review of historical data. Acta Vulcanol. 3, 173–187 (1993).
   Google Scholar
7. Calvari, S. et al. Lava effusion-a slow fuse for paroxysms at Stromboli volcano? Earth Planet. Sci. Lett. 301, 317–323 (2011).
   Article ADS CAS Google Scholar
8. Ripepe, M. et al. Forecasting effusive dynamics and decompression rates by magmastatic model at open-vent volcanoes. Sci. Rep. 7, 1–9 (2017).
   Article CAS Google Scholar
9. Di Traglia, F. et al. The 2014 effusive eruption at stromboli: new insights from in situ and remote-sensing measurements. Remote Sens. 10, 2035 (2018).
   Article ADS Google Scholar
10. Andronico, D. et al. Uncovering the eruptive patterns of the 2019 double paroxysm eruption crisis of Stromboli volcano. Nat. Commun. 12, 4213 (2021).
11. Giordano, G. & De Astis, G. The summer 2019 basaltic Vulcanian eruptions (paroxysms) of Stromboli. Bull. Volcanol. 83, 1 (2021).
12. Calvari, S. et al. Variable magnitude and intensity of strombolian explosions: focus on the eruptive processes for a first classification scheme for stromboli volcano (italy). Remote Sens. 13, 1–30 (2021).
    Article Google Scholar
13. Ripepe, M. et al. Ground deformation reveals the scale-invariant conduit dynamics driving explosive basaltic eruptions. Nat. Commun. 12, 1–8 (2021).
    Article Google Scholar
14. Plank, S. et al. The July/August 2019 Lava Flows at the Sciara del Fuoco, stromboli-analysis from multi-sensor infrared satellite imagery. Remote Sens. 11, 2879 (2019).
15. Stromboli Daily Bulletin. http://lgs.geo.unifi.it/index.php/reports/stromboli-daily.
16. Giudicepietro, F. et al. Geophysical precursors of the July-August 2019 paroxysmal eruptive phase and their implications for Stromboli volcano (Italy) monitoring. Sci. Rep. 10, 1–16 (2020).
    Article Google Scholar
17. Viccaro, M. et al. Shallow conduit dynamics fuel the unexpected paroxysms of Stromboli volcano during the summer 2019. Sci. Rep. 11, 1–15 (2021).
    Article Google Scholar
18. Inguaggiato, S., Vita, F., Cangemi, M., Inguaggiato, C. & Calderone, L. The monitoring of CO2 soil degassing as indicator of increasing volcanic activity: the paroxysmal activity at stromboli volcano in 2019-2021. Geosci 11, 169 (2021).
    Article ADS CAS Google Scholar
19. Aiuppa, A. et al. Volcanic CO2 tracks the incubation period of basaltic paroxysms. Sci. Adv. 7, 191–208 (2021).
    Article ADS Google Scholar
20. Patanè, D. et al. The shallow magma chamber of Stromboli Volcano (Italy). Geophys. Res. Lett. 44, 6589–6596 (2017).
    Article ADS Google Scholar
21. Metrich, N., Bertagnini, A. & Muro, A. D. 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 (2010).
22. Francalanci, L., Tommasini, S. & Conticelli, S. 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. Geotherm. Res. 131, 179–211 (2004).
    Article ADS CAS Google Scholar
23. Felice, S. L. & Landi, P. The 2009 paroxysmal explosions at Stromboli (Italy): magma mixing and eruption dynamics. Bull. Volcanol. 73, 1147–1154 (2011).
    Article ADS Google Scholar
24. Pioli, L., Pistolesi, M. & Rosi, M. Transient explosions at open-vent volcanoes: the case of Stromboli (Italy). Geology 42, 863–866 (2014).
    Article ADS Google Scholar
25. Métrich, N., Bertagnini, A., Landi, P. & Rosi, M. Crystallization driven by decompression and water loss at Stromboli volcano (Aeolian Islands, Italy). J. Pet. 42, 1471–1490 (2001).
    Article ADS Google Scholar
26. Pompilio, M., Bertagnini, A. & Métrich, N. Geochemical heterogeneities and dynamics of magmas within the plumbing system of a persistently active volcano: evidence from Stromboli. Bull. Volcanol. 74, 881–894 (2012).
    Article ADS Google Scholar
27. Francalanci, L., Lucchi, F., Keller, J., De Astis, G. & Tranne, C. A. Eruptive, volcano-tectonic and magmatic history of the Stromboli volcano (north-eastern Aeolian archipelago). Geol. Soc. Mem. 37, 397–471 (2013).
    Article Google Scholar
28. Francalanci, L. et al. Intra-grain Sr isotope evidence for crystal recycling and multiple magma reservoirs in the recent activity of Stromboli volcano, Southern Italy. J. Pet. 46, 1997–2021 (2005).
    Article ADS CAS Google Scholar
29. Landi, P., Métrich, N., Bertagnini, A. & Rosi, M. Dynamics of magma mixing and degassing recorded in plagioclase at Stromboli (Aeolian Archipelago, Italy). Contrib. Mineral. Petrol. 147, 213–227 (2004).
    Article ADS CAS Google Scholar
30. Petrone, C. M., Braschi, E., Francalanci, L., Casalini, M. & Tommasini, S. Rapid mixing and short storage timescale in the magma dynamics of a steady-state volcano. Earth Planet. Sci. Lett. 492, 206–221 (2018).
    Article ADS CAS Google Scholar
31. Blundy, J. & Cashman, K. Petrologic reconstruction of magmatic system variables and processes. Rev. Mineral. Geochem. 69, 179–239 (2008).
    Article CAS Google Scholar
32. Putirka, K. D. Thermometers and barometers for volcanic systems. Rev. Mineral. Geochem. 69, 61–120 (2008).
    Article CAS Google Scholar
33. Mollo, S. et al. An integrated P-T-H2O-lattice strain model to quantify the role of clinopyroxene fractionation on REE+Y and HFSE patterns of mafic alkaline magmas: Application to eruptions at Mt. Etna. Earth-Sci. Rev. 185, 32–56 (2018).
    Article ADS CAS Google Scholar
34. Clocchiatti, R. La transition augite diopside et les liquides silicates intra cristallins dans les pyroclastes de l’activite actuelle du Stromboli: temoignages de la reinjection et du melange magmatiques. Bull. Volcanol. 44, 339–357 (1981).
35. Bertagnini, A., Métrich, N., Landi, P. & Rosi, M. Stromboli volcano (Aeolian Archipelago, Italy): an open window on the deep-feeding system of a steady state basaltic volcano. J. Geophys. Res. Solid Earth 108, 2336 (2003).
36. Bragagni, A., Avanzinelli, R., Freymuth, H. & Francalanci, L. Recycling of crystal mush-derived melts and short magma residence times revealed by U-series disequilibria at Stromboli volcano. Earth Planet. Sci. Lett. 404, 206–219 (2014).
    Article ADS CAS Google Scholar
37. Di Fiore, F. et al. Kinetic partitioning of major and trace cations between clinopyroxene and phonotephritic melt under convective stirring conditions: New insights into clinopyroxene sector zoning and concentric zoning. Chem. Geol. 584, 120531 (2021).
    Article Google Scholar
38. Francalanci, L. et al. Crystal recycling in the steady-state system of the active Stromboli volcano: a 2.5-ka story inferred from in situ Sr-isotope and trace element data. Contrib. Mineral. Petrol. 163, 109–131 (2012).
    Article ADS CAS Google Scholar
39. Landi, P. et al. Magma dynamics during the 2007 Stromboli eruption (Aeolian Islands, Italy): mineralogical, geochemical and isotopic data. J. Volcanol. Geotherm. Res. 182, 255–268 (2009).
    Article ADS CAS Google Scholar
40. Petrone, C. M., Bugatti, G., Braschi, E. & Tommasini, S. Pre-eruptive magmatic processes re-timed using a non-isothermal approach to magma chamber dynamics. Nat. Commun. 7, 12946 (2016).
    Article ADS CAS Google Scholar
41. Métrich, N., Bertagnini, A., Landi, P., Rosi, M. & Belhadj, O. Triggering mechanism at the origin of paroxysms at Stromboli (Aeolian Archipelago, Italy): the 5 April 2003 eruption. Geophys. Res. Lett. 32, 1–4 (2005).
    Article Google Scholar
42. Francalanci, L. et al. Mineralogical, geochemical, and isotopic characteristics of the ejecta from the 5 april 2003 paroxysm at stromboli, italy: inferences on the preeruptive magma dynamics. In Geophysical Monograph Series, 331–345 (2008). https://doi.org/10.1029/182GM27.
43. Speranza, F., Pompilio, M., Caracciolo, F. D. A. & Sagnotti, L. Holocene eruptive history of the Stromboli volcano: constraints from paleomagnetic dating. J. Geophys. Res. Solid Earth 113, 1–23 (2008).
    Article Google Scholar
44. Aiuppa, A. et al. Unusually large magmatic CO2 gas emissions prior to a basaltic paroxysm. Geophys. Res. Lett. 37, 1–5 (2010).
    Article Google Scholar
45. Bertagnini, A. et al. The Stromboli volcano: an integrated study of the. Erupt. Geophys. Monogr. Ser. 182, (2008).
46. Bollettini Multidisciplinari. https://www.ct.ingv.it/index.php/monitoraggio-e-sorveglianza/prodotti-del-monitoraggio/bollettini-settimanali-multidisciplinari.
47. Giuffrida, M. et al. Tracking the summit activity of Mt. Etna volcano between July 2019 and January 2020 by integrating petrological and geophysical data. J. Volcanol. Geotherm. Res. 418, 107350 (2021).
    Article CAS Google Scholar
48. Liu, E. J. et al. Petrologic monitoring at Volcán de Fuego, Guatemala. (2020) https://doi.org/10.1016/j.jvolgeores.2020.107044.
49. Patrick, M. R. et al. The cascading origin of the 2018 Kīlauea eruption and implications for future forecasting. Nat. Commun. 2020 111 11, 1–13 (2020).
    Google Scholar
50. Pouchou, J.-L. & Pichoir, F. Quantitative analysis of homogeneous or stratified microvolumes applying the model “PAP”. In Electron Probe Quantitation pp 31–75 (1991). https://doi.org/10.1007/978-1-4899-2617-3_4.
51. Van Malderen, S. J. M., Acker, T., Van & Vanhaecke, F. Sub-micrometer nanosecond LA-ICP-MS imaging at pixel acquisition rates above 250 Hz via a low-dispersion setup. Anal. Chem. 92, 5756–5764 (2020).
    Article Google Scholar
52. Paton, C., Hellstrom, J., Paul, B., Woodhead, J. & Hergt, J. Iolite: freeware for the visualisation and processing of mass spectrometric data. J. Anal. Spectrom. 26, 2508–2518 (2011).
    Article CAS Google Scholar
53. Guillong, M., Meier, D. L., Allan, M. M., Heinrich, C. A. & Yardley, B. W. D. SILLS: a matlab-based program for the reduction of laser ablation ICP–MS data of homogeneous materials and inclusions. Mineral. Assoc. Can. Short. Course 40, 328–333 (2008).
    CAS Google Scholar
54. Teledyne CETAC Technologies. HDIP LA-ICP-MS Image Processing Software [WWW Document]. https://www.teledynecetac.com/products/laser-ablation/hdip (2020).
55. Van Malderen, S. J. M., van Elteren, J. T., Šelih, V. S. & Vanhaecke, F. Considerations on data acquisition in laser ablation-inductively coupled plasma-mass spectrometry with low-dispersion interfaces. Spectrochim. Acta Part B . Spectrosc. 140, 29–34 (2018).
    Article ADS Google Scholar
56. Šala, M., Šelih, V. S., Stremtan, C. C., Tamas, T. & Van Elteren, J. T. Implications of laser shot dosage on image quality in LA-ICP-QMS imaging. J. Anal. Spectrom. 36, 75–79 (2021).
    Article Google Scholar
57. Jochum, K. P. et al. Determination of reference values for NIST SRM 610–617 glasses following ISO guidelines. Geostand. Geoanalytical Res 35, 397–429 (2011).
    Article CAS Google Scholar
58. Putirka, K. Melting depths and mantle heterogeneity beneath Hawaii and the East Pacific rise: constraints from Na/Ti and rare earth element ratios. J. Geophys. Res. Solid Earth 104, 2817–2829 (1999).
    Article CAS Google Scholar
59. Scarlato, P., Mollo, S., Petrone, C. M., Ubide, T. & Di Stefano, F. Interpreting Magma Dynamics Through a Statistically Refined Thermometer. In _Crustal Magmatic System Evolution: Anatomy, Architecture, and Physico‐Chemical Proc_esses (eds. Masotta, M., Beier, C. & Mollo, S.) 195–212 (AGU, Geophysical Mongraph Series, 2021). https://doi.org/10.1002/9781119564485.ch9
60. Putirka, K., Johnson, M., Kinzler, R., Longhi, J. & Walker, D. Thermobarometry of mafic igneous rocks based on clinopyroxene-liquid equilibria, 0-30 kbar. Contrib. Miner. Pet. 123, 92–108 (1996).
    Article ADS CAS Google Scholar
61. Perinelli, C. et al. An improved clinopyroxene-based hygrometer for Etnean magmas and implications for eruption triggering mechanisms. Am. Mineral. 101, 2774–2777 (2016).
    Article ADS Google Scholar
62. Mollo, S., Putirka, K., Misiti, V., Soligo, M. & Scarlato, P. A new test for equilibrium based on clinopyroxene–melt pairs: clues on the solidification temperatures of Etnean alkaline melts at post-eruptive conditions. Chem. Geol. 352, 92–100 (2013).
    Article ADS CAS Google Scholar
63. Blundy, J. D., Falloon, T. J., Wood, B. J. & Dalton, J. A. Sodium partitioning between clinopyroxene and silicate melts. J. Geophys. Res. Solid Earth 100, 15501–15515 (1995).
    Article CAS Google Scholar
64. Dimanov, A. & Sautter, V. “Average” interdiffusion of (Fe,Mn)-Mg in natural diopside. Eur. J. Mineral. 12, 749–760 (2000).
    Article ADS CAS Google Scholar
65. Petrone, C. M. & Mangler, M. F. Elemental Diffusion Chronostratigraphy. In: _Crustal Magmatic System Evolution: Anatomy, Architecture, and Physico‐Chemical Proc_esses (eds. Masotta, M., Beir, C. & Mollo, S.) 177–193 (AGU, Geophysical Monograph Series, 2021). https://doi.org/10.1002/9781119564485.ch8.
66. Petrone, C. M. Chemical and Mineralogical analysis, SEM images and elemental diffusion timescales. NERC EDS Natl Geosci. Data Cent. (Dataset). https://doi.org/10.5285/93f8a6f4-4f99-41b5-a1bd-107a438b2cd7 (2022).
    Article Google Scholar
67. Mollo, S., Giacomoni, P. P., Andronico, D. & Scarlato, P. Clinopyroxene and titanomagnetite cation redistributions at Mt. Etna volcano (Sicily, Italy): footprints of the final solidification history of lava fountains and lava flows. Chem. Geol. 406, 45–54 (2015).
    Article ADS CAS Google Scholar
68. Mollo, S., Ubide, T., Di Stefano, F., Nazzari, M. & Scarlato, P. Polybaric/polythermal magma transport and trace element partitioning recorded in single crystals: a case study of a zoned clinopyroxene from Mt. Etna. Lithos 356–357, 105382 (2020).
    Article Google Scholar
69. Tommasini, S. et al. Critical assessment of pressure estimates in volcanic plumbing systems: The case study of Popocatépetl volcano, Mexico. Lithos 408–409, 106540 (2022).
    Article Google Scholar

Download references

Acknowledgements

The authors would like to thank Will Brownscombe for his help with the TIMA; Callum Hatch for his help with polished section preparations; Alex Ball for providing knowledge and experience at SEM; John Spratt, Innes Clatworthy and Alex Ball for allowing access to the IAC-NHM laboratories during the Covid-19 pandemic; Hilary Downes for useful discussions. Figure 7a has been reprinted from Di Stefano, F. et al., Mush cannibalism and disruption recorded by clinopyroxene phenocrysts at Stromboli volcano: New insights from recent 2003–2017 activity, Lithos 360–361, 105440 (2020) with permission from Elsevier. Figure 7a has also been also modified from Di Stefano et al.3 with permission from Elsevier. The authors acknowledge the following funding bodies: Natural Environment Research Council UK grant NE/T009292/1 to C.M.P. and R.G.; INGV - Progetti Ricerca Libera 2019 Grant #52/2020 to P.S.; PRIN MIUR Grant #2017J277S9_004 to E.D.B.; INGV Departmental Strategic Project UNO to P.S.; and Swiss National Science Foundation grant 200020-197040 to B.E.

Author information

Authors and Affiliations

1. Natural History Museum, Volcano Petrology Group, Cromwell Road, SW7 5BD, London, UK
   Chiara Maria Petrone
2. Department of Earth Sciences, Sapienza – University of Rome, P.Le Aldo Moro 5, 00185, Roma, Italy
   Silvio Mollo
3. Istituto Nazionale di Geofisica e Vulcanologia – Sezione di Roma 1, Via di Vigna Murata 605, 00143, Roma, Italy
   Silvio Mollo, Piergiorgio Scarlato, Elisabetta Del Bello, Alessio Pontesilli & Gianfilippo De Astis
4. School of Geography, Geology and the Environment, Keele University, Keele, Staffordshire, ST5 5BG, UK
   Ralf Gertisser
5. Natural History Museum, Core Research Laboratory, Cromwell Road, SW7 5BD, London, UK
   Yannick Buret
6. Istituto Nazionale di Geofisica e Vulcanologia, Osservatorio Etneo, Sezione di Catania, Piazza Roma 2, 95123, Catania, Italy
   Daniele Andronico
7. Institute of Geochemistry and Petrology, ETH Zurich, 8092, Zurich, Switzerland
   Ben Ellis
8. Department of Earth Sciences, University of Pisa, Via Santa Maria 53, 56126, Pisa, Italy
   Pier Paolo Giacomoni
9. Department of Physics and Earth Sciences, University of Ferrara, Via Saragat 1, 44122, Ferrara, Italy
   Massimo Coltorti
10. Department of Earth and Environmental Sciences, University of Iowa, Iowa City, IA, 52242, USA
    Mark Reagan

Authors

1. Chiara Maria Petrone
   You can also search for this author inPubMed Google Scholar
2. Silvio Mollo
   You can also search for this author inPubMed Google Scholar
3. Ralf Gertisser
   You can also search for this author inPubMed Google Scholar
4. Yannick Buret
   You can also search for this author inPubMed Google Scholar
5. Piergiorgio Scarlato
   You can also search for this author inPubMed Google Scholar
6. Elisabetta Del Bello
   You can also search for this author inPubMed Google Scholar
7. Daniele Andronico
   You can also search for this author inPubMed Google Scholar
8. Ben Ellis
   You can also search for this author inPubMed Google Scholar
9. Alessio Pontesilli
   You can also search for this author inPubMed Google Scholar
10. Gianfilippo De Astis
    You can also search for this author inPubMed Google Scholar
11. Pier Paolo Giacomoni
    You can also search for this author inPubMed Google Scholar
12. Massimo Coltorti
    You can also search for this author inPubMed Google Scholar
13. Mark Reagan
    You can also search for this author inPubMed Google Scholar

Contributions

C.M.P. conceived, coordinated, wrote, and organised the manuscript, figures and tables with contributions from S.M., R.G. and B.E. Sample collection was carried out by P.S., E.D.B., D.A., G.D.A. and P.G. Y.B., B.E., A.P. and C.M.P. carried out LA-ICPMS trace elements analysis on clinopyroxenes and Y.B. also carried out LA-ICPMS chemical map. C.M.P. carried out clinopyroxene and glass chemical analysis via EPMA, BSE imaging and elemental diffusion modelling. Thermobarometric and hygrometric calculation were carried out by S.M. All authors participated in the interpretation of results and finalisation of the manuscript.

Corresponding author

Correspondence toChiara Maria Petrone.

Ethics declarations

Competing interests

The authors declare no competing interests

Peer review

Peer review information

Nature Communications thanks Eniko Bali and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.

Additional information

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Rights and permissions

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 license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license 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 license, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and permissions

About this article

Cite this article

Petrone, C.M., Mollo, S., Gertisser, R. et al. Magma recharge and mush rejuvenation drive paroxysmal activity at Stromboli volcano.Nat Commun 13, 7717 (2022). https://doi.org/10.1038/s41467-022-35405-z

Download citation

• Received: 27 March 2022
• Accepted: 30 November 2022
• Published: 13 December 2022
• DOI: https://doi.org/10.1038/s41467-022-35405-z