The SO₂ flux results we report on are derived from processing of images taken by the permanent “La Montagnola” (MON) UV camera system (Figure 1). Hardware/software characteristics and the system’s functionalities have been described (D’Aleo et al., 2016; D’Aleo et al., 2019; Delle Donne et al., 2019). In brief, the system employs two JAI CM-140GE-UV cameras sensitive to UV-radiation, fitted with two different band-pass optical filters (Full Width at Half Maximum, FWHM, of 10 nm) with central wavelengths of 310 and 330 nm, respectively. The two co-aligned cameras are programmed to record daily sequences of images (520×676 pixel at 10-bit) at 0.5 Hz rate during optimal illumination conditions (e.g., from 7.30 a.m. to 12.30 p.m. UTC time). Image acquisition is automated by a Fujitsu RX100 Workstation.

The daily images are post-processed by a standard protocol (Kantzas et al., 2010; Delle Donne et al., 2019), in which sets of simultaneously acquired images (from the two cameras) are combined into single absorbance images, applying the Beer–Lambert law. This allows quantification of the relative absorption of scattered UV sunlight by in-plume volcanic SO₂ over the cameras’ field of view (FoV). Sequences of SO₂ absorption images are then converted into SO₂ column density images using the proportionality ratio between absorbance and SO₂ column density, derived from an Ocean Optics USB 2000+ spectrometer pointing in a defined subregion of an absorbance image (Figure 1). Finally, SO₂ column densities are integrated along a profile stretching over the median portion of a rectangular integration sub-area, positioned over the crater terrace (Figure 1). The obtained SO₂ integrated column amounts (kg⋅m^-1), once multiplied by plume speed (m⋅s^-1; derived by applying an optical flow algorithm that tracks gas fronts in consecutive image frames; Delle Donne et al., 2019), yield a 0.5 Hz time step for SO₂ flux (kg⋅s^-1) measurement, examples of which are shown in Figures 2A–C.

The entire operation sequence described above is completed using a MATLAB code that can deliver SO₂ flux time-series in near real-time (Delle Donne et al., 2019). The results reported here (Figure 2; Figure 3) are, however, based on post hoc manual operator check and validation of these automatically processed data. In particular, the image sequences are visually inspected, and any frame of limited (or no) visibility (cloud coverage), or any frame (e.g., C4 and C6 in Figure 2 discussed below) with visible ash in the plume (as indicated by a co-exposed visible camera), are omitted from calculation of SO₂ flux daily averages (illustrated in Figure 3B). This is because the presence of ash in the plume is known to significantly interfere (suppress) the in-plume SO₂ signal (e.g., frames C3 and C5 in Figure 2 discussed below), resulting in artificially low SO₂ fluxes. The implication is that our observations severely underestimate syn-eruptive SO₂ emissions during a paroxysm’s climax, which can be more fully characterized from spaceborne instruments (Corradini et al., 2021). However, as paroxysmal episodes are relatively short-lived (minutes to hours), the under-represented paroxysmal SO₂ contribution is relatively marginal (D’Aleo et al., 2019) if long-term (days to years) SO₂ degassing trends are of interest, as it is the case here.

Artificial suppression of the real SO₂ flux signal by ash is not the only disadvantage of UV Cameras compared with other more established measurement techniques, e.g., Differential Optical Absorption Spectroscopy (DOAS; Platt and Stutz, 2008; Galle et al., 2003; Galle et al., 2010). Targeting the proximal–and therefore often most optical opaque–portion of a plume (Figure 2), UV Cameras are more prone to errors caused by in-plume multiple scattering caused by aerosols and water droplets (Kern et al., 2010). We minimize this effect by selecting (from operator inspection of visible camera images) frames with visually transparent (uncondensed) plumes. More challenging is to account for the effect of radiation scattering (light dilution) in the line of sight between the volcanic plume and the SO₂ camera (Kern et al., 2010). This can attenuate the SO₂ signal, leading to significant underestimation of SO₂ column amounts, depending on the distance between the observation site and the plume. To estimate random and systematic errors in SO₂ column amount retrievals caused by light dilution, we apply the methodology of Campion et al. (2015) assuming an observation distance of 3 km (the distance between MON and Etna’s NSEC crater). For the Etna’s specific conditions, the methodology of Campion et al. (2015) allows derivation of a linear regression between measured and theoretical true scan column densities with a correlation coefficient (R ²) of 0.99, and a slope of ∼0.68. From this regression, we estimate a potential systematic underestimation of the true SO₂ column densities of ∼32%, and an average random error associated with each measurement of ∼20%.

The UV Camera FoV (Figure 1) also limits the range of detectable gas emissions to those produced by the South East Crater (SEC) complex (that includes the NSEC) and from the southern portion (Bocca Nuova) of the Central Craters (D’Aleo et al., 2019) while leaving undetected the emissions of the North East Crater (NEC) and Voragine (the VOR emissions can only be measured during vigorous explosive activity that produce high, sustained plumes; D’Aleo et al., 2019). Therefore, our SO₂ flux measurements are unrepresentative of the bulk (total) plume emissions that, on Etna, are regularly monitored by more distant observations with scanning UV spectrometers (Salerno et al., 2009).

The daily averaged SO₂ fluxes are listed in Supplementary Table S1. The chronology of volcanic activity (see below), and the complete SO₂ flux time-series, are illustrated in Figures 3A, B.

To more fully characterize eruptive activity, the SO₂ flux record is complemented with seismic tremor data and thermal records from both satellite and ground-based instruments (Figure 2 and Figure 3).

