Mafic magma feeds degassing unrest at Vulcano Island, Italy

Abstract

The benign fuming activity of dormant volcanoes is punctuated by phases of escalating degassing activity that, on some occasions, ultimately prelude to eruption. However, understanding the drivers of such unrest is complicated by complex interplay between magmatic and hydrothermal processes. Some of the most comprehensively characterised degassing unrest have recently been observed at La Fossa cone on Vulcano Island, but whether or not these episodes involve new, volatile-rich ascending magma remains debated. Here, we use volcanic gas measurements, in combination with melt inclusion information, to propose that heightened sulphur dioxide flux during the intense fall 2021 La Fossa unrest is sourced by degassing of volatile-rich mafic magma. Calculations using a numerical model indicate observations are consistent with the unrest being triggered by the emplacement of ∼3·10⁶ m³ of mafic magma at ∼4–5 km depth. Degassing of mafic magma is argued as a recurrent driver of unrest at dormant volcanoes worldwide.

Introduction

Of the ~1500 volcanoes with known or postulated eruptive activity¹ in the Holocene, only ~10% (or less) are erupting or persistently degassing at present^2,3. The global volcano population is, therefore, dominated by dormant volcanoes that exhibit mild to feeble hydrothermal fuming activity^3,4,5. A major challenge in volcanology is to detect and interpret the reactivation signs of these dormant volcanoes, many of which are close to densely populated regions¹. However, the periods of volcanic unrest that precede volcano reawakening are often associated with complex, non-linear patterns in volcanological, geophysical, and geochemical monitoring parameters^6,7,8,9,10, posing challenges for the timely recognition of eruption precursors^11,12.

One recurrent sign of dormant volcanoes becoming restless is an escalation in surface degassing activity^13,14, typically manifesting as accelerating outgassing rates^2,15, increasing fumarole temperatures and/or extension of degassing areas^16,17, and increasingly magmatic (and less hydrothermal) chemical characteristics of the emitted gases^{18,19,20,21,22,23}. As not all degassing unrests ultimately culminate into eruption⁸, the fundamental challenge is to distinguish, from data streamed by monitoring networks, between pre-eruption magmatic unrests (those driven by magma ascent, and that will culminate into an eruption^24,25), failed/aborted magmatic unrests (in which degassing magma is involved that will not ultimately erupt²⁶), and hydrothermal unrests (in which no magma is involved but hydrothermal/phreatic explosions are still possible²⁷). A quantitative understanding of the mechanisms underpinning periods of degassing unrest is required to interpret and manage an evolving unrest situation in real-time.

Some of the most puzzling degassing unrests have recently occurred at La Fossa volcano (Fig. 1), a relatively young (<6 kyr old^28,29,30) tuff cone on Vulcano Island (Italy). Since its last³¹ vulcanian eruption in 1880–90, La Fossa has recurrently been in a state of unrest, motivating extensive geochemical work and leading to the acquisition of an uniquely robust gas dataset²³. A common conclusion of these studies^19,21,23 is that intense fumarolic activity is sustained by mixing between CO₂-N₂-He-rich magmatic gases and hydrothermal steam deriving from evaporation of seawater-fed hot and saline fluids. Periods of elevated unrest occurred in the 1920s (ref. ³²) and in 1988-93 (ref. ¹⁹), this latter phase being thought³³ to have started in late 1977. Lower intensity periods of unrest were detected in 1996, 1998, and 2004–2009 (refs. ^23,34,35). These unrest periods have been characterised by volcanic gases becoming increasingly magmatic in composition^19,21,23,36 (with CO₂-N₂-He concentrations increasing relative to H₂O), by visible escalations in outgassing flux^37,38 (Fig. 2a), by increasing fumarole temperatures (up to ~650 °C in 1991; ref. ¹⁹) and by increasing diffuse degassing and thermal output from the volcano’s periphery³⁹. Notably, mechanisms proposed to explain these periodic, abrupt unrests include (i) achievement of some pressure threshold in the hydrothermal system, due to accumulation of exsolved fluids sourced from stationary magma cooling and crystallising at depth^23,36,40, or (ii) increasing magmatic gas supply from intermediate (latitic) magma undergoing decompression and open-system degassing during shallow ascent²¹.

Fig. 1: Escalating degassing activity at La Fossa cone in 2021.

a The La Fossa volcanic plume seen from the Vulcano Porto (harbour) area in November 2021; b to (d), evolution of the fumarolic field, from September 2021, right before onset of the unrest (b), to during the unrest in November 2021 (c) and February 2022 (d) (Photos S.C.).

Fig. 2: Temporal record of gas flux and fumarole composition from La Fossa crater.

a SO₂ fluxes from La Fossa. Green triangles are literature data^{37,49,50,51,52,97,117,118,119,120}. Our 2021 results (Supplementary Table 1) are shown as red-bordered green triangles. The time-averaged 1995–2020 flux is taken as representative of background degassing. b A magnified section of the time series (June-December 2021) shown in (a), showing the acceleration in SO₂ degassing since October 2021. The error bars correspond to 1 standard deviation of the daily averaged flux. The Multi-GAS derived CO₂/S_T ratios of panel (c) are also shown (as circles); (c) Temporal record of fumarole molar CO₂/S_T ratios. Notably, the SO₂ flux refers to the bulk (total) plume, while the reported gas compositions are for individual fumaroles. We integrate results from direct sampling of the fumaroles (red crosses, refs. ^19,36,40; solid blue line, ref. ³⁴; blue crosses, ref. ²³) and from Multi-GAS sensing of near vent plumes (open circles, 2005–2010 results^{55,96,97,118,119,120}; 2015–2021 are from this study, Supplementary Table 2); d Temporal record molar H₂S/SO₂ ratios in fumarole emissions (same sources as panel c).

In mid-September 2021, a new abrupt, intense acceleration in the volcanic gas output from La Fossa was recorded⁴¹ (Fig. 1). Reports^42,43 from Istituto Nazionale di Geofisica and Vulcanologia soon revealed that the vented gases were among the most magmatic (CO₂-rich) ever detected, and that the gas output escalation was also accompanied by shallow long period seismicity and, for the first time on record, by centimetre-scale ground inflation centred on La Fossa crater. In the following weeks, degassing activity remained sustained at the summit fumaroles. At the same time diffuse CO₂ release from the volcano’s periphery also escalated, raising concerns regarding potential asphyxiation hazards in the inhabited areas. Local authorities were forced to prohibit access to the Vulcano Porto village overnight, capturing media attention worldwide⁴⁴. The unrest is still ongoing by the time of writing, with reduced intensity⁴⁵.

During background activity, the relatively low (≤20 tons/day) sulphur dioxide (SO₂) flux^37,41 sustained by the high-temperature (∼300–400 °C) permanent fumaroles suggests pervasive hydrothermal circulation in the volcano’s subsurface. The 2021 unrest, as in 1988–93 (ref. ³⁷), was associated with a remarkable acceleration⁴¹ in the SO₂ flux up to >200 tons/day (Fig. 2a). Although these sin-unrest SO₂ flux levels remain well below those charactering the strongest degassing volcanoes worldwide (where SO₂ is typically discharged at rates of several thousands tons/day; ref. ²), they still require an excess sulphur source. However, the origin of sulphur (S) in La Fossa fumaroles has long been a matter of debate; previous interpretations invoke either a magmatic^46,47, hydrothermal⁴⁸ or mixed¹⁹ magmatic/hydrothermal origin. These interpretations have distinct volcanic hazard implications and therefore interpreting the cause for the elevated SO₂ flux is central to defining a coherent model for the initiation, and evolution, of La Fossa degassing unrests.

Here, we use volcanic gas records (in situ plume compositions and fluxes from remote sensing) to resolve the magmatic versus hydrothermal contributions to the increased S output in 2021. By comparing observations with results of a degassing model, initialised at various P-T-X conditions relevant to Vulcano Island volcanism, we propose that degassing of mafic magma—either due to ascent and decompression at depth or through the release of previously accumulated gas at the roof of a ~5 km deep reservoir underneath La Fossa—triggered the 2021 unrest.

Results and discussion

SO₂ emissions

