Paroxysms at Stromboli Volcano (Italy): Source, Genesis and Dynamics

Deciphering the triggering mechanisms of violent explosive activity is of broad interest for understanding the dynamics of basaltic open-vent volcanic systems. For nearly 1300 years Stromboli has been renowned not only for its continuous degassing activity and mild explosions at the summit craters, but also for short-lived, violent explosive events of variable scale, known as major explosions and paroxysms. Here, we focus on the 1456 and 1930 paroxysms and on the most recent events, in July and August 2019 at Stromboli. We show that shallow phenomena such as flank collapse, lava outpouring through fractures opening, or partial emptying of the shallow conduit, only speed up volatile-rich magma ascent by increasing the decompression rate, whereas pressurization of the crustal system and the deep refilling by magma and its CO2-rich gas phase play a major role in triggering paroxysms. Moreover, we present new data on the geochemistry of the 2019 bulk pumice, along with a compilation of data from the literature, chemical profiles in olivine crystals, and the physical parameters of explosive eruptions of wide ranging magnitude and intensity. For small and large paroxysms, timescales were derived from Fe–Mg diffusion profiles in olivine. In both types of explosion, the last phases of crystallization-diffusion indicate rapid magma ascent rates of two to ten days prior to eruption. Trace element concentrations (Nb, La and Ba) and ratios (Rb/Th) indicate that the 2019 pumice samples plot in the domain of magma batches erupted within the last 20 years at Stromboli. As a whole, there is no correlation between magma geochemistry and magnitude or intensity of explosive eruptions, which span a range of ∼3 orders of magnitude (from major explosions to large paroxysms) based on estimates of erupted tephra volumes. In contrast, olivine compositions are a good proxy for erupted tephra volumes and magma flux. The correlation among physical and chemical parameters, which is valid for the overall spectrum of eruptions, implies that the magmatic source ultimately controls eruptive dynamics.