Volcanic tremor is the most distinctive seismic signal recorded on active volcanoes. Generally associated with physical processes involving interactions between magmatic fluids and surrounding rocks inside the volcano (Chouet and Matoza, 2013), according to prevailing theory volcanic tremor originates from acoustic resonances of fluid-filled cracks triggered by pressure transients (Chouet, 1996). Volcanic tremor is often associated with the upward migration of magma and gases in the plumbing system (McNutt and Nishimura, 2008), and provides useful information on volcanic activity evolution (e.g., Almendros et al., 2001; Patané et al., 2013). At Mt. Etna, these seismovolcanic signals, and their variation over time, have turned out to be strictly related to eruptive activity (e.g., Patanè et al., 2008; Viccaro et al., 2014), and have provided clues on pressurization or depressurization phenomena within the plumbing system (Cannata et al., 2015). Furthermore, the temporal evolution of volcanic tremor amplitude has been found to strongly correlate with the intensity of explosive activity (Andronico et al., 2021) (compare Figures 2A–C). In fact, tremor amplitude gradually increases during the activity ramp-up that leads to a paroxysmal episode, peaking at the eruption climax, and then decreasing rapidly to return to a normal levels at the end of a paroxysmal episode (Alparone et al., 2003; Cannata et al., 2008) (Figure 2C). In order to track volcanic tremor amplitude changes, we calculate the Root-Mean-Square (RMS) of the seismic signal recorded by the vertical component of Serra La Nave seismic station (ESLN), located on the upper southern flank (Figure 1). This station is equipped with a broadband (40 s cut-off period), 3-component Trillium Nanometrics™ seismometer, acquiring at a sampling rate of 100 Hz. The RMS time series is calculated on 1 min long time windows of seismic signal filtered in the band 0.5–5.5 Hz. This frequency band includes most of the volcanic tremor energy (e.g., Cannata et al., 2008). The entire RMS time-series is illustrated in Figure 3C. The seismic station located at Cratere del Piano (ECPN), closer to the crater, exhibits a similar trend (Supplementary Figure S1).

Etna’s paroxysmal episodes are also characterised by thermal observations from both ground and space. In order to derive high-rate information on thermal activity during a paroxysmal episode (data in Figure 2C), we use ground-based long-wave infrared (LWIR) records derived from a FLIR Thermovision A40M infrared camera located in Nicolosi, 15 km south of the summit craters (ENT). The camera is installed on the roof of a school in front of the INGV headquarters in Nicolosi. Thermal images are converted on board by the camera into a single output value that corresponds to the maximum temperature (Figure 2C) recorded in a region of interest, centred on the zone potentially affected by eruptive activity. This temperature value is transmitted at 1 Hz, through a WiFi system, to the acquisition centre in Nicolosi, where thermal images are processed in real time by the NewSaraterm software (Andò and Pecora, 2006), developed with LabVIEW^TM 8.0, to detect ash-rich eruptive columns, explosive and effusive activity. This software calculates, with a median filter, the daily thermal offset and extracts, for each measure, the difference between the detected value and the dynamic thermal offset. NewSaraterm also sends messages to the INGV-OE Control Room in Catania when a set of thresholds is exceeded, to alert on-duty staff for early warning procedures. The maximum apparent temperatures acquired by this camera are also published in real time on the INGV-OE Intranet page. The system has demonstrated its ability to recognize thermal anomalies over background. The temperatures displayed in the video frames are much lower than the true temperatures at the target owing to distance between the thermal camera site and the target (the summit craters), and the effect of environmental factors such as air humidity, pressure, air temperature, and the presence of protective lenses on the camera (Ball and Pinkerton, 2006; Sawyer and Burton, 2006; Spampinato et al., 2011). The temperature values illustrated in Figure 2C are thus apparent, not absolute, as they have not been corrected for the aforementioned factors. The chosen temperature range on the infrared camera is from —40°C to 120°C, and temperatures exceeding 130°C are not shown in Figure 2C.

Longer-term fluctuations in radiant heat (Figure 3D) are characterised using radiant heat flux derived from space-borne images acquired by the Spinning Enhanced Visible and InfraRed Imager (SEVIRI) on board Meteosat Second Generation (MSG) satellite. MSG-SEVIRI multispectral images are retrieved via the CL-HOTSAT system, a new version of the HOTSAT system (Ganci et al., 2011; Ganci et al., 2016; Ganci et al., 2018) implemented in the OpenCL standard to obtain near-real-time monitoring of volcano-related thermal activity. The SEVIRI sensor is especially useful to observe and track rapid eruptive events such as Etna’s paroxysmal explosions (Calvari et al., 2021; Calvari et al., 2022), since it allows the acquisition of images in the thermal and medium infrared region of the electromagnetic spectrum at very high temporal resolution (every 15 min; up to 5 min using the rapid scanning service at medium latitudes). The CL-HOTSAT system developed at INGV-OE allows for automatic detection of thermal anomalies and quantification of thermal activity through calculation of the radiant heat flux. The radiant heat flux time-series in 2020–2021 is illustrated in Figure 3D.