We quantify the SO₂ flux from La Fossa cone during the 2021 unrest from ground-based remote dual UV Camera measurements (Supplementary Fig. 1; see Methods). Although discontinuous, our UV Camera data (Supplementary Table 1) record an escalation in SO₂ release during fall 2021 (Fig. 2a, b). Our measurements do not constrain the exact onset of the degassing unrest (fixed at September 13 based on geodetic evidence⁴³) but instead bracket the transition period from background to unrest regimes. We demonstrate that, in June 2021, the SO₂ flux was low at 10–16 tons/day (Fig. 2b); this is very close to the background SO₂ emission levels from La Fossa during 2008–2016, as inferred from the mean (18 tons/day) and median (22 tons/day) fluxes derived from regular observations with a permanent scanning spectrometer^49,50 (Fig. 2a). We find that, in mid-October 2021, the SO₂ flux increased by a factor >10 (190–256 tons/day on October 12), and then stabilised at a flux well above previous background levels (73–282 tons/day) in late October-December 2021 (Fig. 2b). A similar SO₂ escalation has recently been reported from results of a network of permanent UV scanning spectrometers⁴¹. Our results constrain the average SO₂ flux during the unrest (October-December 2021) at 143 ± 60 tons/day. For comparison, the mean SO₂ flux during the 1988–93 unrest was ~88 tons/day (inferred from averaging data in refs. ^37,51,52). Similarly, the peak SO₂ emission was 155–240 tons/day in June 1991 (ref. ⁵¹) and subsequently declined since 1992, with background emissions being restored by 1995 (Fig. 2a). We therefore consider the median SO₂ flux of ∼23 ± 15 tons/day for the 1995–2020 period, in which the La Fossa showed only mild and short-lived degassing unrests^19,34,35,53, as representative of the background La Fossa emissions (Fig. 2a).

Changes in gas composition

Changes in the H₂O-CO₂-SO₂-H₂S composition of the near-source atmospheric plume emitted by the La Fossa crater fumaroles, as measured by a portable Multi-Component Gas Analyser System (Multi-GAS) (Supplementary Fig. 2; see Methods), indicate that - in addition to the increasing SO₂ flux (Fig. 2a, b)—the major element composition of the atmospheric plume changed markedly in October-December 2021, relative to the 1995-early September 2021 period (Supplementary Table 2; Fig. 2b–d).

We find that fumarolic emissions in October-November 2021 have anomalously low molar CO₂/S_T ratios (where S_T = SO₂ + H₂S; Fig. 2b, c) compared to the typical background composition. The range in molar CO₂/S_T in October-November 2021 is large, from 4.2 to 38.6, reflecting both temporal variability (Supplementary Table 2; Fig. 2b, c) and, more importantly, the compositional heterogeneity of the fumarolic field (Fig. 3). Walking Multi-GAS surveys conducted in November 16–18 demonstrate the spatial diversity within the La Fossa fumarolic field, in which both S-rich (CO₂/SO₂ ratio <10) and S-poor (CO₂/SO₂ ratio >30) fumaroles coexist (Fig. 3). Figure 3 also highlights, however, a core of low CO₂/SO₂ ratio (and CO₂/S_T ratio) fumaroles in the central segment of the outer rim fumarolic field, extending along the F0 fumarole and partly including the FA-F5AT-F11 fumaroles (Fig. 3). In the context of previous Multi-GAS (circles) and direct sampling (crosses) observations (Fig. 2c), these October-November 2021 compositions correspond to some of the lowest CO₂/S_T ratios measured at La Fossa, comparable only to those seen during the 1988-93 unrest (Figs. 2c, 4a; refs. ^19,36).

Fig. 3: Mapping the chemical diversity of the La Fossa fumarolic field.

The map shows the spatial heterogeneity of the fumarolic CO₂/SO₂ ratio in the fumarolic field, and is obtained by post-processing the results of the 16–18 November 2021 Multi-GAS walking traverses using the methodology of ref. ¹⁰⁰. This processing routine allows resolving the spatial heterogeneity of the fumarolic field with high temporal/spatial resolution, and is thus complementary to the standard protocol used to derive the ratios listed in Supplementary Table 2 (in which all data acquired nearby any single fumarole are jointly analysed). Numbers in the map identify the position of the fumaroles for which compositional results are shown in Supplementary Table 2. The base map used in the figure is from sources: Esri, DigitalGlobe, GeoEye, i-cubed, USDA FSA, USGS, AEX, Getmapping, Aerogrid, IGN, IGP, swisstopo, and the GIS User Community.

Fig. 4: Volcanic gas composition at La Fossa crater.

a Scatter plot comparing the October-December 2021 Multi-GAS results with directly sampled fumarole data taken in 1988–91, at the peak of the 1988–93 unrest¹⁹. The 1988–91 fumarole data identify two linear trends (for inner crater and crater rim fumaroles, respectively) that converge toward a magmatic pole (purple star labelled M) characterised by low H₂O/CO₂ and S_T/CO₂ ratios (the magmatic end-member of ref. ¹⁹). The October-December 2021 data (this study) match this magmatic signature; b H₂O/10-CO₂-5·S_T triangular plot, comparing the composition of the La Fossa October-December 2021 fumarolic plume samples (Supplementary Table 2) with previous Multi-GAS (Supplementary Table 2) and direct sampling (see legend for data provenance) results. The fall 2021 data cluster in a relatively narrow area of unusually high S_T, and low H₂O, compositions. The curves show the model-predicted pressure-dependent evolution of the magmatic gas phase at equilibrium with silicate melts, obtained during the degassing model runs (see Supplementary Table 4; curve B6 not shown for clarity but see Supplementary Fig. 4). Stars identify equilibrium gas at initial pressure conditions. The B1-B5 model curves are for a basaltic magma with 1–2 wt.% initial CO₂ content and redox conditions from NNO-1 to NNO + 1 (see Supplementary Table 4). Models L1 and T1 are for latitic and K-trachytic magmas, and R1 for a rhyolitic magma composition (see Supplementary Table 4); c 0.5CO₂-3SO₂-5·H₂S triangular plot for the La Fossa fumaroles (symbols and curves as in (b)). The compositional fields of typical volcanic gas types/processes are from ref. ²⁷). The fall 2021 La Fossa gases vary compositionally along a relatively narrow area stretching between the fields of deep magmatic gases and hydrothermal gases (having suffered S loss via hydrothermal scrubbing processes). The bars at the bottom of the diagram illustrate the calculated equilibrium gas H₂S/SO₂ ratios for hydrothermal gases buffered at the Hematite-Magnetite (HM) buffer at 0.1 to 20 MPa pressure (ref. ⁴⁸), and for steam-driven pyrite oxidation¹⁹.

Our Multi-GAS results also indicate that SO₂ prevails over H₂S at all fumaroles during October-December 2021, with relatively uniform H₂S/SO₂ ratios averaging 0.09 ± 0.03 (1σ; range, 0.05–0.18; Supplementary Table 2 and Fig. 2d). Thus, the decrease in CO₂/S_T ratio (Fig. 2c) essentially reflects increasing SO₂ in the fumarolic emissions, consistent with the SO₂ flux increase observed (Fig. 2a). In addition to a low CO₂/S_T signature, the October-November 2021 emissions are also particularly H₂O-poor, with measured molar H₂O/CO₂ ratios (1.2–4.1) being the lowest on record (Fig. 4a). Multi-GAS data from December 2021 exhibit more variable composition than in previous months (Fig. 4b), together with a gradual progression toward higher CO₂/S_T ratios (11.6-70; Supplementary Table 2, Figs. 2b, 4a, b) at similarly low H₂O/CO₂ ratios (of 1.5–4.7).

Overall, the October-December 2021 results identify a compositional field that overlap the magmatic gas compositional end-member identified by directly sampled 1988–91 fumaroles¹⁹ (Fig. 4a).

An especially SO₂-rich volcanic gas in fall 2021

Our gas observations provide evidence for a particularly SO₂-rich fumarolic composition in fall 2021. The October-December 2021 plume samples plot at the high SO₂-low H₂S end (Fig. 4c) of the range exhibited by the 1979–2004 La Fossa fumaroles^19,36,40, which encompasses the 1988–93 unrest. They are also substantially more SO₂-rich than the typical background composition during the 2005–2020 quiescent phase, measured using the same Multi-GAS technique used here, and with direct fumarole sampling²³ (see Supplementary Table 2 for data sources; Fig. 4c). Further, the fall 2021 plume data cluster in a relatively narrow compositional domain of unusually high S_T and low H₂O (Fig. 4b). Notably, the October–November 2021 fumarolic emissions are similar, or even less hydrous than, the compositions observed at the climax of the intense 1988–93 unrest^19,36 (Fig. 4a).