INTRODUCTION
Long-lasting, basaltic open-vent volcanoes are characterized by eruptive plumes and persistent gas emissions, explosive activity at the craters and/or emissions of lava flows resulting in hazards of variable intensity (Rose et al., 2013). Paroxysmal eruptions may occur unexpectedly and represent one of the major threats (e.g., the 2015 paroxysm at Villarrica volcano; Aiuppa et al., 2017;Johnson et al., 2018). Among basaltic open-vent systems, Stromboli is persistently active since the eighth century CE, and its Strombolian regime is occasionally interrupted by short-lived, violent explosive eruptions Rosi et al., 2013). Hence, understanding the mechanisms driving the paroxysmal eruptions of variable magnitude at Stromboli is of broad interest.
Low-intensity Strombolian activity at Stromboli volcano consists of mild explosions which occur at a rate of several explosions per hour, and is accompanied by nearly continuous active degassing (puffing) Ripepe et al., 2008;Rosi et al., 2013). This activity is fed by crystal-rich, almost degassed magma residing in shallow conduits, a structure seismically identified down to 1 km below sea level (b.s.l.) (Chouet et al., 2008;Patanè et al., 2017). Violent explosions that occurred at Stromboli, called "major" and "paroxysmal" eruptions depending on their size, cover a wide spectrum in terms of recurrence, intensity and magnitude. Such short-lived events involving more than one vent were first defined as Strombolian Paroxysms by Mercalli (1907). Paroxysms have significantly higher intensities (>10 6 kg/s) than major explosions (10 4 kg/s) and differ in the characteristics (crystal content, vesicularity, proportions between recycled and juvenile crystals) of the erupted tephra (Bertagnini et al., 2008). These observations suggest that the two eruption types are fed from different parts of the magma column and involve variable gas/melt ratios. Major explosions typically affect the upper part of the edifice and are unable to reach the inhabited areas of the island. In contrast, paroxysms generate eruptive columns several km in height, ballistics ejecta and pyroclastic flows. Along with tsunamigenic landslides, paroxysms represent a major hazard at Stromboli volcano, posing a serious threat to inhabitants, tourists and scientists (see Bertagnini et al., 2008 andRosi et al., 2013 for a review). These highly energetic events involve the deepest magma ponding zones at the crust-mantle boundary (e.g., Allard, 2010;Métrich et al., 2010;Aiuppa et al., 2011). Preferential pathways allowing magma ascent from deep crustal levels along the main regional NE-SW fault system, and multistep magma ponding zones were proposed by Mattia et al. (2008). The recent records of volcano-tectonic events have allowed the identification of two seismogenic sectors located at 3.5-6 km and 10-12 km b.s.l. (Bonaccorso et al., 2012;Gambino and Scaltrito, 2018). The latter depth coincides with that of the deep magma ponding zone hypothesized on the basis of geochemistry data (Bertagnini et al., 2008;Aiuppa et al., 2010a;Métrich et al., 2010).
The possible mechanisms causing very short-lived paroxysmal eruptions at Stromboli have been widely discussed in the literature. It has been suggested that the intrusion of magma batches with high amounts of dissolved volatile species (H 2 O, CO 2 , S, Cl) and the formation of bubbly basalt blobs able to rise under closed-system conditions (e.g., Métrich et al., 2001Métrich et al., , 2005Bertagnini et al., 2003Bertagnini et al., , 2008Francalanci et al., 2004) led to the pressurization of the crustal magma storage zone (at ∼6-9 km depth b.s.l.). Other authors have assigned a dominant role to the CO 2 -rich gas phase considered to be the driving force of paroxysms. It either involves the rapid transfer of CO 2 -rich gas slugs from the sub-volcano plumbing system (Allard, 2010;Aiuppa et al., 2011) or is related to explosive degassing and fragmentation of CO 2 -oversaturated magma batches due to differential diffusion and disequilibrium degassing of CO 2 relative to H 2 O at high magma ascent rates of 3 to ∼1 m/s (Pichavant et al., 2009(Pichavant et al., , 2013Le Gall and Pichavant, 2016).
Alternative mechanisms discussed in the literature are magma decompression induced by the opening of new pathways in the deep portion of the volcano (Rittmann, 1931), depressurization of the shallow system following conduit obstruction (Falsaperla and Spampinato, 2003). More recently, it was proposed that morphological and structural changes in the summit crater zone, as observed after the 2007 eruption, would allow the ponding and degassing of magma batches in the uppermost conduits (Calvari et al., 2014). These conditions would favor magma mingling and degassing under open-system conditions, decreasing the likelihood of paroxysms (e.g., Calvari et al., 2006Calvari et al., , 2010. Two out of four paroxysms (2003 and 2007) in the last decades occurred during effusive crises from lateral vents along the volcano's flank. Some authors therefore suggested that paroxysms were possibly triggered by the emptying of a critical volume from the shallow system and the downward propagation of the decompression wave (Aiuppa et al., 2010b;Calvari et al., 2011). Based on estimates of magmastatic decompression and decompression rates obtained for all the post-1985 effusive events, Ripepe et al. (2015), Ripepe et al. (2017) proposed that, in the case of lava flows, the decompression rate induced in the deep crustal system by the drainage of the shallow reservoir is mostly controlled by the location of the effusive vents, which regulates effusion rate trends. According to Ripepe et al. (2015), Ripepe et al. (2017), the effusion rate (i.e., the drainage velocity of the shallow reservoir) may control the equilibrium among the different reservoirs in terms of the system's decompression rate more than the total volume of erupted lava during effusive crises.
The triggering mechanisms of paroxysms at Stromboli, to what extent the decompression of the shallowest part of the system propagates downward and how the quantity of deep, gas-richer magma may affect external dynamics (i.e., magnitude/intensity) are highly relevant in terms of volcanic hazard. The aim of the present work is to investigate possible mechanisms driving the paroxysmal eruptions of varying magnitude at Stromboli in the light of new data on samples emplaced during small (2019 CE) and large (1456 CE and1930 CE) paroxysms. We present: (i) new geochemical data for the 2019 paroxysms combined with literature data; (ii) systematic analyses of mineral phases (olivine and clinopyroxene) and Fe-Mg diffusion profiles acquired in olivine in order to discern the timescale of deep magmatic processes and final magma ascent; and (iii) the physical parameters of the eruptions (intensity, magnitude). The compiled geochemical dataset for paroxysms at Stromboli is the largest presented to date. Moreover, we review historical information on the characteristics of paroxysmal eruptions and combine it with observations from the last 20 years to infer their "last-mile" dynamics.

MATERIALS AND METHODS
The present work focuses on tephra emplaced during the two most recent paroxysms occurring on 3 July and 28 August 2019  Figure 1 and Supplementary Table S1 for comparison, whereas the chronology and main characteristics of the investigated paroxysms are detailed in the next section. The analyzed pumice clasts were selected carefully and show no mingling structures (indicative of syn-eruptive mingling with the highly crystallized magma filling the upper portions of the volcano conduits), which could have biased the acquisition and interpretation of geochemical data. Bulk pumice clasts ST07-19 and ST8-19 were analyzed for major and trace elements at the Service d'Analyse des Roches et Minéraux (CRPG-CNRS, Nancy, France) using Inductively Coupled Plasma Optical Emission Spectroscopy (Thermo Fisher ICap 6500) and inductively coupled plasma mass spectrometry, respectively (Carignan et al. 2001). Bulk rock analyses were duplicated (Supplementary Table S2A), and analytical precisions are reported in Pompilio et al. (2011). All samples from previous eruptions were analyzed in the same laboratory, allowing full comparison.
Olivine and clinopyroxene (nearly 700 crystals) from samples ST207, PST205, ST07-19 and ST08-19 were hand-picked systematically from the ¼-½ mm grain-size fraction and the ½-1 mm fraction in the 2019 pumice clasts. They were embedded in Epoxy resin (1 fraction 1 mount). Olivine crystals were positioned parallel to the crystallographic face [100], when clearly identified, and polished until the center was at surface in order to measure core-to-rim profiles. Following the visual orientation of the crystals, the crystallographic faces were subsequently recognized on the backscattered images recorded with a Zeiss EVO 304 MA 10 scanning electron microscope (SEM) at the Istituto Nazionale di Geofisica e Vulcanologia (Pisa). A first round of analyses was acquired using an ISIS-Oxford micro-analytical system linked to the SEM-EDS, with at 15 KeV, using a 300 pA current probe and an acquisition time of 100 s. Concentrations were obtained after ZAF correction using natural minerals or pure oxides for calibration. Crystals showing resorption features, as previously described at Stromboli (e.g., Bertagnini et al., 2003), were discarded. Olivine (∼70 crystals) and clinopyroxene (∼20 crystals) chemical profiles were acquired using the SXFive electron microprobe (Camparis, Sorbonne University, France) with a focused beam current of 20-40 nA, and peak counting times of 10-100 s for Ni, Ca and Mn (Supplementary Tables S2B-F). Analytical accuracy and reproducibility on standards for both SXFive and SEM-EDS are reported in Supplementary Table  S2F. A total of 24 Fe-Mg diffusion profiles were selected to calculate timescales using the DIPRA software (Girona and Costa, 2013), at temperatures of 1,150 and 1,190°C and pressures of 250 and 400 MPa for small and large paroxysms, respectively (Pichavant et al., 2009(Pichavant et al., , 2011Métrich et al., 2010).
Tephra volumes were estimated on the basis of field data (thickness or loading per unit area of fallout deposits vs. isopach or isomass area) for the 1456 and 3 July 2019 events ( Figure 1) and using a single exponential decay law according to Pyle (1989). For the 23 August 1998 eruption, the volume was estimated from the average thickness and deposit dispersal (based on maps from Bertagnini et al., 1999). For the 28 August 2019 event, the total mass was estimated from the observed height of the eruptive column (∼6 km; Giordano and De Astis, 2021) using the relationship of Parfitt and Wilson (1999) for transient events, which links the total mass of erupted material (i.e., total amount of heat release) with plume height, as previously discussed in Pistolesi et al. (2011) for the 2007 event. All the new estimates of tephra volumes obtained from tephra mass are calculated with an average deposit density of 900 kg/m 3 measured in the lab from the FIGURE 1 | Shaded relief map of Stromboli showing sampling sites (this work) for tephra emplaced during paroxysms occurring between the 15th century and 2019. Mass loading (empty red circles, in kg/m 2 ) and isomass curves (red lines), and thickness data (empty yellow circles, in cm) and isopachs (yellow curves) for the 3 July 2019 and 1456 events, respectively, are also shown and were used to estimate volumes of fallout deposits.
Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 9 | Article 593339 volume occupied by a known weight of bulk sample; these are summarized in Supplementary Table S3 and compared with the literature data (e.g., 1930, 2003, 2007 and 2009 events).

Magnitude and Intensity
Paroxysmal eruptions at Stromboli are powerful impulsive explosions lasting a few minutes; they eject pyroclasts up to a few kilometers from the vents and simultaneously produce a convective ash-gas plume of up to ten km, ending with minutes of vigorous degassing and ash emission (Ponte, 1919;Rittmann, 1931;Barberi et al., 1993;Rosi et al., 2006;Pistolesi et al., 2011). The first phase is usually accompanied by the emission of metersized, ballistic scoria bombs and blocks, which are dispersed up to 2 km from the craters along the NE-SW direction of alignment of the feeder dyke and vents. The non-symmetrical, west-southwest and northeast distribution of the blocks suggests that the ejection direction is generally controlled by the geometry of the shallow conduit system feeding the summit conduits. Ballistic blocks are ejected with velocities of up to 200 m/s and likely represent the brittle failure of the conduit wall rocks and of the partially crystallized shallow plumbing system due to the rise of the over-pressurized magma slug (Rittmann, 1931;Rosi et al., 2006;Renzulli et al., 2009;Pistolesi et al., 2011). Capaldi et al. (1978), Barberi et al. (1993), Rosi et al. (2013) and Bevilacqua et al. (2020) have compiled a detailed review of all recent events for which information is available. Barberi et al. (1993) and Rosi et al. (2013) point out that, in addition to the emission of large quantities of blocks and ash, practically all the historical events have also involved the ejection of highly vesicular pumice varying widely in density and with a variable degree of mingling between shallow and deep magmas. Studies on the 2003 and 2007 paroxysms indicate that the bulk volume of tephra fallout deposits was of the order of 10 4 -10 5 m 3 (10 7 -10 8 kg in mass), the typical peak mass discharge rate was 10 6 -10 7 kg/s and gas-ash jet velocities exceeded 200 m/s during the first seconds of the event following the system's rapid decompression (Rosi et al., 2013 and references therein).
If one considers a wider interval of time, the range of erupted volumes increases, suggesting a higher variability in eruption magnitude/intensity. We highlight that the 2003 and 2007 paroxysms (like those in 2019) were smaller (Rosi et al., 2006, Pistolesi et al., 2011 this work) than the larger events occurring in previous centuries (e.g., 1930 and 1456 CE), for which the estimated magma volume was at least one order of magnitude greater this work). There are also other paroxysms (15th-17th centuries or undated) which, based on tephra thickness and dispersal, can be considered large in scale Métrich et al., 2010;Pichavant et al., 2011) and shall serve for comparison in the next section.

Inferred Source Depth of Paroxysms
The compiled published data on melt inclusions reveal that sulfur degassing of magma batches that fed small paroxysms followed the same path ( Figure 2). Sulfur starts degassing at a pressure of 150 MPa, according to sulfur solubility experiments conducted on Stromboli magma (Lesne et al., 2011). This would place the depth of magma decompression at nearly 5 km b.s.l. The depth of the ponding source magma prior to the small 15 March 2007 paroxysm is reported at 150-250 MPa on the basis of experiments (e.g., Pichavant et al., 2009Pichavant et al., , 2011, and at 210 ± 20 MPa (∼7 km b.s.l.) on the basis of the total fluid pressure (P H2O +P CO2 ) recorded by the olivine-hosted melt inclusions (Métrich et al., 2010). A deeper origin (up to 15 km) has been also proposed and discussed for large paroxysms (e.g., Bertagnini et al., 2008;Aiuppa et al., 2010a;Allard 2010;Métrich et al., 2010).
Notwithstanding a possible deeper source, the final (∼7 km b.s.l.) ascent of the magma batch under closed-system decompression has been considered fast enough to inhibit crystal nucleation (e.g., Bertagnini et al., 2008;Allard 2010;Aiuppa et al., 2010b;Métrich et al., 2010), in agreement with the increasing acceleration of the gas-melt mixture upon decompression and final fragmentation (e.g., Parfitt and Wilson, 1999). This interpretation is in very good agreement with: (i) the high vesicularity of the pumice (∼80%; Pioli et al., 2014) and, (ii) the volcano-tectonic events that preceded the 2007 eruption. The latter were localized at a depth between 3.5 and 6 km b.s.l. and interpreted as a possible brittle response of the crust to magma pressurization and ascent from deeper zones (Gambino and Scaltrito, 2018).

Targeted Eruptions of Variable Magnitude
The four targeted paroxysms of variable magnitude occurred in 2019, 1930 and 1456 CE. During summer 2019, Stromboli  Bertagnini et al. (2003). The pressure of 150 MPa (gray area) at which sulfur is exsolving in the gas phase was determined experimentally for Stromboli magmas (Lesne et al., 2011). Abbreviations are: LP and SP for large-and small-scale paroxysms, respectively, and ME for major explosions.
Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 9 | Article 593339 entered a phase of intense explosivity, with two paroxysms separated by two months of intense, violent explosive activity punctuated by lava overflows from the summit craters. The first paroxysm occurred in the early afternoon of 3 July 2019, with the formation of an explosive plume almost 8 km high (Giordano and De Astis, 2021). The violent explosion covered the village of Ginostra with ash, lapilli and bombs that immediately ignited fires when the hot pyroclastic fragments reached the dry vegetation ( Figure 3). The partial collapse of the convective column generated two pyroclastic flows which, impacting the sea surface, produced a two m-high tsunami that reached the cost of southern Italy half an hour later (Fornaciai et al., 2019). The formation of density currents during paroxysms has been typically associated with primary, directed ejection of pyroclastic material and/or partial collapse of the vertical jet and usually remains confined within Sciara del Fuoco, as already observed for the small 2003 event (Rosi et al., 2006). Nobody was directly injured during the 3 July explosion, but one tourist eventually died after trying to escape the rain of hot lapilli and bombs that hit the western sector of the island. Observations around the crater area revealed a high thickness of spatter bombs and lapilli (as high as 1 m west of the crater area); mass loading carried out along the west coast allowed a volume estimate of 3.5 × 10 5 m 3 for the fallout phase. After the 3 July paroxysm, Strombolian activity remained very high (i.e., number and intensity of the Strombolian explosions), and a phase of continuous lava outpouring from the summit craters started. On 28 August 2019, a second paroxysm occurred, which covered the summit areas with bombs, lapilli and minor lithic blocks and the village of Stromboli with several centimeters of ash, mostly reproducing the sequence of events registered on 3 July. The column was slightly lower, with a fallout volume of  Table S3); again, a pyroclastic flow along Sciara del Fuoco triggered a tsunami. Effusive activity ceased definitively after a few days, on 31 August, although higher degassing from the crater area persisted for several weeks after the eruption ended. The 1930 event, very likely the most energetic of the post-1900 paroxysms, is the eruption that produced the highest number of casualties (at least eight). It has been studied in great detail by Rittmann (1931), who gathered residents' accounts, observed the deposits and conducted a topographic survey of changes in the crater produced by the explosion. The peak phase of the eruption was preceded by the ejection of lithic materials and lasted about 10 min; it generated a column several km high and produced m-sized spatter bombs, dm-sized scoria, lapilli and ash with a minimum volume of 2.5 × 10 6 m 3 ) and a tsunami of uncertain origin. The high accumulation rate of pyroclastic material along the cone's slopes in 1930 resulted in the formation of pyroclastic density currents that developed outside the Sciara, as also reported for other large paroxysms (i.e., 1906, 1944Riccò, 1907;Rittmann, 1931;Ponte, 1948).
The paroxysmal eruption in 1456 CE was selected because it represents one of the largest events occurring in historical times, as suggested by field observations and by the dispersal and thickness of tephra deposits (Rosi et al., 2019;Pistolesi et al. 2020; this work). The fallout deposit occurring in several outcrops on the eastern sector of the island was used to roughly estimate a tephra volume of the order of 5 × 10 6 m 3 . In addition to its high magnitude, it has been suggested that this event was linked to a tsunami-triggering flank collapse of Sciara del Fuoco. The stratigraphic sequences cropping out on the NE sector of the island reveal a tsunamiite that was immediately covered by the paroxysm deposits (Rosi et al., 2019;Pistolesi et al., 2020).
As a whole, the magnitude/intensity of explosive activity at Stromboli varies considerably. Both the erupted tephra mass and mass flux vary by several orders of magnitude and are perfectly scaled to the height of the ejected products ( Figure 4). In particular, the scale of paroxysms differs by two orders of magnitude, with both small-(e.g., 2003, 2007 and 2019 CE events) and large-scale paroxysms (e.g., 1456 and 1930 CE).

Geochemical Characteristics of the Studied Pumice Clasts
The bulk compositions of the 3 July and 28 August 2019 samples ( Supplementary Table S2A), with 6.3-6.4 wt% MgO, plot in the domain of Stromboli pumice clasts ( Figure 5). As a whole, CaO and MgO concentrations vary widely but without distinction between small and large paroxysms ( Figure 5), indicating that the magmas that supplied these eruptions reached comparable degrees of differentiation.  Ripepe et al. (1993) and Patrick et al. (2007); major explosions refer to the 2009 events (small green dots; Andronico and Pistolesi, 2010;Pioli et al., 2014) and to their average [large green dot, also used in (B)]. Paroxysms refer to 1919,1930,2003,2007 and 2019 events (Ponte, 1919;Rosi et al., 2006;Bertagnini et al., 2011;Pistolesi et al., 2011;Pioli et al., 2014;Giordano and De Astis, 2021;this work). In the case of 1919, 1930 and August 2019, only one parameter was known, with the other obtained by using the relation of Parfitt and Wilson (1999). They are displayed anyway to show the linear correlation within the whole dataset. In   (2011). Sample PST9 was used as the starting material for high pressure experiments (Pichavant et al., 2009; According to experimental works (Pichavant et al., 2009(Pichavant et al., , 2011, the 2019 samples would plot in the stability domain of clinopyroxene. However, high-pressure experiments were performed on one pumice sample (PST9) which is relatively enriched in MgO with respect to the others, possibly shifting the cotectic line to the right side of the diagram. We therefore suggest that the 2019 samples plot on the olivine-clinopyroxene cotectic line, as do samples ST207 and PST205 and most bulk rocks of present-day activity at Stromboli. This conclusion is in agreement with the MELTS simulations of magma crystallization paths at Stromboli (Métrich et al., 2010).
There is an overall positive correlation between La and Nb ( Figure 6A), and among incompatible elements in general. In contrast to the latter, the concentrations of elements such as Co and Yb vary within very limited ranges ( Figure 6B). This implies that trace elements track processes in the mantle source, where the removal of olivine and clinopyroxene (which should have depleted the residual melt in compatible elements) has little impact. These results support previous conclusions reporting magma source effects (source enrichment, degrees of melting) on Stromboli magma geochemistry, with variations in Rb/Th, Ba/ La and Nb/Ce ratios and Sr isotopes over the last 1800 years (e.g., Francalanci et al., 2004;Tommasini et al., 2007;Pompilio et al., 2011).
In sample from the 1456 CE paroxysm (ST207), olivine crystals display Fo 89 -Fo 90.5 cores and relatively thick Fo 85 -Fo 86 rims. They show normal zoning and systematic . Data are comparable as they were provided by the same laboratory; we also verified reproducibility, as reported in Pompilio et al. (2011). Symbols as in Figure 5.
Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 9 | Article 593339 swallow tails, indicating ongoing crystallization ( Figure 7A). These olivine crystals represent nearly 15% of the population in a set of 100 hand-picked crystals and contain no or very rare glass inclusions but tiny Cr-spinels ( Figure 7A). They coexist with olivine crystals with hopper-like texture ( Figure 7B), which display compositions in the range Fo 85.8-87.9 and two-step crystallization. Olivine starts crystallizing (Fo 87.9 ) as internally hollow crystals (hopper olivine) and the central cavity, containing a wealth of melt inclusions, is subsequently sealed as Fo 85.8 . These features are also described for other paroxysms (e.g., undated ST82 paroxysm and 2007 CE event; Métrich et al., 2010). They result from cooling-heating cycles during crystal growth (Faure and Schiano, 2004) and thermal disequilibrium events associated with magma mixing or magma convection (Mourey and Shea, 2019 Figure 7C). We found rare crystals with a swallow tail and texture indicating ongoing crystallization ( Figure 7D). Slightly reversely zoned clinopyroxene (Mg# 0.72-0.76; Fs 16.5-13.7 ) is found in association with olivine Fo 80-84 . Clinopyroxene microcrystals (Mg# 0.70, Fs 18 up to Mg# 0.88, Fs 6.5 ) with surface dissolution textures, highlighted by the green-light yellow color, are ubiquitous and inherited from the surrounding.
The new mineralogical data for pumice samples of the large 1456 CE (ST207) and 1930 CE (PST205) paroxysms fully confirm the presence of Fo 90 magnesian olivine, containing Cr-spinel inclusions, with diopsidic clinopyroxene, which represent the liquidus phase of deep-seated magma at Stromboli, in accordance with data from the literature (e.g., Bertagnini et al., 2008). Moreover, further mixing between this CO 2 -rich magma batch and the overlying and resident magma could be regarded as a general process at Stromboli as deduced from our present dataset and previous observations (Bertagnini et al., 2008;Pichavant et al., 2009;Métrich et al., 2010).
The small 2019 paroxysms provide new and interesting information. In July 2019 pumice clasts (ST07-19), olivine crystals from the ¼-½ mm (microphenocrysts) and ½-1 mm fractions (phenocrysts) display the same limited range of compositions, from Fo 86.7 to Fo 83.5 and Fo 86.4 to Fo 83.6 , respectively (Supplementary Figure S1). However, among a total of 25 analysed crystals, the composition of 60% of the phenocrysts ranges from Fo 85.7 to Fo 86.3 and that of 61% of the microphenocrysts between Fo 85.4 and Fo 86.0 . In both Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 9 | Article 593339 8 fractions, olivine commonly shows a central portion richer in melt inclusions and empty gas bubbles, with no regular shape, surrounded by a relatively thick rim ( Figure 7E): a texture similar to that of olivine in 1456 CE sample. As for clinopyroxene, Mg# varies from 0.87 to 0.84 (7.3-8.3 Fs) in phenocrysts and from 0.88 to 0.84 (6.3-8.2 Fs) in microphenocrysts. For both olivine and clinopyroxene, compositional ranges for the two grain-size fractions are therefore comparable. This points to similar P-T crystallization conditions. The Fe-Mg partition coefficient (K d olivine/melt) between the olivine and melt in the 2019 pumice samples calculated at ΔNNO+1 ( Table S2G).
For each sample, Fe-Mg diffusion profiles (Supplementary Table S2H) calculated using the DIPRA model (Girona and Costa, 2013) were used to estimate the timescales of magmatic processes. For large-scale paroxysms, timescales range between 56 and 2 days in 1456 CE sample (ST207) and between 58 and 8 days in 1930 CE sample (PST205). For small-scale paroxysms, the interval is longer, from 110 to 5 days in July 2019 and 115-10 days in August 2019. In 1930 CE sample, several profiles were unusable due to melt inclusions. Moreover, in two crystals (e.g., Ol56 in Figure 7C), profiles measured along the c and b axes provided different timescales ranging from 1 to 2 months along the c axis and 5-8 days along the b axis. This feature could be related to the orientation of the olivine and/or to the crystallization dynamics. Mourey and Shea (2019) suggested that olivine growth along the c-axis could be favored under lower undercooling conditions. However, brief (i.e., days) timescales are in agreement with rapid-growth morphologies such as dendritic overgrowths ( Figure 7D), indicating fast rates during late-stage crystallization (Faure and Schiano, 2004). Even more importantly, all samples show final crystallization after 2-10 days, an indication of rapid ascent rates, as discussed in the following section.

Is there a relationship between magma geochemistry and eruption dynamics?
By looking at the frequency record of paroxysms occurring in the last 120 years, we observed that the frequency during the 20th century was one event every 3-4 years in the first 60 years, whereas the first paroxysm after the 1959 event occurred in 2003.
Since this year, there have been four events in 15 years Bevilacqua et al., 2020).
Trace element data reveal that magma batches extruded from 1998 to 2019 were characterized by relatively high average Rb/Th (4.7 ± 0.2) and Ba/Th (70 ± 2) ratios but low Ta/Th (11.6 ± 0.2), like the 1456 CE samples. These results support previous conclusions on the FME-enriched signature of the most recent magmas at Stromboli . The available dataset, which is not exhaustive because some tephra deposits have not been dated, also reveals that the intense paroxysmal activity of the 16th and 17th centuries produced 1930-type magmas.
According to Tommasini et al. (2007), the mantle source underlying Stromboli experienced two main periods of metasomatism, as deduced from Th-U-Sr-Nd-Pb isotopes. The first was associated with supercritical fluids derived from the subducting Ionian slab at high pressures (>5 GPa), the second with aqueous fluids enriched in fluid-mobile elements (FME; such as Ba, Rb, K) that were released at relatively lower pressures (<5 GPa). Mixtures of these fluids could explain the wide compositional range of magmas erupted over Stromboli's lifetime. From the geochemical dataset reported here, the magma batches refilling the plumbing system in the past ∼120 years record short-lived variations. The latter are of low amplitude but significant and may reflect slightly variable contributions from the "uppermost" FME-enriched mantle source.
The main point here is that there is no relationship between the trace element composition of magmas at Stromboli and eruptive dynamics, since small (i.e., 2003, 2007, 2019) and large paroxysms (i.e., 1456 CE), or even major explosions and August 2019 pumice and reported in previous works (e.g., Bertagnini et al., 2003;Francalanci et al., 2004;Pompilio et al., 2011), is the dominant process governing the long-lasting steady state at Stromboli.
FIGURE 9 | Timescales representative of magmatic processes calculated from Fe-Mg diffusion modeling compared to Fo (mol%) concentration profiles in olivine crystals from samples ST07-19 (July 2019) and ST207 (1456). Modeling was performed using the DIPRA software (Girona and Costa, 2013). Natural data (black) with error bars and calculated diffusion curves for minimum (blue) and maximum time (green) are shown. Best fit to Fo natural data is represented in red. The crystallographic axes (in square brackets) of olivine used for crystal orientation are also indicated.

Crystal Records and Magmatic Timescales
The log-log distribution of olivine compositions (Figure 8, reporting comparable systematic data only) can be used to distinguish paroxysms of differing magnitude.
Data from 1456 to 1930 CE samples (ST207 and PST205, respectively) confirm that olivine Fo 90, with Cr-spinel, is present as euhedral microcrystals. In the former, Mg-rich cores with swallow tails rims (∼Fo 86 ) are ubiquitous and suggest that olivine grew in equilibrium with a melt of similar composition to the bulk pumice (or matrix glass). The normal Fe-Mg zoning of the ST207 olivine, attributed to mixing, is well preserved (minimum Fe-Mg diffusion) because of high decompression rates and quenching. These observations are also in agreement with the presence of pure glassy matrices in the ST207 pumice clasts. A second crystal population consists of irregularly shaped olivine (Fo 87.9-85.8 ) containing numerous melt inclusions and empty gas bubbles at the time of entrapment (e.g., Figure 7B). Both the episode of rapid growth, favoring melt and gas entrapment, and the presence of sector zoning in clinopyroxene testify to a rapid crystallization event.
In the case of the 1930 CE eruption (PST205), olivine compositions range widely and crystals show resorption features and reverse zoning, features that are consistent with magma ascent through fractures in crystal mush (es). Mg-rich microcrystals with swallow tails ( Figure 7D) exist but are rare.
The mineralogy of the pumice emplaced in 1456 CE and 1930 CE, together with the compiled published data on paroxysmal explosions, suggests a scenario in which paroxysmal explosions are controlled by the infracrustal magma ponding zone. Deepderived magma batches in equilibrium with Fo 90 olivine (and diopsidic clinopyroxene) mix with the resident magma ponding at 6-9 km and crystallize rapidly, reaching equilibrium conditions before erupting explosively and quenching.
The timescales recorded by Fe-Mg diffusion profiles (Supplementary Table S2H; Figure 9) indicate fast ascent rates of 2-8 days from 6 to 9 km to the surface for the 1456 CE paroxysms (ST207). We consider these estimates to be reliable, as they were measured in five crystals (out of nine). Timescales for the 1930 CE paroxysm are fully comparable (5-8 days). In both cases, we also estimated longer times (i.e., 2 months), which may represent the interval between refilling of the reservoir at 6-9 km and final magma ascent.
In paroxysms of smaller magnitude such as the July 2019 event, the olivine composition is predominantly Fo 85.4-86 , as also reported for the 15 March 2007 paroxysm (Fo 85.5-86.5 ;Métrich et al., 2010). Olivine-melt equilibrium in both phenocrysts and microcrysts of the July 2019 samples indicates remarkably stable crystallization conditions through time. In other words, the magma always ponds at the same level of the feeding system. In July 2019, steady-state conditions prevailed prior to the paroxysm, with no clear evidence of significant arrival of deeper-derived magma (melt). This strongly suggests that the system was prevalently destabilized by intrusion of a CO 2 -rich gas phase, which induced active convection. In contrast, the chemical spectrum of the 28 August 2019 samples is wider, and reverse zoning in olivine testifies to the arrival of a new magma batch. Given that the July 2019 paroxysm was preceded by a major explosion occurring on 25 June 2019, we propose the following sequence of events: (i) pressurization of the 6-9 km system by deep refilling of magma and its CO 2 -rich gas phase characterized by a high gas/melt ratio; (ii) major explosion, accompanied by discharge of a CO 2 -rich gas phase; (iii) gas-induced destabilization of the system, magma decompression and a small paroxysm on 3 July (8 days after the major explosion); (iv) a second small paroxysm on 28 August, with evidence of magma mixing in the erupted products. For the 28 August paroxysm, a contribution to the decompression of the feeding system due to the continuous lava overflows during July and August cannot be ruled out.
Equilibrium crystallization may have started up to 4 months (110 days) prior to the July paroxysm, as testified by euhedral phenocrysts in pumice clasts, but it took no more than 5-10 days prior to explosion to complete the last stage of crystallization and decompression. We stress that the ascent rate of 6-10 days for the July crystals agrees with the interval between the major explosion on 25 June and the paroxysm on 3 July.
In summary, Mg-Fe diffusion profiles in olivine crystals suggest that the timescale of final magma ascent (i.e., days) prior to eruption is remarkably comparable among large-and small-scale paroxysms.

Eruption magnitude and crystal cargo: what are the magma volumes?
Our dataset indicates that olivine compositions can be used as a good proxy for the volume of erupted magma. We show that the volume of magma emitted during known eruptions at Stromboli in the last 1300 years of activity increases exponentially from major eruptions to small and large paroxysms (Figure 4). The volatile-rich magma that feeds the Stromboli system contains a small amount of crystals (nearly 3 vol%) and commonly mingles with a crystal-rich, shallow-seated magma during major eruptions and paroxysms (Pioli et al., 2014). Olivine is thus rare but ubiquitous, and can be studied as a proxy for the magma involved in explosive eruptions. The amount of volatilerich magma in violent explosions can be derived from total tephra volumes and componentry analysis of deposits (Supplementary Table S3). We thus consider only the deep magma contribution for paroxysms and major eruptions Pioli et al., 2014).
Olivine Fo 89-90 is systematically associated with the largest tephra discharge (≥10 6 m 3 ; Figure 10). The lower volume (10 4 -10 5 m 3 ) of tephra emplaced during small-scale paroxysms correlates with the lower Fo content of the associated olivine crystals. This trend is further confirmed by the progressive decrease in erupted tephra and in Mg contents in olivine (Fo mol%) from paroxysms to ordinary activity, where the volatile-rich component is absent. Lastly, based on the vesicularity of the tephra and the average 80% porous volume, we estimate that the amount of melt feeding the eruptions is roughly 10 1 -10 5 m 3 . This implies that the volumes of gas-rich melts at play can be small.

Decompression Events Prior to Paroxysms
One of the questions that arises from the investigation of historical paroxysms is the relative role of deep magma  (Ripepe et al., 2017). In most cases, paroxysms occur suddenly without any precursor or are in some cases heralded by seismic events days before, as reported in 1916 (Ponte, 1921). Although the 1456 CE paroxysm was tentatively linked to a lateral, tsunamigenic flank collapse (Pistolesi et al., 2020), this represents the only known case of paroxysm preceded by a lateral collapse. We therefore suggest that pressurization of the deep magmatic system is the usual trigger of paroxysms. In this view, flank collapse and/or lava flows may enhance the decompression dynamics of the eruptible, gas-rich magma in the deep reservoir.
Paroxysms thus testify to a complex change in magma recharge rate and to the likely opening of new pathways in the volcanic structure, enabling the direct, rapid ascent of magma. Given the complexity of the triggering mechanism, interpretation of the degassing history, particularly that of CO 2 , is not straightforward; nevertheless, it has been proved to be a precursor of rapid magma ascent and used as an effective tool for monitoring activities, even major explosions FIGURE 10 | Relationship between volume of erupted tephra and olivine composition expressed as Forsterite (mol%) [100×Mg/(Mg + Fe)] (A) tephra volumes are reported for an ordinary explosion (typical Strombolian activity; Rosi et al., 2013); for major eruptions and paroxysms we considered only the deep magma contribution based on the ratio between volatile-rich and volatile-poor-related tephra (for high-porphyritic/low porphyritic ratios see Pioli et al. (2014events, and Bertagnini et al. (1999 and Bertagnini et al. (2011Bertagnini et al. ( ) for 1998Bertagnini et al. ( and 1930. For 1456 and 2019 events we used an average ratio of 0.7 based on observed deposits and tephra. (B) Only major eruptions and paroxysms are considered here. The gray area represents the transitional domain between large-scale paroxysms (closed-system degassing and magma ascent) and ordinary explosive activity (open-system degassing).
Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 9 | Article 593339 (Aiuppa et al., 2010b(Aiuppa et al., , 2011. Published data on melt inclusions evidence a cluster in the CO 2 -H 2 O diagram in the domain of 250-300 MPa, as reviewed by Pichavant et al. (2018). The highest dissolved amounts of CO 2 and H 2 O recorded by melt inclusions may reach 1500-1600 ppm and 3 wt%, respectively (e.g., Métrich et al., 2010). Moreover, it is well demonstrated now that a large proportion (≥80%) of the CO 2 migrated into the contraction bubble associated with the melt inclusion (e.g., Wallace et al., 2015;Moore et al., 2015Moore et al., , 2018. As expected, a significant proportion of CO 2 is still present in the gas phase at ≤500 MPa in basaltic systems (e.g., Wallace, 2005;Lesne et al., 2011). In contrast, the total amount of CO 2 in the system at the time of magnesian olivine crystallization could have reached ≥0.8-1 wt%, a value lower than that estimated for the gas fluxes of parental magma at Stromboli (∼2 wt%; Allard, 2010). This suggests that at the time of mixing between the deep-derived magma and the magma residing at 6-9 km, a large proportion of CO 2 had already exsolved as a gas phase, a good trigger for paroxysms. If this holds true for the large paroxysms, a larger proportion of CO 2 must have been injected within the system together with the melt. Note that although oversaturated CO 2 magma increases the decompression rate (Pichavant et al., 2018), Stromboli melt still contains a high proportion of H 2 O (3 wt%) at the time of decompression. Given the high diffusion rate of H 2 O and the increase in its molar volume with pressure with respect to CO 2 , the effect of water during the decompression event cannot be neglected. The amount of H 2 O dissolved in Stromboli magma reaches ∼2.7-3.5 wt% prior to outgassing (i.e., Métrich et al., 2010). Upon magma decompression and ascent, plagioclase starts to crystallize as H 2 O is exsolved in the gas phase and becomes the dominant phase, along with clinopyroxene, as the H 2 O content of the melt reaches ∼1.2 wt% (e.g., Di Carlo et al., 2006). Slow magma ascent and ponding in the shallow system would promote water degassing and plagioclase crystallization. Plagioclase is scarce in Stromboli pumice, and observed specimens result from very late-stage degassing (e.g., Bertagnini et al., 2008). These observations contradict earlier hypotheses (e.g., Calvari et al., 2005), suggesting a low-pressure ponding and degassing of the deep-derived magma before paroxysm. Irrespective of the deep trigger, the rapid ascent of a CO 2 /H 2 O-rich slug enveloped by melt always involves a final fragmentation event at very shallow levels. As observed for the 2003, 2007 and 2019 events, fragmentation in the conduit is eventually followed by the downward propagation of the fragmentation level, resulting in an initial spherical tephra ejection which progressively shift toward a vertical jet. As for the recurring major explosions at Stromboli, the CO 2 -rich gas was considered to be the trigger (Aiuppa et al., 2011). Since the compositions of the glassy matrices do not vary significantly from large paroxysms to major eruptions (Supplementary Table  S4), the higher crystallinity reflects the high proportion of inherited crystals. The latter are entrained from the volcanic pile and from the overlying magma that is strongly degassed and rich in large euhedral crystals of plagioclase, olivine and clinopyroxene.

Response of the Shallow System to the Ascent of Deep Magma
The feeding system at Stromboli is traditionally considered to be in steady-state equilibrium, but temporal variations in the eruptive regime are recorded during the continuous acquisition of geophysical and geochemical parameters (e.g., Aiuppa et al., 2010a;Ripepe et al., 2017). These 'cyclic' variations in volcanic activity occur over the mid-term (months to years). During periods of low-level activity, the magma level tends to drop in the conduits, and Strombolian explosions become weaker and fewer (i.e., <10 per hour; Ripepe et al., 2008). More intense activity at the craters is typified by a rise in magma level, in the number (up to 20 per hour) and energy of explosions, the height of gas jets and their load of incandescent fragments (Ripepe et al., 2017). Lava overflowing at the summit craters and/or dyke opening along Sciara del Fuoco are heralded by an increase in activity, as repeatedly reported in the last decades (Calvari et al., 2005(Calvari et al., , 2010Lodato et al., 2007;Valade et al., 2016;Ripepe et al., 2007Ripepe et al., , 2017, and only in some cases have preceded paroxysmal events. When paroxysms occur, differences in external dynamics may be related to the different state of the shallow conduit prior to eruption. For example, lava effusion can drain the magma in the conduit and induce final fragmentation at very shallow levels. When paroxysms occur during lava outpouring, once the conduit is partially emptied by lateral magma drainage, the fragmentation depth is limited by the magma level and deepens. This results in external jets of pyroclasts that travel at length along the conduit and are eventually constrained by the conduit walls. The deepening of the final fragmentation phase also results in the emission of light gray, crystal-rich to holocrystalline, ballistic subvolcanic rocks representing the deeper levels of the conduit that are broken through rapid decompression, as observed during the 2003 and 2007 events, when lava effusion was in progress and the conduit was partially emptied (Rosi et al., 2006;Pistolesi et al., 2011;Del Moro et al., 2013;Ripepe et al., 2017). Ballistic ejection of crystalrich subvolcanic rocks was also reported during the 1919 and 1930 paroxysms (Ponte, 1919;Rittmann, 1931), although no lava emission preceded the events. In these particular cases, because of the higher volume of gas-rich, deep-derived magma, the process of fragmentation may have involved a longer (hundreds of meters) portion of the conduit, affecting the section at depth from which rocks forming the crater walls of the volcano are usually ripped.
When the conduit is full of magma (no ongoing lava effusion), fragmentation is not limited by the magma level and becomes shallower, developing a giant gas bubble that explodes at the surface. In the latter cases, as observed in July and August 2019, no fresh lithic blocks are ejected because fragmentation takes place in a region of high-temperature, ductile rocks, and crystalrich magma residing in the conduit is forced upward by the ascending gas-rich magma. As a result, the shallowest magma has to be ejected first; after the 2019 events, crystal-rich scoriae were Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 9 | Article 593339 found to be covered by (and thus ejected before) the late, highvesicular spatters and lapilli forming the eruptive column. When the conduit is not emptied by lava flows, the ratio among crystalrich, shallow magma and deep-derived magma is higher with respect to paroxysms that occur during effusive crises, which contribute to the discharge of shallow magma prior to the paroxysm.

CONCLUDING REMARKS
Historical paroxysms highlight a very wide range of conditions in the shallow part of the volcano. While some occurred during lava effusion or were heralded by a lateral collapse, the majority (including those of 2019) occurred during the ordinary Strombolian regime, with no significant precursory signs in the long term. This variability suggests that paroxysms are controlled mostly by the deep magmatic system, with the shallow dynamics (e.g., effusive crises, lateral flank collapses) acting as possible catalysts by favoring the ascent of deep magma through decompression of the shallow reservoir. The geochemistry of bulk tephra samples, along with crystal records and physical parameters, allowed us to shed light on the source and dynamics of paroxysms at Stromboli volcano. Although trace elements show no significant correlation with the magnitude or intensity of paroxysmal eruptions, crystal data prove that primitive (Fo 89-90 ) olivine is ubiquitous in the crystal population of large-scale paroxysms and represents the refilling of the 6-9 km-deep reservoir by the deeper magma most likely ponding in the lower crust or at the crust-mantle transition. In contrast, melt-olivine equilibrium testifies to steady-state conditions prior to small-scale paroxysms (e.g., March 2007 and July 2019), suggesting that the 6-9 km reservoir was pressurized and then destabilized by the transfer of a CO 2 -rich gas phase. The associated melt fraction is recorded by the crystal zoning of August 2019. The occurrence of large-or small-scale paroxysms (and major eruptions) can therefore be considered to be determined by the volumes of magma involved and by the gas/ melt ratio (Figure 11).
The high magma ascent rates, calculated from Fe-Mg diffusion profiles for both small-and large-scale paroxysms, range between 2 and 10 days, suggesting that the "last-mile" magma ascent time from 6 to 9 km to the surface is a common feature irrespective of eruption intensity. Longer timescales (months) correspond to crystallization within the 6-9 km b.s.l. reservoir induced by the arrival of melt and its CO 2 -rich gas phase in variable proportions, convection and magma mixing ( Figure 11).
We also show that olivine composition correlates perfectly with eruption intensity, suggesting that the magma source controls the eruption dynamics (magnitude/intensity). The volume of magma involved in the large eruptions is in fact at least one order of magnitude greater than that erupted during FIGURE 11 | Schematic representation of olivine crystallization and growth during magma ascent for (A) small-scale and (B) large-scale paroxysms. For smallscale paroxysms, we propose a process of deep recharge (1-Input) of predominantly CO 2 -rich gas. As an example, in July 2019, olivine (Fo 86 ) shows slight normal zoning (2-Reservoir) suggesting that destabilization is induced only by a CO 2 -rich gas phase (with almost no melt). In contrast, olivine reverse zoning in August 2019 points to a process of magma mixing (2-Reservoir). For large-scale eruptions, the deep-recharge of CO 2 -rich gas, melt and olivine (Fo 90 ) enters the 6-9 km reservoir (1-Input), and induces olivine growth with variable textures (2-Reservoir). In both cases (A) and (B), steady-state conditions are reached prior to eruption and timescales of final magma ascent (3-Ascent), represented by the growth of outer rims (dashed line) on crystals of variable textures, are in the order of 2-10 days. Numbers in italic refer to Fo (mol%).
small-scale paroxysms, and three orders of magnitude greater than that erupted during major explosions, showing a good correlation with the chemistry of ubiquitous olivine crystals. As a whole, the correlation among chemical (e.g., degassing extent, crystal composition) and physical (mass eruption rate) parameters is valid for the overall spectrum of eruptions (from mild Strombolian activity to major and paroxysmal eruptions). This brings strong evidence that the magmatic source ultimately controls eruptive dynamics.

DATA AVAILABILITY STATEMENT
All datasets presented in this study are included in the article as Supplementary Material.

AUTHOR CONTRIBUTIONS
NM performed electron microprobe analyses, treated mineral chemistry and bulk rocks data, compiled data from the literature, and wrote the manuscript with contributions from AB and MP. AB collected 1930 and 1456 pumice samples, separated olivine and clinopyroxene crystals from pumice, took SEM images and helped write the manuscript. MP collected 1456 and 2019 samples, as well as volcanological data in the field, reviewed the physical parameters of paroxysms, performed diffusion profile modeling and helped write the manuscript and prepare figures.

FUNDING
This work was supported by the Italian Civil Protection in the framework of the DEVNET project, under an agreement between the Universities of Florence and Pisa.