In the case of effusive eruptions, the radiant heat flux can be converted into a time average discharge rate (TADR) (Harris et al., 1998) assuming a thermal balance between the heat within the erupted flow volume and that lost from it. This balance implies that i) the volumetric flux is proportional to flow area, and ii) lava surfaces cool exponentially with time and distance from the source (Harris and Baloga, 2009). These assumptions cannot be verified in the case of transient events such as Etnean paroxysmal events, during which steady-state, and hence thermal equilibrium, cannot be reached (Garel et al., 2013). However, for these short and transient eruptive events, a satellite-derived volume estimation can be derived from modelling the radiant heat flux curve measured during the cooling-dominated phase of activity (Ganci et al., 2012; Ganci et al., 2013). Our satellite-derived radiant heat flux (Figure 3D) is converted into magma volumetric fluxes (Figure 3D) by modeling the time t) required for lava flow to cool following equation: t = T h − T a ε σ ( T s u r f 4 − T a 4 ) + h c T s u r f − T a − R c R r R c + R r k 1000 1 2 λ 2 1 κ (1) where T_surf is the surface temperature for a stagnant, stable, cooling lava surface, T_h is the initial temperature, T_a is ambient temperature, ε is emissivity, σ is the Stefan-Boltzmann constant, h_c is the convective heat transfer coefficient, k is thermal conductivity, λ is a dimensionless scaling value, κ is thermal diffusivity, and R_c and R_r are thermal resistance terms for convection and radiation, respectively. Converting the surface temperature in radiant heat flux, hence including the radiating area, and fixing a minimum and maximum thickness expected for these short-lived events, an estimation of the emplaced volume can be obtained. Part of the thermal dataset used here is also discussed in Ganci et al. (2023).

A chronology of volcanic activity during our observational interval (January 12,020 to December 31 2021) is illustrated in Figure 3A, as derived from INGV-OE monitoring reports (INGV-OE, 2020). As evident, volcanic activity has taken a variety of forms during the study period, with the four summit craters being intermittently (and at time all concurrently) active. As our primary goal is to analyze degassing dynamics at NSEC (the most active crater, but also the one that can be most successfully observed from our UV Camera position, Figure 1), we classify volcanic activity based on the characteristics exhibited by activity at this specific crater. Based on this, we categorize our observational interval into seven categories (activity periods), as discussed below.

The Precursory Phase PPa (January-April 2020) is distinguished for the lack of any visible eruptive activity at the NSEC (Figure 3A). Volcanic activity is therefore confined at NEC and VOR that exhibit some relatively continuous intra-crateric explosive (ash expulsion and mild Strombolian events) and effusive (at VOR) activity.

The Precursory Phase (PPb) (May 2020 to 12 December 2020) marks resumption of volcanic activity at the NSEC. Activity progressively intensifies (becomes increasingly more frequent; Figure 3A) at the NSEC during PPb, but remains limited to mild Strombolian activity and intermittent ash expulsion events. No paroxysmal episode takes place in this phase. Sub-period PPb’, the last 2 months of phase PPb (Figure 3A), is identified based on a visible acceleration in degassing activity at NSEC (based on the SO₂ flux record; see below).

Onset of paroxysmal activity on 13 December marks onset of PS1, the first of two especially powerful paroxysmal sequences (dark orange in Figure 3A) occurring in 2020–2021 (INGV-OE, 2020; Andronico et al., 2021; Bonaccorso et al., 2021; Calvari et al., 2021; De Gori et al., 2021; Marchese et al., 2021; Corsaro and Miraglia, 2022). PS1 is opened by 4 “transitional” (following the terminology of Andronico et al. (2021): “strong Strombolian activity alternating with short periods of lava fountains”) eruptive episodes on 13, 21 and 22 December 2020, and on 18 January 2021 (Andronico et al., 2021; Bonaccorso et al., 2021; Corsaro and Miraglia, 2022) (Figure 3A). Since 16 February 2021, paroxysmal episodes become more rhythmic and closely spaced (Figure 3A). In total, PS1 includes 21 paroxysms, 17 of which are powerful lava fountaining episodes between February 16 and 1 April 2021 (Figure 3A). Although all episodes occur at NSEC, PS1 is notable for that all summit craters exhibit some milder intra-crateric explosive activity (Figure 3A).

The PS1 paroxysmal sequence closes abruptly with a final paroxysm on 1 April 2020, and is followed by circa 2 months of relative quiescence (no paroxysmal activity) in Post Paroxysmal Phase 1 (PPP1).

Paroxysmal activity resumes again at NSEC on 19 May 2021, marking the onset of a second paroxysmal sequence (PS2) that terminates on October 23. PS2 includes 34 episodes in total (INGV-OE, 2020) (Figure 3A), and is followed by a period of relative quiescence (Post Paroxysmal Phase 2; PPP2) until the end of the year (Figure 3A).

Figure 2 shows examples of UV-camera derived, 0.5 Hz SO₂ flux time-series obtained in three distinct days, taken respectively as representative of a) passive degassing during a quiescence period, b) passive degassing during a paroxysmal sequence (e.g., during a repose interval in-between two consecutive paroxysmal episodes), and c) a temporal interval encompassing a paroxysmal episode. These examples, although not exhaustive, illustrate the variety of Etna’s degassing modes in 2020–2021, and clearly identify a progressive SO₂ flux increase from degassing regime a) to c). During quiescent periods (e.g., Figure 2A), the SO₂ flux exhibits the typical “pulsed” degassing behavior seen in previous work (Tamburello et al., 2013), in which nearly periodic SO₂ flux fluctuations repeat at periods of tens to hundreds of seconds. These are interpreted to reflect either atmospheric/dispersion processes, or the rhythmic release of sequences of concentrated SO₂ clouds (or puffs) from the open vents (Tamburello et al., 2013). In this specific example, the daily averaged SO₂ flux is typically at <1000 Mg/day level, and seismic tremor amplitude is at low levels, too (Figure 2A; Figure 3). Degassing regime b) (passive degassing during repose intervals between two consecutive paroxysmal episodes) is similar to regime a), but the SO₂ flux fluctuates at around higher (>>1000 Mg/day) time-averaged levels, and seismic tremor amplitude is also more intense (Figure 2B; Figure 3). The SO₂ flux reaches its highest levels (>5000 Mg/day) in the run-up to a paroxysmal episode, during which seismic tremor and ground-based thermal activity also escalate (Figure 2C). Seismic tremor and thermal activity peak during the paroxysmal climax (Figure 2C), during which our SO₂ time-series are abnormally low. This apparent sin-paroxysm SO₂ drop is a technical artifact due to the optically opaque nature of the ash-laden eruptive cloud (Figure 2C) that precludes quantitative SO₂ detection (D’Aleo et al., 2019).