Together with the coeval SO₂ flux increase⁴¹, these atypical gas compositions highlight an escalation in SO₂ transport within—and release from—the La Fossa fumarolic system, and require an additional source of S in excess to that sustaining the background SO₂ emission of ∼23 ± 15 tons/day (Fig. 2a). The fumarolic gases released from La Fossa have been long been attributed to binary mixing between a CO₂-N₂-He-rich magmatic endmember and H₂O-rich hydrothermal vapours deriving from evaporation of hot (300–400 °C), seawater-derived brines^19,21,23,36. The origin of S has been more difficult to recognize but geochemical models suggest a mixed magmatic-hydrothermal derivation^19,48. CO₂ and S_T concentrations were found⁴⁸ to be highly correlated in contemporaneously sampled fumaroles during 1988-91, indicating that S was primarily of magmatic origin during that time. However, during the final stage of the same degassing unrest between 1992 and 1995, modelling instead supports a prevalent hydrothermal S origin⁴⁸. From a hazard perspective, it is imperative to determine whether the enhanced transport and release of S as SO₂ during fall 2021 is driven by hydrothermal or magmatic processes.

Hydrothermal sulphur remobilization

Sulfides, especially pyrite, are common in the sub-soil of La Fossa cone⁵⁴. Thus, when hot CO₂-rich magmatic fluids are supplied during an unrest phase¹⁹, it is possible that S is remobilized into the gas phase by either steam-driven oxidation¹⁹ or brine-mediated dissolution⁴⁸ of previously-deposited pyrite. Remobilization of mineral- (pyrite-) hosted hydrothermal S is, at least in principle, a viable mechanism to explain the SO₂ escalation seen during degassing unrests (e.g., in 1988–91 and in fall 2021; Fig. 2). Here, we evaluate this potential S source in the context of our gas compositional data.

The gas molar H₂S/SO₂ ratio can in principle be used²⁷ to resolve the S fraction produced in the hydrothermal envelope. This is because, in the prevailingly reducing hydrothermal conditions (relative to more oxidised magmas) H₂S is favoured over SO₂ (Fig. 4c; ref. ²⁷). The H₂S/SO₂ ratio has shown large fluctuations at La Fossa over the last 45 years (Fig. 2d). However, using this ratio to resolve magmatic vs. hydrothermal gas provenance is hampered at La Fossa by the relatively oxidized state of the hydrothermal system (thought to be buffered by coexisting hematite and magnetite; i.e., at the Hematite-Magnetite [HM] buffer⁴⁸). In dry (steam-only) conditions, and at such high redox condition, calculations show that the hydrothermal gas H₂S/SO₂ ratio, buffered by pyrite oxidation, can lie anywhere between ~160 (at 10 MPa and ~650 °C) and ~0.1 (at 0.1 MPa and ~700 °C) (ref. ¹⁹) (Fig. 4c). Similarly, in models assuming the La Fossa hydrothermal reservoir is a bi-phase (brine+steam) system at ~20 MPa and 300–400 °C (ref. ⁴⁸), the equilibrium (HM-buffered) H₂S/SO₂ gas ratio would vary between 150 (300 °C) and 0.4 (400 °C), well above the ratios measured in the surface emissions during the 2021 unrest (0.05–0.18; Fig. 4c). However, if hydrothermal HM buffering prevails until the vapour phase is discharged at the surface (i.e., at 0.1 MPa) then the equilibrium H₂S/SO₂ gas molar ratios vary from 0.77 (300 °C) to 0.02 (400 °C), which encompass the observed fumarolic range (Fig. 4c). Therefore, in both models (dry or two-phase), an oxidised La Fossa hydrothermal system can buffer H₂S/SO₂ in the gas at values comparable to those obtained for a magmatic gas buffered by shallow (<50 MPa) magma at the Nickel-Nickel Oxide buffer (NNO; ref. ⁴⁸) (Fig. 4c). Consequently, as previously noted⁵⁵, the H₂S/SO₂ molar gas is an unsuitable tracer to resolve hydrothermal from magmatic contributions at La Fossa, and the relatively low ratios measured in October–December 2021 (Fig. 2c) cannot be taken in favour for the prevalence of any of the two potential S sources.

The S yield of pyrite hydrolysis reactions is strongly temperature-dependent¹⁹. Therefore, if increasing gas-phase S remobilization during magmatic-gas driven heating of the hydrothermal system is underpinning the enhanced S transport observed during the 2021 unrest (Figs. 2, 4) then a relationship between temperature and fumarolic S content is expected. However, an inspection of published fumarole compositions from La Fossa (1978–2010; Fig. 5) shows that neither the fumarolic molar CO₂/S_T, nor the total S content (S_T), correlate with fumarole discharge temperature (Fig. 5 and Supplementary Fig. 3). Rather, we highlight a clear temporal lag between the S_T peak (occurring in the 1988–90) and the temperature peak in 1993 (at FA fumarole), when S_T was already decreasing (Fig. 5). Similarly, the SO₂ emission flux reaches its maximum nearly two years before (in June 1991) the 1993 FA temperature peak (Fig. 5).

Fig. 5: Chronology of sulphur release from La Fossa crater.

A review of literature data on fumarole temperatures, SO₂ fluxes, and fumarole total S content (S_T). SO₂ flux results (triangles) are from same source as Fig. 2. Data from this study (Supplementary Table 1) are identified by red bordered triangles. Total S content (S_T) in directly sampled fumaroles (red crosses, refs. ^19,36,40; solid blue line, ref. ³⁴; blue crosses, ref. ²³). Fumarole temperatures are from ref. ^19,36 (circles, F1 and F5 fumaroles) and from ref. ³⁴ (FA and F11 fumaroles). A clear temporal shift is observed between the peak in sulphur release (both SO₂ fluxes and S_T) and the temperature peak (the latter lagging behind >2 years in time).

In summary, while we do not question that hydrothermal reactions can source some, if not all, the background SO₂ flux from La Fossa (when the CO₂/S_T is also high; Fig. 2a, b), we see no evidence in the published volcanic gas record that the elevated SO₂ transport during the 1988-93 unrest was linked to heating of (and S remobilization from) the hydrothermal system. Ultimately, considering that the fumarolic temperature increase in fall 2021 is more limited (<100 °C; refs. ^42,43) than in 1988–93, and that temperatures stabilised, or even started to decrease, by December 2021 (ref. ⁴³) while the SO₂ flux remained elevated, we consider it unlikely that the enhanced transport of SO₂ (and corresponding decrease in CO₂/S_T) are controlled by temperature-dependent hydrothermal reactions to any major extent.

SO₂ release from La Fossa magmas

We next test the hypothesis that the elevated SO₂ emissions at La Fossa during fall 2021 are sourced from magmatic degassing. We model S partitioning from silicate melt into the magmatic gas phase using a H-O-C-S volatile saturation model^56,57 (see Methods) initialised at P-T-X conditions relevant to the La Fossa magmatic system. As it is becoming increasingly evident at arc volcanoes worldwide⁵⁸, the Vulcano Island magmatic plumbing system (Fig. 6a) is complex and vertically-elongated, with multiple magma ponding zones distributed over a range of depths/pressures^{31,35,59,60,61}. This polybaric system is thought to have evolved in time during the >100 kyr history of the volcano²⁸ and can be described schematically by (i) a deep storage zone(s) at mid (17–12 km) to lower (18–21 km, Moho limit) crustal depths, occupied by mafic (basaltic) magma; (ii) an intermediate storage reservoir at ∼3–5 km (the level of neutral buoyancy of ascending magmas), where intermediate (latitic) magma develops from basalt differentiation (fractional crystallization and assimilation), and (iii) a very shallow (1–2 km depth), heterogeneous (trachytic to rhyolitic) storage zone beneath La Fossa (Fig. 6a).

Fig. 6: Pressure-dependent volatile evolution in the La Fossa plumbing system.

a A schematic cross-section of the La Fossa magmatic plumbing system, showing the inferred storage depths of the various magma types. The deepest reservoir (basaltic to shoshonitic) is located between 21 and 18 km, at the mantle-crust boundary^59,121. A shallower shoshonitic to latitic storage zone is hypothesized between ~13 and 8 km (refs. ^28,59,108). The shallow crustal system is characterized by latitic to trachytic and rhyolitic magmas at a depth comprised between 5.5 and 1.3 km (refs. ^28,31,59,60). b Model results showing the pressure-dependent evolution of dissolved S in the melt, for all model runs; stars indicate the initial conditions of each model run (Supplementary Table 4). The model curves illustrate the declining dissolved S contents predicted by isothermal decompression. c Model-predicted pressure-dependence of the magmatic gas CO₂/S_T ratio at equilibrium with La Fossa melts in the T-X compositional domain explored. The magmatic gas CO₂/S_T decreases during the decompression path in each model run. The vertical grey bar identifies the observed CO₂/S_T ratio range (7.5 to 11.5) of the most H₂O-poor (most magmatic) La Fossa fumarolic plumes. Intersection between the average CO₂/S_T ratio of 9.9 (indicated by the vertical dashed black line) and degassing model curves identifies the possible pressures (depths) of gas-melt separation for the different model scenarios. These range from ~55 MPa (B1 run) to ~20 MPa (B2 and B4 runs) and ~10 MPa (B5 run) for the basaltic magma scenarios, and from ~20 MPa (L1), <5 MPa (T1) and ~22 MPa (R1) in the other magma type scenarios.

To account for this trans-crustal magmatic system, we initialise the model runs by ad-hoc sampling the Vulcano Island melt inclusion (MI) population (Supplementary Table 3) to cover the full spectrum of magma compositions (basaltic-latitic-trachytic-rhyolitic; Fig. 7). The model run conditions are summarised in Supplementary Table 4. Runs B1–B6 (Supplementary Table 4 and Supplementary Fig. 4) simulate the decompression-driven, closed-system degassing of basaltic magma (represented by the La Sommata scoriae erupted at ~42 ka; ref. ⁶²) on ascent from the deep Vulcano Island magma storage zone (initial pressure of 300 MPa). We investigate isothermal decompression paths from a range of starting conditions, including redox conditions from NNO-1 to NNO + 1 and initial CO₂ contents between 1 and 2 wt. % (Figs. 6, 7, Supplementary Fig. 4). The L1 and T1 runs refer to latitic (represented by Vulcanello I tephra, ~X century AD; ref. ³⁰) and K-trachytic (PAL-D eruption, ~XIII century AD; ref. ³⁰) magmas, respectively (each run with distinct initial S and H₂O contents to cover the spread of volatile contents in MIs; Fig. 7), degassing in open-system conditions²³ upon their ascent from the intermediate magma storage level (initial pressure, 50 MPa). Lastly, the R1 run shows the degassing path for rhyolitic magma (represented by PAL-B eruption occurred at ~XI century AD; ref. ³⁰) undergoing open-system degassing upon decompression from 40 to 0.1 MPa (see Supplementary Table 4 and Methods). An open-system degassing mode is favoured for these runs as implied by the relatively shallow pressure (depth) conditions of magma emplacement, as in previous study²³.

Fig. 7: Melt inclusion data compared to degassing models.

a SiO₂ vs. S scatter plot for the melt inclusion Vulcano island populations (data from Supplementary Table 3). We observe a drop in MI S contents from basaltic to intermediate (latitic) and silicic (trachytic to rhyolitic) melt compositions. The large stars labelled B1–B6, L1, T1 and R1 indicate the initial compositions of the degassing model runs (Supplementary Table 4). The solid red lines indicate the range of S contents predicted in the basaltic model runs (B1–B6) upon decompression degassing in the 300–0.1 MPa pressure range. Small stars indicate the model-predicted dissolved S content in basaltic melts at 50 MPa at the B1 (upper star) and B6 (bottom star) run conditions. b H₂O vs. S scatter plot for the melt inclusion Vulcano island populations (data from Supplementary Table 3). The solid lines indicate the model-predicted dissolved H₂O and S contents in the melt along the modelled degassing paths. Stars as in (a) with the addition of a small blue star indicating the modelled dissolved S content in the L1 model latitic melt at 20 MPa.

Modelled decompression paths predict a strong pressure-dependence for the dissolved S content in magmas (Fig. 6b) and the gas CO₂/S_T ratio in equilibrium with the melt (Fig. 6c). Dissolved volatile compositions in melts during degassing runs B1/B4 (at NNO-1) and B2/B5 (at NNO) reproduce the H₂O-S degassing trends recorded in basaltic MIs (see model curves in Fig. 7b) more closely than the B3/B6 runs (at NNO + 1). Further, the model runs initiated under more oxidising conditions (B3/B6) predict limited S partitioning to the gas phase at high pressures (Fig. 6b) that is not reflected by the MI dataset (Fig. 7b). These results suggest that redox conditions between NNO (refs. ^48,63) and NNO-1 (Fig. 7b), rather than NNO + 1, prevail in the magmatic system. The results of the B3/B6 model runs are thus considered for illustrative purposes only in the following discussion (Figs. 4, 8; Supplementary Fig. 4). The L1, T1 and R1 model runs (Fig. 6b) match well the low H₂O-S end of the MI populations (Fig. 7b).

Fig. 8: Resolving hydrothermal vs. magmatic gas contributions.

A scatter plot of H₂O/S_T vs. CO₂/S_T ratios in La Fossa fumaroles (from direct sampling; crosses, symbols as in Fig. 4) and in their near-vent atmospheric plumes (from Multi-GAS, data from Supplementary Table 2; circles, symbols as in Fig. 4). The October–December 2021 La Fossa fumarolic plumes cluster in a compositional band overlapping the magmatic gas field (red coloured area) delimited by the B1/B4, B2/B5 and L1 model degassing trends (curves as in Fig. 4). The positive H₂O/S_T vs. CO₂/S_T relationship, exhibited by the October-December 2021 datapoint array (Fig. 8), implies some extent of mixing with hydrous (and S-poor) hydrothermal fluids. The magmatic gas mixing proportion (green dashed lines) are estimated at 30–100% by resolving mixing (mixing curves shown as dashed back lines) relationships between a magmatic gas end-member (white star M; calculated averaging the composition of the most H₂O-poor Multi-GAS data) and a possible set of hydrothermal end-members (orange squares; identified using the gas samples with the most extreme H₂O-rich compositions). Grey dashed lines illustrate mixing between the magmatic gas end-member M and the 1970s–1990s hydrothermal gas end-members identified by ref. ³⁶ (purple squares) and ref. ⁴⁸ (golden squares). The blue-coloured fields thus identify the potential compositional fields of the hydrothermal component(s) in the 1970s–1990s and 2005–present.

The theoretical magmatic gas compositions coexisting in equilibrium with the modelled melts—for all melt compositions explored—evolve from being CO₂-rich at initial pressures (300 to 50–40 MPa, depending on the runs) to H₂O- and S-rich at 0.1 MPa (Fig. 4b, Supplementary Fig. 4), reflecting the pressure-dependent solubilities of the major volatile species in silicate melts^64,65. This compositional evolution is reflected in the decreasing CO₂/S_T ratio in the fluid phase with decreasing pressure (Fig. 6c).

Magmatic versus hydrothermal mixing proportions, and a pressure estimate for the magmatic source

The predicted magmatic gas compositional range circumscribed by the B1/B4, B2/B5 and L1 modelled degassing paths overlap with the range of volcanic gas compositions measured between October and December 2021 (Figs. 4b, 8 and Supplementary Fig. 4). The T1 and R1 runs produce, instead, modelled gas compositions that are more hydrous than the observed fumaroles (Figs. 4b, 8).

This agreement between modelled (B1/B4, B2/B5 and L1) and measured gas compositions is strong evidence that the unusually H₂O-poor (and S-rich) gas composition in October–December 2021 reflects venting of predominantly magma-derived fluids (see also Fig. 4a). However, the positive relationship between H₂O/S_T and CO₂/S_T defined by the October–December 2021 array (Fig. 8) also implies some degree of mixing between magmatic and hydrous (and S-poor) hydrothermal fluids^19,21,23,36.