The ash-free (cfr. 2.1) SO₂ flux results acquired during each observation day (examples in Figure 2) are averaged to derive the daily mean SO₂ fluxes illustrated in Figure 3B. The plot clearly highlights large temporal oscillations in the daily averaged SO₂ flux. In the relatively calm PPa phase, the SO₂ flux averages at 1977 ± 1531 Mg/day (Figure 3B), but nearly doubles during the subsequent PPb phase of intensifying activity at the NSEC (mean, 3351 ± 2515 Mg/day) (here and below, we quote mean ±1 standard deviation). Onset of PS1 on December 13 marks an abrupt SO₂ flux increase; the mean flux for this period is 4836 ± 3667. A large inter-daily variability is observed (range, 880–19850 Mg/day) during PS1, and the peak SO₂ flux in our time-series (19850 Mg/day) is observed on December 14 (a few hours after the first paroxysmal episode on December 13). End of PS1 on April 1 marks an abrupt SO₂ emission drop, and the SO₂ flux remains at low levels of ∼2600 Mg/day (range, 580–8500 Mg/day) during the entire PPP1, to manifestly increase again during PS2 (mean, 4290 Mg/day; range 760–12160 Mg/day). The lowest fluxes are recorded in the PPP2 (mean, 1484 ± 1059 Mg/day).

Tremor amplitude (Figure 3C) exhibits long-term (weeks-to months-long) fluctuations (as exemplified by the 20 days moving average; red curve), over which sharp and transient (hours) increases are superimposed that mark the repeated Strombolian and lava fountaining episodes in 2020–2021. The two paroxysmal sequences (PS1 and PS2) are therefore well evident in the seismic tremor record as clusters of repeated, abrupt tremor amplitude escalations during intense explosive intensity. The longer-term trend (red curve) identifies a progressive (thought irregular), mild tremor amplitude increase from PPa to PPb. This tremor amplitude escalation culminates in phase PS1, during which seismic tremor amplitude remains sustainedly high (background tremor, red curve, peaking on January 27). Seismic tremor amplitude then abruptly drops during phase PPP1, to resume again in phase PS2, though at reduced levels relative to phase PS1. Tremor amplitude is overall at low levels in phase PPP2, with some transient increases in December 2021 (Figure 3C).

Thermal anomalies are observed in SEVIRI images in correspondence to almost all the paroxysmal eruptions (Figure 3D). Images recorded during the 17 March, 27 May and 16 June paroxysmal episodes are affected by consistent cloud coverage that hampers derivation of erupted volumes. The paroxysmal episodes of the PS1 paroxysmal sequence are the most energetic, with heat flux peaks above 20 GW (Figure 3D). These convert into volumes of up to ∼2.5 ×10⁶ m³ per episode (Figure 3D), providing an average output rate of about 7 m³/s. The second paroxysmal sequence (PS2) is instead characterized by prevalently lower heat fluxes (<20 GW), and by an average output rate of 2.8 m³/s. The median volume per episode is also lower (<1.6·10⁶ m³), and <1·10⁶ m³ until 22 June (Figure 3D). The last two episodes (21 September and 23 October) are exceptional for their heat flux (>20 GW) and estimated volumes (of ∼1.5·10⁶ m³) (Figure 3D).

Recent advances in geophysical (Alparone et al., 2007; Aloisi et al., 2011; 2018; Patanè et al., 2013; Bonaccorso and Calvari, 2017; Cannata et al., 2018; De Gori et al., 2021) and geochemical (Allard et al., 2005; Aiuppa et al., 2010; La Spina et al., 2015; Paonita et al., 2021) monitoring, combined with careful volcanological, petrological and remote sensing observations (Corsaro et al., 2017; Calvari et al., 2018), have led to enormous progress in understanding, interpreting and modelling the drivers of Etna’s paroxysmal sequences. Consensus has been reached that the transition from ordinary to paroxysmal activity at NSEC is a response to escalating supply of volatile-rich, mafic magma into a vertically elongated, shallow (∼2–4 km depth below the summit, bvs) magma storage zone (Figure 5A). Mixing of this mafic magma with more evolved, shallow resident magma (Kahl et al., 2015; Corsaro et al., 2017; Giuffrida and Viccaro, 2017) would cause gas-driven pressure build-up over timescales of days to months (Giuffrida and Viccaro, 2017), and ultimately eruption (Bonaccorso et al., 2013; Patanè et al., 2013; Corsaro et al., 2017; Cannata et al., 2018; Bonaccorso et al., 2021). For example, the PS1 sequence in late 2020 to early 2021 has recently been proposed to represent the culmination of a magma resupply event of the ∼2–4 km deep magma storage zone, started in late 2020, as also indicated by a seismic swarm on 19 December 2020 (Bonaccorso et al., 2021). The escalating tremor amplitude since December 13, as seen in our results (Figure 3C), is entirely consistent with this hypothesis.