We estimate the composition of the magmatic gas end-member by averaging the most H₂O-poor (i.e., those with the smallest hydrothermal component, if any) plume samples in our dataset. This end-member magmatic gas (white star labelled M in Fig. 8) contains ~68.4 mol.% H₂O, ~28.7 mol.% CO₂, and ~2.9 mol.% S_T, and is less hydrous than estimated previously for the 1988–93 unrest (see purple star labelled M in Fig. 8; ref. ¹⁹).

Identifying the composition of the hydrothermal component is more challenging, as previous work has found this component to vary in composition over time⁴⁸. The inferred hydrothermal compositional range for the 1970s–1990s (refs. ^19,48) is illustrated in Fig. 8; however, mixing (see grey dashed lines in Fig. 8) between such hydrothermal fluids and our magmatic end-member cannot reproduce the October–December 2021 data array. Instead, the 2021 array, together with the 1999–2010 fumarole compositions²³, suggest the involvement of a more S-poor hydrothermal component (e.g., with higher H₂O/S_T and CO₂/S_T ratios than in the 1970s–1990s) between 1999 to present. Although determining the exact hydrothermal fluid signature is beyond the scope of this study, we identify several possible end-member compositions based on the samples exhibiting the most hydrothermal (H₂O-rich, S-poor) signatures and derive the corresponding magmatic-hydrothermal mixing curves (see black dashed curves in Fig. 8). From these mixing curves, and considering the compositional variation of fumarolic plume samples, we infer that between 30% and 100% of the S released by La Fossa fumaroles during the October–December 2021 unrest is of magmatic origin (Fig. 8). These calculations suggest that the fumarole-specific magmatic vs. hydrothermal mixing proportions control the spatial heterogeneity of the fumarolic field (e.g., the CO₂/S_T variability seen in each survey; see Figs. 2a, b and 3), as previously demostrated^21,23. Further, based on average plume compositions, we estimate that the magmatic sulphur contribution reduces from ~80% during October-November 2021 to ~50% in December 2021 (Fig. 8).

If, as the above modelling suggests, fumarolic emissions from the October–November 2021 degassing unrest are indeed prevalently magmatic, then their CO₂/S_T ratio—at least that of the H₂O-poorest samples, in which the magmatic gas proportion is estimated to be >80%—can be used to infer the pressure of final gas-melt equilibration (i.e., the pressure at which gas separates from melt) based on the distinct pressure-dependent solubilities of the major volatile species in silicate melts^64,65,66. We show how the CO₂/S_T ratio in the exsolved fluid phase is predicted to vary as a function of pressure, for the various model starting conditions for La Fossa (Fig. 6c). The least hydrous (most magmatic) La Fossa plume samples yield CO₂/S_T ratios between 7.5 and 11.5, with an average of 9.9 ± 1.0 (see Fig. 6c). This mean ratio translates into a range of possible saturation pressures depending on the chosen model run conditions shown in Table 1. For a basaltic magma starting composition, we estimate equilibrium pressures ranging from ~55 ± 5 MPa (B1) to ~20 ± 2 MPa (B2 and B4) and ~10 ± 1 MPa (B5) (Fig. 6c). For a latitic magma, the magmatic gas CO₂/S_T ratio of 9.9 suggests gas-melt equilibrium pressures of ~19 ± 2 MPa (L1). Lastly, decompression of a K-trachyte and rhyolite magma would imply degassing pressures of <5 MPa and ~22 ± 2 MPa, respectively (Fig. 6c).

Table 1 Results of magma budget calculations.

Constraints on the magmatic source from magma budget calculations

Our comparison between modelled and observed gas compositions, although clearly supporting the involvement of magmatic volatiles in the unrest, does not allow us to unambiguously resolve the nature (composition) and storage conditions (emplacement/storage depth) of the magmatic source. The T1 and R1 modelled gases exhibit hydrous compositions, and therefore do not match the observed gas compositions in fall 2021 (Figs. 4, 8). In addition, trachytic to rhyolitic magmas, such as those erupted during the ∼1200 AD and 1888–90 La Fossa eruptions, are volatile-poor (Fig. 7) and have more radiogenic (³He/⁴He ratios (R) normalised to the air ratio (Ra) of ∼3.2 to ∼4.8 Ra/Ra) He isotopic signature⁶⁷ than present-day fumaroles (∼6 R/Ra), making them unsuitable candidates as potential sources of the 1978-to-present La Fossa fumarolic activity^21,23,68,69. Instead, fluid inclusion studies demonstrate that only more deeply-stored magmas, such basalts to latites (Fig. 6), have He isotopic signatures consistent with those of present-day gas emissions⁶⁷; we therefore focus on these two magma compositions in the subsequent discussion.

The question of whether mafic or intermediate (latitic) magma is currently feeding recent (1970s to present) degassing activity at La Fossa is critical to future volcanic hazard assessment³⁵; however, no consensus has thus far been reached. Volcanic gas geochemical modelling suggests that the temporal evolution of fumarole composition (both chemical and isotopic) is best explained by degassing of a shallow (<40 MPa), vertically-elongated latitic magma reservoir²³. However, volatile contents within melt inclusions indicate that shallow-stored (latitic to rhyolitic) magma underneath La Fossa is typically S-poor (see Fig. 7) and therefore that degassing of this magma source is insufficient to sustain the long-term La Fossa SO₂ (and CO₂) outgassing budget⁴⁶. Consequently, ascent and degassing of deep-derived (>50 MPa) mafic (basaltic) melt has been invoked to explain the increased SO₂ output during the 1988–93 degassing unrest⁴⁶.

Together with the results of model calculations (Figs. 6–8), our 2021 SO₂ flux time series (Fig. 2a) offers insight into the source characteristics of the parent magma. The time-averaged SO₂ flux during October–December 2021 is 143 ± 60 tons/day. At this rate, we estimate that a total SO₂ mass of 1.6 × 10⁴ tons was outgassed between the unrest onset (September 13, 2021; ref. ⁴³) and December 31. If we assume, conservatively, that 80% (September–November) and 50% (December) of the degassed SO₂ is magmatic (Fig. 8), then the total magmatic SO₂ mass reduces to 1.1 × 10⁴ tons (Table 1). Using the typical S yield (mg of S released per kg of magma) for the various magma types, as inferred from difference between initial (from melt inclusions) and final (at emplacement depth/pressure; from model run calculations; Fig. 6b) melt S contents, we show that the volume of degassing magma required to explain the surface SO₂ emissions during the period of unrest ranges from 1.8·10⁶ (B4) to 7.2·10⁶ (L1) m³ (see Table 1 for full mass balance calculations).

The mass balance calculations presented above assume that the increased S emissions, and ultimately the degassing unrest itself, are triggered by decompression-driven volatile exsolution and degassing of a magma batch, ascending from reservoir conditions (Fig. 6a) to final emplacement pressure (Fig. 6c), where gas-melt separation finally occurs. In this interpretation, we can use a simple Mogi modelling approach⁷⁰ to calculate the expected uplift at barycentre (e.g., at La Fossa crater) from magma emplacement volume and depth (see Table 1). Although simplistic in that magma compressibility is not considered in the calculations⁷¹ and that the deformation is assumed to be a simple point source, the calculated uplift ranges from ≥50 cm (B5, L1) to ~3.6 cm (B1). The measured uplift during the unrest was ~3 cm (refs. ^42,43), thus suggesting that only the B1 scenario is geologically plausible.

In summary, we conclude that our gas observations, combined with geodetic information (refs. ^42,43), can consistently be explained by the emplacement of ∼3·10⁶ m³ of mafic magma at ∼4–5 km depth (Model B1; Table 1). We stress that the calculations above are in no way meant to quantitatively model the deformation pattern at La Fossa – they should merely be intended as back-of-the-envelope estimates of the expected order-of-magnitude uplift in the different magma scenarios explored.

Migrating or stationary magma?