Indirect evidence for recharge/over-pressuring of the shallow Etna’s magma ponding zone comes from the petrology/geochemistry of erupted products (e.g., Pompilio et al., 2017), from modelling of seismic and geodetic signals (e.g., Bonaccorso et al., 2021), and from the chemistry of volcanic gases released prior to (Aiuppa et al., 2010; Paonita et al., 2021) and during (Allard et al., 2005) the lava fountaining episodes. Inversion of geodetic data (Palano et al., 2008; Aloisi et al., 2011; Palano et al., 2017; Aloisi et al., 2018; Cannata et al., 2018), and petrological information (Métrich et al., 2004; Corsaro et al., 2017; Giuffrida and Viccaro, 2017; Cannata et al., 2018), have been used to suggest the mafic magma feeding paroxysmal sequences is sourced by an intermediate (at ∼ 8–11 km depth below sea level, bsv) magma storage zone (Figure 5A), part of a vertically elongated, complex transcrustal magmatic feeding systems extending down to the Moho (Giuffrida and Viccaro, 2017) (Figure 5A). However, a fundamental challenge remains the identification of the magma input rate threshold (if any) at which the volcano behaviour transitions to a paroxysmal sequence.

SO₂ fluxes are especially useful in this regard as they can be converted into magma input (degassing) rates with the petrological method based on the equation: Q i n = ϕ S O 2 2 ∆ X S ∙ 1 − x ∙ ρ m ∙ 86.4 (2) where Q_in is the magma input rate (in m³/s) above the SO₂ exsolution level (the level at which effective sulphur S degassing from melt into a fluid phase starts), ΦSO₂ is the SO₂ flux (in Mg/day), ΔX_S is S loss upon degassing (as inferred from the difference between parental melt and residual—post-eruption—S contents in Etna’s magmas; Figure 5B), and x and ρ_m are melt crystallinity and density (in kg/m³), respectively. In brief, the equation calculates the volumetric magma input rate required to justify the surface-measured SO₂ output, based on knowledge of initial (pre-degassing) and final (post surface emplacement) S content in melt, and has extensively been used at Etna to constrain degassing magma budgets over decadal (Allard, 1997; Allard et al., 2006), yearly (Coppola et al., 2019), or even daily (D’Aleo et al., 2019) timescales. Estimates point to a general decline of time-averaged Q_in values, from 5.4 to 8.9 m³/s in the 1975–1995 period (Allard, 1997), to ∼3 m³/s in 1995–2005 (Allard et al., 2006), and ∼1.6 m³/s in 2005–2010 (Coppola et al., 2019), reflecting an SO₂ flux drop following the 2001 Etna’s flank eruption (Catalbiano et al., 2003). The common conclusion is, however, that such SO₂-derived, time-averaged magma input rates are larger (1.8–4 times higher) than the mass eruption rate, implying long-term endogenous volcano growth (Allard, 1997) and deformation (Coppola et al., 2019).

Application of Equation (2) requires a model that suitably describes S behaviour during magmatic degassing of Etna’s alkali basalts. This model, starting from knowledge of the initial S content in Etna’s parental magmas, should be able to predict the fraction of S lost to vapor phase at each step of the magma decompression path, and up to surface pressure (0.1 MPa) (Allard, 1997). Deeply trapped, primitive melt inclusions are typically used to infer the initial (pre-degassing) parental melt S contents (Wallace and Edmonds, 2011). At Etna, melt inclusion results for the 2001 and 2002–2003 eruptions (Métrich et al., 2004; Spilliaert et al., 2006) indicate a parental melt S content of ∼3000 ppm (0.3 wt%), and a systematic decrease of dissolved S with pressure decrease (with S loss to vapour starting at pressures ≤150 MPa; Figure 5B). This 3000 ppm estimate, and 95% S extraction from melt upon degassing at ≤150 MPa (e.g., ΔS = 0.00285 in Eq. (2)), have been used in previous magma input rate estimates (Allard, 1997; Allard et al., 2006; Coppola et al., 2019). However, more recent work (Gennaro et al., 2019) has found more variable parental melt S contents in magmas, and less systematic (than found for the 2002–2003 melt inclusions; Spilliaert et al., 2006) pressure-dependence of melt S contents (Figure 5B). The authors suggest that variable redox conditions exert a strong control on dissolved S in Etna’s magmas.

In view of the above, we use the model of Moretti et al. (2003) to simulate S behaviour during de-compressional degassing of Etna’s magmas, using two distinct scenarios for the initial S contents (2001 runs: 3000 ppm, same run conditions as in Aiuppa et al., 2007; FS runs: 2200 ppm, melt composition as in Gennaro et al., 2019). We also explore the role of magma redox by running degassing model simulations over a range of redox conditions (from −2 to +2 log units below/above the Nickel-Nickel Oxide buffer, NNO) (see Table 1). Our model results (Figure 5B) output pressure-dependent model dissolved S contents matching the compositional range of Etna’s melt inclusions (Figure 5B). Importantly our model calculations demonstrate that, independently on model input conditions, S loss to vapor is negligible for pressures ≥200 MPa, and only becomes significant for pressures ≤150 MPa (Figure 5C). Therefore, in spite the S degassing path is manifestly redox-dependent (with more reducing conditions - NNO-2—leading to some larger “deeper/earlier” S degassing), it remains that S loss to vapor is minimal in the deep Etna’s plumbing systems (e.g., at conditions deeper than the 8–11 km bsv storage zone of Figure 5A). This makes of the surface measured SO₂ flux a key tracer of the rates of magma supply, from the intermediate magma storage zone to the shallow portions of the Etna’s plumbing system, where paroxysmal sequences are initiated (Figure 5).