Models involving the movement of magma masses²¹ beneath La Fossa, as implied by our calculations (Table 1), have recently been challenged in light of the lack of supporting geophysical evidence^23,72,73,74. During hydrothermal background degassing activity, such as between 1999 and 2013, modelling⁷⁴ of ground deflation and relatively minor (infrequent) volcano-tectonic (VT) seismicity indicate the absence of magma intrusion events, and instead suggest a deflating magma reservoir at ∼3–5 km depth. In contrast, during degassing unrests the frequency of VT seismicity increases and is typically accompanied by ground uplift⁷⁴, as observed in the 1980s–1990s (refs. ^75,76) and in fall 2021 (refs. ^42,43). However, even during these unrest periods, seismicity and deformation are confined mostly within the shallow edifice and are therefore interpreted as pressurisation caused by the accumulation of magmatic gas within the shallow hydrothermal system^{74,75,76,77,78}. Similarly, pressure build-up in the ∼3–5 km magma storage zone in the late 1980s has been attributed to increased magmatic fluid pressure rather than to magma emplacement⁷⁴. In summary, geophysical evidence points to a key role for the migration of magmatic fluids (likely controlled by changes in the regional stress field), rather than magma itself, in generating the observed geophysical monitoring signals⁷⁴.

To reconcile geophysical observations with geochemical evidence, a plumbing system model has been proposed²³ in which temporal changes in magmatic gas composition are explained by time-varying gas contributions from multiple stationary, compositionally-distinct magma storage reservoirs, located at various depths within the vertically-elongated La Fossa plumbing system (Fig. 6a). In this stationary magma scenario, degassing unrest reflects gas accumulation in the roof zone of a magma reservoir (with gas bubble deriving from deeper in the magmatic system, and/or through exsolution driven by cooling and crystallisation^65,79), until such a time as this volatile-rich cap is destabilised, either by some external forcing or once some threshold pressure is reached. Following this, accumulated fluids are released abruptly and ascend through the shallow magmatic-hydrothermal system, causing unrest²³.

With our data alone, we cannot unambiguously resolve between the ascending magma or stationary magma hypotheses. Indeed, our modelling results would be entirely consistent with a model of stationary magma, in which escalating fumarolic gas fluxes in 2021 are fed by a pulse of magmatic gas released from either (i) basaltic magma at ∼4–5 km or (ii) latitic magma at ∼2 km (Fig. 6c and Table 1). We argue, however, that latitic magma is relatively poor in volatiles (Fig. 7) and is therefore unlikely to have sustained S (ref. ⁴⁶) and CO₂ (ref. ⁵³) discharges over a degassing unrest that, although varying in intensity, has lasted for more than 40 years (ref. ³⁵). Using the time-averaged La Fossa SO₂ flux of 23 ± 15 tons/day as characteristic of background periods (1978–1987 and 1994–12 September 2021), and considering the increased SO₂ flux during the two major unrest periods (1988–1993 and fall 2021), we derive a total SO₂ release of ~0.5 Mtons over 44 years (1978 to 2021). If we assume conservatively that, during background degassing, only ~10% of the degassed SO₂ is magmatic in origin and that this fraction increases to ~80% during unrest (Fig. 8), then the cumulative magmatic SO₂ release reduces to ~0.19 Mtons. Using this total mass, and assuming the extreme case of complete S loss from a latitic magma with 0.1 wt.% to 0.05 wt.% (Fig. 7b) initial S, we find that ~5·10⁷ to 10⁸ m³ of latitic magma is required to have degassed and ultimately crystallised, since 1978. Notably, these estimated volumes are minimal as they do not consider degassing during previous, poorly-documented unrest periods (e.g., 1920s) since the previous 1888–90 eruption. Although tomographic studies at La Fossa have been limited to the upper (<1.5 km) portion of the volcano⁸⁰, it is unlikely these relatively large melt volumes would have remained undetected by geophysics⁷⁴. Therefore, we conclude that, even in a stationary magma scenario, deep mafic magma must be implicated as the magmatic gas source fuelling the unrest.

Implications for the unrest

The La Fossa chronostratigraphic record is surprisingly difficult to interpret, in view of the high number of tephra layers of similar characteristics, limited thickness, and fine grain-size (favouring erosion by wind and water). This led to controversies and different interpretations in eruptive history reconstructions^28,29. However, recent integrated analysis of tephrochronological and radiometric-paleomagnetic results³⁰ has found that, since medieval time, volcanic activity has been concentrated within two main eruptive clusters (~850–1300 AD and 1420–1890 AD), in which eruptions occurred every 35 years on average (at least between AD 850 and 1550). Volcanic hazard assessment exercises³⁵ may thus need to be scaled up to consider that historical eruptive activity has been more frequent than previously thought²⁸.

The geological history of La Fossa provides copious evidence that the ascent of mafic magma has been a recurrent eruption trigger^31,46,62. Thus, tracking any mafic recharge event within the magmatic plumbing system is critically important to future hazard assessment scenarios. Our results presented here indicate that a large fraction of the elevated SO₂ emissions in October–December 2021 is magmatic in origin, and requires a mafic magma source. SO₂ flux emissions waned since January 2022 (Fig. 9) but, by the time of writing, remain well above background, with a novel phase of increase in May–June 2022 (Fig. 9), in correspondence with an abrupt phase of hydrothermal fluid venting in the Baia di Levante area⁸¹. Likewise, fumaroles’ CO₂/S_T ratios in 2022 are higher (e.g., more hydrothermally controlled; Fig. 8) than in October–November 2021, but still remain at lower levels than during the background (June-September 2021) phase (Fig. 9). These observations suggest that the unrest has waned in 2022 but not stopped, and that the mafic magma source may therefore still be actively degassing (even if at reduced rate) below Vulcano Island. Although a stationary magma scenario (in which gas is sourced by destabilisation of previously accumulated foam layer in the roof zones of a magma reservoir) cannot be excluded, modelling (Fig. 6) and calculations (Table 1) indicate that the case of mafic magma ascent (in the 300 to 50 MPa pressure range) and emplacement (at ∼4–5 km depth) should at least be considered as a possibility. For this ascending magma case, basaltic magma is likely to be relatively reduced (NNO-1) and moderately CO₂-rich (1 wt.%) (the conditions of B1 model run). We emphasise that, according to calculations, the degassing magma volume required to explain the S budget of the 2021 unrest period (e.g., that would have emplaced at ∼4–5 km depth) is relatively small (~3.6·10⁶ m³; Table 1) and may therefore be difficult to resolve from seismic-geodetic records, especially if basement rocks at 4–5 km depth are hot and plastic. Therefore, although any shallower (<4 km) magma ascent is arguably excluded by our calculations (see the unacceptably large inflation predicted by model scenarios B2, B4 and B5; Table 1), the likely presence of mafic magma at depth, combined with the pervasive hydrothermal circulation existing, make the volcano potentially unstable²⁷. For example, the reduced SO₂ flux in 2022 may reflect a combination of reduced gas supply from depth and reduced permeability (sealing)—a risky conditions at hydrothermal volcanoes²⁷. Ultimately, volcano monitoring efforts should prioritise the detection of any (further?) basaltic magma ascent and, for this, volcanic gas monitoring should remain central to detect any further escalation in magmatic S release.

Fig. 9: Evolution of the unrest through early 2022.

Temporal record of bulk plume SO₂ fluxes from La Fossa (red-bordered green triangles, data from Supplementary Table 1; error bars correspond to 1 standard deviation of the daily averaged flux) and fumarole molar CO₂/S_T ratios (circles) in the period June 2021 to June 2022. A drop of SO₂ flux emissions is visible since January 2022, followed by some moderate increase in May to June 2022. The CO₂/S_T ratios in 2022 are higher than in October–November 2021, but still below background (see June to September 2021 data for comparison).

Broader implications