We caution that our degassing models offer a simplified view of real degassing conditions of Etna’s magma. For example, our models are run in isothermal conditions, and therefore do not consider the effect of crystallization upon cooling, and the consequent role crystals may play in favouring bubble nucleation (Woitischek et al., 2020). Likewise, the model assumes gas-melt equilibrium at each step of the decompression path, and therefore cannot account for any disequilibrium generated during gas bubble expansion within the melt (Oppenheimer et al., 2018). This latter effect, however, mostly impacts gas redox and sulphur speciation (e.g., the SO₂/H₂S ratio) rather than total S degassing behaviour. Tests made by initialising runs at different initial volatile conditions indicate that varying the parental melt H₂O content leads to some minor shift in the derived pressures of sulphur degassing initiation, but does not alter our overall conclusion of a shallow (≤150 MPa) condition of sulphur loss to the gas phase. Therefore, although simplified, our models offer a useful representation of the portion of the Etna’s plumbing system in which sulphur degassing actively takes place.

Stemming from the above, we use Equation 2 to convert our daily averaged SO₂ fluxes (ΦSO_2; data from Figure 3B) into daily averaged magma input rates (Q_in) by taking into account the possible full spectrum of ΔS and crystallinity x values in Etna’s magmas (Figures 4B, C).

We illustrate the possible range of magma input rates Q_in as a white shaded area in Figure 4C (in the diagram, the Q_in values are expressed in the form of 7 days mobile averages). This area is limited by an upper range estimate (obtained using ΔS = 0.00205 and crystallinity x = 0.3) and a lower range estimate (obtained fixing ΔS at 0.00285 and crystallinity x at 0). These upper and lower range scenarios correspond to initial S contents of 2200 and 3000 ppm respectively, and 95% S extraction from melt upon degassing (see above). Density is fixed at 2700 kg/m³ in all calculations. The blue solid line in Figure 4C represents an intermediate scenario, obtained by averaging results of the upper and lower range scenarios. The same upper/lower/intermediate range estimates are illustrated in Figure 4B in the form of cumulative degassing magma volume, e.g., the cumulative volume of degassing magma required to explain the surface SO₂ emissions. This cumulative volume is obtained from the summation of daily Q_in values.

Our results constrain the time-averaged magma input rate Q_in for the entire temporal interval (2020–2021) at 3.6 ± 2.7 m³/s (the mean values quoted here and below are based on the intermediate scenario; blue curve), within the range of previous estimates (Allard, 1997; Allard et al., 2006; Coppola et al., 2019). However, we observe large Q_in temporal fluctuations (Figure 4C), with a progressive increase from PPa (mean: 2.0 m³/s; lower-upper scenario averages: 1.4–2.7 m³/s) to PPb (mean: 3.3 m³/s; lower-upper scenario averages: 2.2–4.3 m³/s) and finally PS1 (mean: 4.7 m³/s; lower-upper scenario averages: 3.5–6.9 m³/s) (Figure 4C). Then, after a sharp drop in PPP1 (mean: 2.7 m³/s; lower-upper scenario averages: 1.8–3.6 m³/s), Q_in escalates again during PS2 (mean: 4.7 m³/s; lower-upper scenario averages: 3.2–6.3 m³/s). The lowest estimated Q_in values are for PPP2 (mean: 1.8 m³/s; lower-upper scenario averages: 1.2–2.4 m³/s) (Figure 4C).

The calculations indicate that the magma input rate is 2-3 higher during a paroxysmal sequence (PS1 and PS2) than during the ensuing Post Paroxysmal phases (PPP1 and PPP2). Also, paroxysmal sequences exhibit on average a factor 30%–45% Q_in increase relative to the precursory Phase (PP). At this point, we focus our attention on the last 2 months of the precursory phase (PPb’; identified in dark orange in Figure 3 and Figure 4), in which the time averaged SO₂ flux increases (to ∼4260 Mg/dayay), and the corresponding Q_in accelerates to values (mean: 4.7 m³/s; lower-upper scenario averages: 3.2–6.2 m³/s) in all similar to those characteristic of PS1 and PS2. Manifestly, the Q_in rump-up pre-dates (by 1–2 months) the onset of PS1.

Ultimately, our results (Figure 4) suggest that onset and duration of a paroxysmal sequence are controlled by the rates of magma transport into the shallow magma plumbing system that feeds the NSEC (Figure 3 and Figure 4). The key aspect to consider is that our SO₂ flux measurements are unrepresentative of the bulk (total) plume, but rather only capture the SO₂ emissions from the South East Crater (SEC) complex (that also comprises the NSEC) and from the southern segment of the Central Craters (D’Aleo et al., 2019). Notably, our SO₂ fluxes for the 2020–2021 paroxysmal sequences (PS1-PS2; Figure 3B) are in the same range of the SO₂ flux time-averages derived by the INGV-OE network of UV scanners (INGV-OE, 2020; Salerno et al., 2022) that sense the distal, bulk (total) plume (Salerno et al., 2009). This fact strongly suggests that the southern portion of the crater terrace (SEC-Bocca Nuova) becomes the Etna’s dominant source of gas (of SO₂, at least) during a paroxysmal sequence. The corollary is that our inferred magma input rates (Figure 4C) are intended to represent supply rates to the SEC-Bocca Nuova feeding system–these can only approach the total magma input rates only when/if activity concentrates at NSEC (e.g., shortly prior, or during, a paroxysmal sequence).