Mafic magma as the S source, and trigger of unrest, is not unique to Vulcano Island but is instead common to many volcanoes worldwide. However, a remarkable aspect of unrest at Vulcano Island is its periodic nature and sudden onset: several degassing escalations have repeated in pulses over the last 45 years, each starting abruptly (as demonstrated by the 2021 example, in which degassing activity was still at background levels in early September, only a few days prior to the 13 September start of the degassing crisis). Such pulsatory degassing behaviour is perhaps only comparable to that observed⁸² at White Island volcano, New Zealand, prior to its 1976–2000 eruptive sequence, in which repeated injection of mafic magma into shallow-stored andesitic/dacitic magma fuelled the volatile output and provided the trigger for the subsequent eruption⁸³. Comparison with unrest phases elsewhere, interrupting decades to centuries of volcano dormancy, is complicated by different volcanic setting, but still brings evidence for mafic magma involvement. For example at Turrialba volcano in Costa Rica, the hydrothermal-to-magmatic gas transition (as highlighted by increasing SO₂ proportions in fumarolic emissions combined with escalating fluxes^22,24) has been more gradual, occurring over timescales of decades rather than weeks/months (as observed at Vulcano), but still fuelled by deeply rising magma²². It should be cautioned, however, that the Vulcano Island gas record is unique for its longevity and continuity; any rapid magmatic gas pulse may have remained undetected at the many quiescent volcanoes (including Turrialba) where gas observations have remained sporadic and discontinuous, at least until the advent of near-continuous instrumental monitoring²⁴. For example, mafic magma has been strongly implicated in sustaining excess S degassing during the long-lived (1995–2010) Soufrière Hills eruption in Montserrat island^84,85 and yet the volcano’s fumarolic fields had been sampled on only five occasions in the 30 years preceding eruption onset, potentially precluding detection of any mafic magma-related precursory signal⁸⁶. An important corollary of our study is that, in addition to regular monitoring, modelling S degassing behaviour in magmas⁸⁷ is critically important to detect the presence, and chemical identity, of any magmatic contribution to degassing unrests at dormant volcanoes. Combining field and petrological observations (i.e., volcanic gas and melt inclusions^88,89), modelling⁶⁵ and thermodynamic/experimental constraints⁹⁰ in a multi-parametric geochemical approach is therefore recommended to interpret the evolution of magmatic-hydrothermal unrest at volcanoes globally.

Methods

UV Camera

We use the dual UV Camera technique^91,92,93 to quantify the SO₂ flux release from the La Fossa fumarolic system. Measurements were initially conducted during campaigns (Supplementary Fig. 1a) with a portable UV Camera unit similar to the one described in ref. ⁹⁴. The instrument was powered by a 12 V battery and commanded via a Toughbook PC using the Vulcamera software⁹⁵. We acquired data for several consecutive hours per day (at 0.5 Hz rate) and the mean daily recorded SO₂ flux is listed in Supplementary Table 1 (data illustrated in Fig. 2a). On November 23, a permanent UV Camera system was field deployed in the Vulcanello area, ~2 km away from La Fossa (see Supplementary Fig. 1a for location). This fully autonomous UV camera system, similar to those in operation at nearby Stromboli⁹⁶, is powered by 12 V battery kit and solar panels, and is commanded by an embedded mini-PC (Jetway NF531M) using a set of ad-hoc designed algorithms⁹³ that allow the system to acquire 5 h per day and store data in external hard-drives. The data listed in Supplementary Table 1 correspond to daily averages calculated over temporal windows with optimal viewing conditions (clear sky, high plume visibility, un-condensed plume), while unsuitable periods (cloudy, rainy, condensed plume) are discarded. Both units (portable and fixed) are equipped with two JAI CM-140GE-UV cameras sensible to UV-radiation, fitted with two distinct band-pass optical filters (both of 10 nm Full Width at Half Maximum) with central wavelengths of 310 (SO₂ absorption) and 330 nm (no SO₂ absorption) (refs. ^93,95). The cameras acquire images (520 × 676 pixels at 10-bit resolution) at 0.5 Hz that are later post-processed using standard techniques^93,94,95,96. In brief, sets of co-acquired images are combined into sequences of absorbance images (Supplementary Fig. 1c) that are then converted into slant column amount (SCA) images using calibrations derived from either calibration cells or a co-aligned Ocean-Optics USB2000 + Spectrometer. Finally, Integrated Column Amounts (ICA) timeseries are obtained by integrating the sequences of SCA images along a cross-section perpendicular to plume transport (Supplementary Fig. 1c). A SO₂ flux time-series is finally derived by multiplying ICA by plume speed (derived by tracking the motion of plume gas fronts in image sequences, using an optical flow algorithm⁹³). Uncertainty in the derived SO₂ fluxes is evaluated at ~ ±30% (ref. ⁹¹) at the low SO₂/ash in-plume contents observed at Vulcano.

Multi-GAS

The compositions of the La Fossa fumaroles’ atmospheric plumes are characterised from periodic field surveys using a portable Multi-Component Gas Analyser System⁹⁷ (Multi-GAS). As in previous work at La Fossa^97,98,99, the Multi-GAS is used in the walking traverse mode. During operations, a light (<3 kg) portable unit is moved by an operator through the field while continuously acquiring (at 1 Hz rate), and with its position recorded by a synchronised GPS. During the measurements, the operator performs a sequence (typically 2) of walking transects through the fumarolic field with the Multi-GAS handheld, and with its inlet exposed to the fumarolic gas effluents at ~50 cm distance above each vent (in order to avoid sensors’ signal saturation). Each fumarolic vent is monitored for a few minutes, and the various fumaroles are sequentially analysed to cover the entire field in circa 1 h. During October-December 2021, degassing activity was concentrated on the La Fossa outer crater rim, while the inner crater fumaroles were only weakly degassing (and relatively colder). Our walking traverses, therefore, concentrated on the former.

In October, the very vigorous degassing activity prevented us for completing a walking traverse, and only some distal bulk plume observations were obtained (Supplementary Table 2). Walking traverses of the fumarolic field were realised in the November-December surveys. On November 7, two transects were performed in the outer crater rim area stretching between F0 and F11 fumaroles, two of the most systematically monitored fumaroles on Vulcano^23,34. The various fumaroles analyses all exhibited a pretty homogeneous gas signature on that occasion (Supplementary Fig. 2). On November 16–18 and December 14, a wider area was investigated, instead (see Fig. 3), also including a number of newly opened (fall 2021) fumarolic vents to the north of the F0-F11 fumarole clusters. The fumarolic field was found more compositionally heterogeneous in these occasions (Fig. 3 and Supplementary Fig. 2), as indicated by the ratio results for the individual fumarolic vents (GPS coordinates listed in Supplementary Table 2). Our Multi-GAS dataset is also complemented by a series of measurements carried out during routine monitoring surveys in 2009 to June 2021. These unpublished results are also listed in Supplementary Table 2. Interpretation of our Multi-GAS dataset is assisted by integration and comparison (Figs. 2, 4, 5, 8) with results of direct sampling of the fumaroles from the literature (Supplementary Fig. 3).

The Multi-GAS records in-plume gas concentrations using a sensor kit consisting of (i) one Gascard nondispersive infrared (NDIR) spectrometer from Edinburgh Sensors (for CO₂) (ii) two specific electrochemical sensors (both from City Technology) for SO₂ and H₂S, and (iii) one KVM3/5 Galltec-Mela sensor for temperature (T) and relative humidity (Rh), which were converted into in-plume H₂O concentrations⁹⁸ using the Arden Buck equation. The most recent configuration of the Multi-GAS, including details on sensor’s accuracy and repeatability, calibration ranges and gas standards, is described in ref. ¹⁰⁰). Our Multi-GAS version, assembled at UniPa, uses a Campbell CR6 datalogger to log simultaneously at 1 Hz the signals outputted from all sensors. We calculate ratios between SO₂ and other volatiles in the fumarolic plumes of individual fumaroles (listed in Supplementary Table 2) by post-processing the acquired dataset using a routine processing protocol implemented in the Ratiocalc software⁹⁹. To this aim, sub-sets of data obtained at individual fumaroles are initially selected, and gas ratios are then derived from the gradients of the best-fit regression lines in gas (either H₂O, CO₂ or H₂S) vs. SO₂ scatter plots. The obtained ratios are listed in Supplementary Table 2. At the high gas concentrations observed at La Fossa (tens to hundreds of ppm SO₂; Supplementary Fig. 2), error on the ratios is typically <10%. This traditional processing routine⁹⁷ is complemented, in the case of the 16–18 November 2021 Multi-GAS walking traverse results, by a supplementary post-processing methodology¹⁰¹ that is aimed at resolving (and mapping; Fig. 3) the spatial heterogeneity of the fumarolic field with high temporal/spatial resolution. This processing routine explores the field acquired time series of GPS-localised gas concentrations, searching for periods of high temporal coherence between CO₂ and SO₂ concentrations in 200 s long moving windows. Within each of the sequentially scanned 200 s long temporal intervals, the routine outputs a CO₂/SO₂ ratio if (i) the Pearson correlation is ≥0.7, (ii) the SO₂ concentration threshold of 20 ppm is reached. The obtained GPS-localised CO₂/SO₂ ratios are interpolated with multi-level B-splines¹⁰² to obtain the contour map of Fig. 3. This figure explores the spatial heterogeneity of the fumarolic field at finer spatial resolution, and is thus complementary to the individual fumarole results listed in Supplementary Table 2.