Additional clues can be derived from a comparative analysis of the SO₂ flux record (Figure 3B) with independent geophysical evidence (Figure 3C; Figure 3D).

We infer that, in the 2 years covered by our study, ∼15 Tg of SO₂ have cumulatively been released (Figure 4B). When the intermediate scenario is considered (blue curve labelled DMV), this cumulative SO₂ mass output requires degassing of ∼260·10⁶ m³ of magma (from a minimum of 180·10⁶ m³ in the lower range scenario to a maximum of 354·10⁶ m³ in the upper range scenario). Of these, we infer that roughly half (∼120·10⁶ m³) have sustained degassing during the 266 days interval covering the two paroxysmal sequences (PS1-PS2). Our MSG-SEVIRI results (Figure 3D), owing to the relatively coarse (1–3 km) spatial resolution of the multispectral images used, capture only feeble and inconstant thermal anomalies during open-vent ordinary activity (PPa, PPb, PPP1 and PP2 phases) (Figure 3D). In contrast, very intense radiant heat fluxes are detected by SEVIRI associated to each of the PS1/PS2 paroxysmal episodes (Figure 3D). These radiant heat fluxes, once converted into erupted volumes (Figure 3D), constrain the total (cumulative) magma output (required to explain the satellite-measured thermal emissions) at ∼70·10⁶ m³ (Figure 4B). These calculations imply that, during the 2020–2021 paroxysmal sequence, the output magma volume (∼70·10⁶ m³) corresponds to a relatively large (∼58%) proportion of the magma supplied to (and releasing SO₂ from) the shallow Etna’s magma plumbing system. Such relatively large output/input fraction (0.58) is well above that typically observed during ordinary activity, when a large output vs input unbalance is typically observed (input/output ratios of ∼0.1, Allard et al., 2006). The implication is that, during a paroxysmal sequence, the transfer of magma to the surface is facilitated and, as previously found (Coppola et al., 2019), the output/input ratio approaches 1 (e.g., all magma supplied to the shallow magmatic system is erupted).

Comparison between SO₂ flux and seismic tremor amplitude is also informative (Figs, 3, 4, 6). The seismic tremor record (Figure 3C) indicates clear high-frequency fluctuations (resulting from transient, intense seismic tremor increases during explosive activity) that are superimposed on a more continuous, longer-period trend. Focusing on this longer-term component (expressed by the 20 days moving average, exemplified by the red curve in Figure 3C), it is evident that seismic tremor is higher during a paroxysmal sequence (both PS1 and PS2) than during ordinary activity (see background seismicity in phases PPa, and PPP1 and PPP2). However, while the two paroxysmal sequences exhibit similarly elevated SO₂ degassing regimes (see Figure 3B; Figure 4B), there is a marked contrast in their seismic tremor signatures, with PS2 being associated to less intense tremor (compare red trend in PS1 and PS2; Figure 3C). Also, while the SO₂ flux escalation manifestly starts in phase PPb’, well before onset of paroxysmal activity (cfr. 4.2), the tremor increase in phases PPb and PPb’ is modest, with a notable escalation being only observed during PS1 (background tremor peaks on January 27). This complex SO₂ flux vs seismic tremor relationship is additionally indicated by Figure 6 that compares the time-averaged SO₂ flux vs seismic tremor signatures for the 7 activity phases. This plot confirms that the precursory phase (PPb’) is essentially manifested by temporally increasing SO₂ fluxes at nearly invariant tremor (relative to PPb), and that the PS1 and PS2 sequences have similar (within uncertainties) time-averaged SO₂ fluxes but diverse (> factor 2) seismic tremor amplitudes. Ordinary activity (phases PPa, PPP1 and PPP2) is characterised by consistently low SO₂ fluxes (<4000 Mg/day) and seismic tremor amplitude (<2 mV). The conclusion is that, although SO₂ flux and tremor are broadly correlated in the long-term (Figure 6; see also Salerno et al., 2018), their relationship is complicated by temporal swings in between the multiple tremor sources recognised on Etna (Cannata et al., 2008; Cannata et al., 2015; Zuccarello et al., 2022).

Although precise seismic tremor source location is beyond the scope of this study, we use the ratio of the seismic tremor amplitudes measured at the ECPN and ESLN seismic stations (e.g., the ECPN/ESLN tremor amplitude ratio; Figure 7) to investigate the temporal stability vs variability of the seismic tremor source in 2020–2021. This approach stands on the basis that temporal shifts in the ratio of tremor amplitudes recorded at summit (ECPN) vs intermediate height (ESLN) seismic station (see location in inset of Figure 7B and in Figure 1D) are informative of spatio-temporal changes in seismic tremor source (Alparone et al., 2003; Alparone et al., 2007; Cannata et al., 2008). We identify large oscillations in the ECPN/ESLN seismic tremor ratio in 2020–2021 (Figure 7) that are clearly indicative of a non-stationary seismic tremor source. These ratio fluctuations, in combination with information on seismic tremor location (extracted from INGV-OE reports; see Figure 7), support distinct PS1 vs PS2 tremor signatures; in which, during phase PS2, seismic tremor source is irradiated deeper (in May-June 2021; as indicted by low ECPN/ESLN ratio) and/or underneath BN/VOR crater sector (during July-October 2021), while seismic tremor source is more stably underneath the SEC in phase PS1 (Figure 7). A manifest upward migration (increasing ECPN/ESLN ratio) of the seismic tremor source, to cluster underneath the SEC, is also detectable in phases PPb and PPb’ (Figure 7).