Melt inclusions

The pre-eruptive volatile content of La Fossa magmas is inferred from information derived from analysis of silicate melt inclusions and by the application of mineral-melt equilibria. We rely on an extensive dataset (reported in Supplementary Table 3) assembled from collection of published major and volatile elements melt inclusion analyses^{31,61,62,103,104,105,106}. The compilation is complemented with unpublished and new results obtained in this work (data all listed in Supplementary Table 3).

Silicate melt inclusions (MIs) are hosted in phenocrysts picked from pumiceous lapilli of pyroclastic deposits, which guarantee quenching of the trapped melt upon eruption. The sampled eruptions span the last 1000 years of activity of La Fossa and Vulcanello (Vulcanello 1, Palizzi eruptive period, Lower Pietre Cotte, 1888–90 eruption; ref. ³⁰; Eruptive Epoch 8 in ref. ²⁸). The dataset covers the entire basalt-to-rhyolite range (Fig. 7a). Since only few basaltic melt inclusions have been found in recent La Fossa products (in Fo90 olivine and diopside antecrysts in 1888–90 eruption⁶²), the melt inclusions of La Sommata (about 42 ky; Eruptive Epoch 5 in ref. ²⁸) are taken as proxy for the basalts feeding the La Fossa plumbing system, based on the similar composition and the petrogenetic reconstructions⁶². Available data from other Vulcano eruptions (Punta Luccia shoshonite, Mastro Minico scoriae latite) of the last 50 ky are also reported^28,30.

The Vulcano Island basalts have measured H₂O contents of 3–4 wt%, and contain 1500–3000 ppm S (Fig. 7), and 2300–4000 ppm Cl; an H₂O content of up to 6 wt% at the time of MI trapping has been postulated considering H diffusion through the host olivine¹⁰⁷. Our B1-B6 model runs are therefore initialized at 3.9 wt.% H₂O and 2500 ppm S (Supplementary Table 4; Figs. 6, 7). A marked decrease in H₂O and S, and a moderate Cl decrease, is registered at the transition from basalt to shoshonite composition (Supplementary Table 3 and Fig. 7). The shoshonitic and latitic MIs exhibit ≼3 wt% H₂O, and the majority of them have H₂O contents between 1 and 2 wt.% (Fig. 7b). S is typically at ≤1000 ppm and Cl at 2000–3000 ppm. The K-trachyte melt inclusions have H₂O contents of 1.8–2.5 wt%, S mostly below 500 ppm, and Cl in the 3000–4000 ppm range; the maximum water content of K-trachytes is confirmed by the application of plagioclase-liquid hygrometer (up to 2.54 ± 0.57 wt%, ref. ¹⁰⁸). Our L1 and T1 model runs are designed as to cover the compositional range of latites and K-trachytes, respectively, and they are initialized at 2 and 2.5 wt% H₂O and 1000 and 400 ppm S, respectively. Rhyolitic melt inclusions are only found in the rhyolitic pumiceous lapilli of the Palizzi period and show maximum H₂O contents of 2.5–3 wt% (1.99 ± 0.30 according to the application of pl-lq hygrometer¹⁰⁸), S < 100 ppm and Cl up to 5000–7000 ppm. Our R1 model run is therefore initialized at 2.5 wt.% H₂O and 80 ppm S. CO₂ is not detected during FTIR analysis of MIs of any composition (50 ppm is the FTIR detection limit), and low CO₂ values are detected by SIMS in La Sommata basaltic melt inclusions (<184 ppm, ref. ¹⁰⁵). No CO₂ partitioning into the shrinkage bubble is found using Raman spectroscopy, suggesting CO₂ loss before MI entrapment. Our model runs are therefore initialized at 1 to 2 wt.% CO₂ using information¹⁰⁹ from nearby Stromboli volcano. The rhyolitic run R1 is initialized at lower CO₂ contents (0.1 wt.%), but the effect of higher CO₂ (1 wt.%) is also explored (in run R2, see Supplementary Table 4 and Supplementary Fig. 4).

Literature information is complemented with novel melt inclusion information for selected compositions (rhyolitic melt inclusions) that have received little or no attention in previous studies (Supplementary Table 3). Major elements, S and Cl were measured by EMPA at the LaMA Laboratory of the Department of Earth Sciences (University) and CNR-IGG of Florence, with a JEOL Superprobe JXA-8230. The instrument is equipped with 5 wavelength dispersion spectrometers (WDS) operating with different analyser crystals and different detectors (gas flow and sealed Xe-type). The working conditions for glasses are 15 kV of accelerating voltage, with variable beam current intensity depending on the different matrix material under analysis and 5 nA beam current. The beam diameter varies from 5 μm to 10 μm for glasses to reduce the diffusion effect on alkali. Variable counting times are set for major and minor elements: 15 s on peak and 7 s on both backgrounds for major elements; 30–40 s on peak and 15–20 s on both backgrounds for minor elements, except for Na that is measured for 10 s on peak and 5 on both backgrounds. The acquisition order starts with Si, Na, K, Cl and Fe, followed by Mg, Ca, Mn, F, S and Al, Ti and P and, eventually Sr and Ba. A selection of natural phases are used as primary standards for the elemental calibration (Astimex Albite for Si and Na, Plagioclase for Al, Olivine for Mg, Diopside for Ca, Sanidine for K, Bustamite for Mn, Apatite for P, Celestine for Sr and S, Tugtupite for Cl, Barite for Ba, Fluorite for F and Smithsonian ilmenite for Ti and Fe). H₂O and CO₂ contents are measured by FTIR using a Nicolet Magna-IR560 spectrometer coupled with a NicPlan infrared microscope (at Dipartimento di Scienze della Terra, Pisa, Italy). At least three spectra for different areas were collected for each MI. The concentrations (c, wt.%) were calculated according to the Beer–Lambert law: c = 100 × A × M/[ε × ρ × d], where A (absorbance) is the height of the IR absorption peak, M is the molar mass (g·mol⁻¹, 18.02 for total H₂O, molecular H₂O, and OH and 44.00 for CO₃²⁻), ε is molar absorptivity (l · mol⁻¹ cm⁻¹), ρ is the density of glass (g·cm⁻³), and d is the thickness of the inclusion. Total dissolved H₂O is measured from the intensity of the broad, asymmetric peak at 3550 cm⁻¹, which corresponds to the fundamental OH⁻ stretching vibration¹¹⁰. We did not use the peaks at 4500 and 5200 cm⁻¹ because they were strongly superimposed by interferences fringes, which make interpretation of the spectra difficult. Absorbance, A, of the band of interest is calculated from the spectra by averaging at least two to three different height measurements of the absorption band peak. Glass density (ρ) is estimated based on MI composition at 2400 g/cm³. Thickness (d) of each inclusion is calculated by averaging results of three distinct measurements made using a Mitutoyo digital micrometer in three different points before the analysis, with an error of 2–3 μm. We use a molar absorptivity of 56 for total dissolved H₂O at the 3550 cm⁻¹ band. No C peak is observable above the background, indicating that CO₂ is always below the detection limit (~50 ppm).

Magma degassing model calculations

We use the volatile (C-O-H-S system) saturation model of Moretti et al.⁵⁷ to calculate (i) dissolved volatile contents (Fig. 7) and (ii) equilibrium composition of magmatic fluids (Figs. 4, 6, 8; Supplementary Fig. 4) at T-P-X conditions relevant to Vulcano Island magmas. Full details and model theoretical background and functionalities can be found in refs. ^{56,57,111,112,113}. The model input conditions are summarised in Supplementary Table 4. Model runs are initialised using representative parental melt compositions (major and volatile elements) determined from analysis of the Vulcano Island melt inclusion record (Supplementary Table 3). Runs are performed by simulating magma decompression from an initial pressure (see table) down to atmospheric conditions (0.1 MPa), under isothermal conditions, and in either closed-system (i.e., gas and melt continuously re-equilibrating during the decompression path) or open-system (gas-melt separation at each decompression step) conditions.