Drawing on all the evidence, we propose a mechanism in which, during ordinary activity (e.g., phase PPa), volcanic activity is concentrated at the Central Crater (VOR) and NEC (Figure 2A; Figure 3A). The VOR/NEC vent area has long been recognised as the most stable (persistent) magma degassing pathway during ordinary activity, while the SEC is thought to represent the surface expression of a secondary lateral branching (connected to the VOR/NEC in a branched feeding conduit geometry) (Marchetti et al., 2009; Aiuppa et al., 2010; Cannata et al., 2015). Our derived low SO₂ fluxes (<2000 Mg/day) and magma input rates (∼1–3 m³/s) during the PPa phase (Figure 3; Figure 4; Figure 6) imply limited magma transport toward the south-eastern branching system (that feeds the SEC), and rather suggest preferential magma channelling in the “stable” central feeding conduit system. Seismic tremor locations clustering in the VOR/NEC area confirm this (INGV-OE, 2020; Figure 7). The radiant heat flux is low (Figure 3D), implying limited surface magma emplacement. This is typical of ordinary activity phases (e.g., Coppola et al., 2019), when a high input/output magma unbalance is typically observed that implicates shallow magma circulating in the central conduits evolves in open-system degassing conditions (Aiuppa et al., 2010; Edmonds et al., 2022). Open-system degassing implies efficient (low-pressure) separation of gas bubbles from melt to sustain quiescent open-vent degassing (and minor intra-crateric explosive activity) at VOR/NEC, and leads to recycling of dense/degassed magma deep into the plumbing system (Allard, 1997). Seismic tremor amplitude remains at low (background) levels (Figure 3; Figure 6).

Increasing SO₂ fluxes (Figure 3; Figure 4; Figure 6) and magma input rates (Figure 4C) in phase PPb imply preferential magma migration toward the southern part of the crater terrace, and its channelling into the structurally/gravitationally unstable (Bonforte et al., 2017) Etna’s south-eastern upper flank, where the SEC branching system insists (Marchetti et al., 2009; Aiuppa et al., 2010; Patanè et al., 2013). This increased magma supply rate causes progressive re-awakening of the NSEC (Figure 3A), and some moderate seismic tremor increase (Figure 3C; Figure 6). Notably, the seismic tremor source consistently migrates upward and toward the SEC in this period (Figure 7; INGV-OE, 2020). The process intensifies in phase PPb’ when (at relatively invariant tremor amplitude; Figure 6) the SO₂ flux ramps up to > 4000 Mg/day, suggesting accelerating magma transport (and degassing) toward the SEC. We propose that, as the magma input rate into the shallow SEC magma plumbing system increases to > 4 m³/s (PPb’ to PS1), such larger volumes of feeding mafic magma cannot be accommodated by quiescent (passive) open-system degassing, but rather lead to magma eruption in a paroxysmal sequence. We argue, based on previous work, that ascent of faster moving bubbly magma may lead to foam accumulation on-top of the SEC shallow magma reservoir (Allard et al., 2005), ultimately leading to pressure-build up (Bonaccorso et al., 2013; Bonaccorso et al., 2021), and finally eruption. Elevated seismic tremor and radiant heat at the PS1 onset mark explosive magma outburst in repeated paroxysmal episodes (Figures 3C,D). In this view, the sequence of paroxysmal episodes would perdure until fresh magma feeding into reservoir that feeds the NSEC vanishes down to the quiescent (persistent activity) values (in phase PPP1; Figure 4). Ultimately the process repeats again in May to December 2021, with phases PS1 and PPP2.

Our SO₂ flux record (Figure 3B) is evidence for the utility of spatially resolved UV Camera observations in detecting gas/magma migration prior to/during clustering of volcanic activity at one specific crater. In particular, the manifest SO₂ flux increase in phase PPb’ (Figures 3, 6), well before onset of paroxysmal sequence PS1, suggests that UV Camera observations have the potential to trace increasing magma supply weeks to months in advance to resumption of paroxysmal activity. However, comparison between the two recorded paroxysmal sequences (Figure 3; Figure 4; Figure 6) clearly shows that the transition from ordinary to paroxysmal activity can not only happen progressively (as in the case of PS1), but also suddenly (PS2). SO₂, tremor and thermal records all increase abruptly at the PS2 onset, and we see no detectable precursory SO₂ flux increase beforehand (Figure 3B). The diversity between PS1 and PS2 sequences is additionally demonstrated by the excess cumulative SO₂ mass, shown in Figure 4B (red curve labelled “Ex”). This is calculated by subtracting the average SO₂ flux for phase PPa (taken as representative of background degassing) from the cumulative SO₂ mass. This excess SO₂ mass progressively ramps up (during phases PPb to PPb’) prior to PS1, implying a gradual, protracted (several months long) acceleration in the SO₂ degassing regime. In contrast, the excess SO₂ mass does not show any precursory increase prior to onset of PS2, and rises contemporaneously with the start of paroxysmal activity. Tremor amplitude also rises during PS2, but to a lesser extent than in PS1 (Figure 3; Figure 6), implying distinct tremor source area (Figure 7). In particular, the prevailingly low ECPN/ESLN tremor amplitude ratios (Figure 7) suggest deepening of seismic source location in phase PPP1, and at the onset of phase PS2 (INGV-OE, 2020). This comparison demonstrates the challenges in anticipating paroxysmal behaviour at Etna, especially when the driving processes are rapid (hours) and deeply initiated (e.g., PS2).