Monitoring endogenous growth of open-vent volcanoes by balancing thermal and SO2 emissions data derived from space

Measuring the amount of magma intruding in a volcano represents one of the main challenges of modern volcanology. While in closed-vent volcanoes this parameter is generally assessed by the inversion of deformation data, in open-vent volcanoes its measurement is more complicated and results from the balance between the magma entering and leaving the storage system. In this work we used thermal and SO2 flux data, derived from satellite measurements, to calculate the magma input and output rates of Mt. Etna between 2004 and 2010. We found that during the analysed period more magma was supplied than erupted, resulting into an endogenous growth rate equal to 22.9 ± 13.7 × 106 m3 y−1. Notably, this unbalance was not constant in time, but showed phases of major magma accumulation and drainage acting within a compressible magma chamber. The excellent correlation with the inflation/deflation cycles measured by ground-based GPS network suggests the thermal, SO2 flux and deformation data, can be combined to provide a quantitative analysis of magma transport inside the shallow plumbing system of Mt Etna. Given the global coverage of satellite data and the continuous improvement of sensors in orbit, we anticipate that this approach will have sufficient detail to monitor, in real time, the endogenous growth associated to other world-wide open-vent volcanoes.

where ρ m is the dense rock equivalent (DRE) magma density (kg m −3 ) and 86.4 is the unit's conversion factor (from tonnes day −1 to kg s −1 ). The magma output rate (Q out ) is instead calculated from repeated measurements of volumes erupted over discrete time intervals, or through a space-based thermal approach 5 . The latter, also called "thermal proxy", derives from an energy balance of lava flows, whereby there exists a proportionality between the area, the surface temperature, and the discharge rate that feeds lava flows [6][7][8] . Hence, the DRE magma output rate (Q out ), can be calculated from the satellite-derived Volcanic Radiative Power (VRP) as 8 : (2) out lava m where ρ lava (kg m −3 ) is the bulk densities of the lava flow and c rad (J m −3 ) is an empirical best-fit parameter that relates the lava discharge rate to the thermal radiation for any given rheological case 8 . To these two rates of magma transport, a third parameter of vital importance needs to be considered to fully constrain the magma budget of a volcano. This is represented by the rate of magma accumulation, or withdrawal (dQ m = Q in − Q out ), that determines if a volcano is inflating or deflating 9,10 . However, because magma is compressible, only a fraction of the net magma volume entering or leaving the chamber (ΔV m = V in − V out ) will be accommodated by the expansion/contraction of the chamber itself (ΔV c ), so that 10,11 : where β m and β c are the compressibilities of magma and chamber, respectively. In turn, subsurface volume changes (ΔV c ) are generally constrained from the surface displacements observed by geodetic techniques, such as GPS and InSAR 12 , although source geometry and material properties can also play a fundamental role in modulating the ground deformation 11,13 . It is therefore clear that there must be a direct correlation between the rates of degassing, eruption and deformation, through which the magma budget is satisfied (eq. 3).
At closed-vent volcanoes (i.e. volcanoes with sealed conduits that may exhibit edifice-wide long-term ground uplift due to reservoir pressurization 14 ), the classic model of the "volcano deformation cycle", is generally operative 2 . According to this model, during inter-eruptive periods, the magma output rate (Q out ) is equal to zero, so that all the magma supply (Q in ) gradually inflates a chamber beneath the volcanic edifice, causing uplift. Once the eruption starts (i.e. once a critical pressure is reached at which point the chamber ruptures), Q out becomes much higher than Q in , so the volcano deflates, because of withdrawal from the magma chamber 9,10 . In some cases, this long term inflation/deflation cycle has shown predictable behaviours, allowing successful forecast of eruptions [15][16][17] .
On the other hand, it is now well established that open-vent volcanoes (i.e. volcanoes characterized by persistent degassing and continual emission of magmatic-related products directly to the atmosphere 18 ), emit an excess of gas to that contained in the erupted magma 19,20 . This unbalance, commonly known as "excess degassing" 19,20 , has been ascribed to phenomena of endogenous or cryptic growth, in which some unerupted magma is continuously intruded in the volcano edifice, or ascends, degasses and then sinks back down the conduit to mix with a deep reservoir 1,4,[21][22][23][24][25] . In most open-vent volcanoes this large unbalance does not translate into significant long-term surface deformation, but only into short-term, localized deformation related to shallow conduit processes 14 . For example, short-term deformation cycles have been observed in association with Strombolian or Vulcanian eruptions 26,27 as well as with variations in the lava lake level 28 . Also, is not uncommon to observe deflation phenomena during major eruption phases, that temporary drain the magma chambers 9,29 or the upper portion of the plumbing system [30][31][32] . According to some authors the lack of long-term inflation at these open-vent volcanoes cannot be explained by magma reservoirs being too deep to create detectable uplift, but rather by the possibility of magma to rise toward the surface without pressurizing the magma reservoirs 14 .
Therefore, is still unclear how the large amount of unerupted magma ΔV m, characterizing open-vent volcanoes is related to the subsurface volume changes ΔV c inferred from ground deformation, although the two volumes are clearly linked (eq. 3).
Mt. Etna has been object of several works focused on the balance between supplied and erupted magma, all showing that only a small fraction of magma entering the degassing zone (3-4 km b.s.l.), is actually extruded 23,29,33 . Most studies have examined this open-vent behaviour over periods of decades, suggesting an overall endogenous growth by intrusion and the creation of a cryptic plutonic complex within Etna's sedimentary basement 23,29,33,34 . Some detailed magma budget was also presented for three individual eruptions of Mt. Etna 35 , suggesting that over timescale of weeks, the balance between input and output rates is variable and much more complex than previously thought. However, this study did not consider the magma budget during the periods between eruptions, thus lacking one half of the whole history. First clues on Etna's magma budget during a sequence of short-lived paroxysms 36 suggests that the majority of the magma erupted during these events was likely stored within the conduit, having gone through extensive degassing for days to weeks prior to the paroxysms. Notably, all these works 23,29,[33][34][35][36] focused on the magma budget of Mt. Etna over different timescale (from decades to days), but none looked at the relationships with the volcano deformation over the same timescales.
On the other hand, ground and satellite measurements of deformation suggest that Mt. Etna is affected by an overall long-term inflation trend, periodically interrupted by more or less important deflation periods 29,37-40 , a pattern somehow similar to the classic closed-vent volcano deformation cycle. One work first correlated this deformation pattern with the long-term bulk magma budget characterizing the 1993-2005 period 29 , but it does not have the sufficient temporal detail to explore the validity of the above model in the mid-, or short-term.
In this work we used satellite thermal data, derived from Moderate Resolution Imaging Spectroradiometer -MODIS (Fig. 1a) and SO 2 flux data, derived from Ozone Monitoring Instrument -OMI (Fig. 1b), to constrain the magma supply and output rates (methods) characterizing the activity of Mt. Etna between 2004 and 2010 (Fig. 1c). These two independent datasets provide for the first time a clear, quantitative view of the magma accumulation www.nature.com/scientificreports www.nature.com/scientificreports/ and withdrawal processes affecting the volcano plumbing system at different timescales. These results are finally compared with the deformation data measured by the CGPS network during the same period 38,39 and highlight the great potential of using combined satellite data to monitor quantitatively the endogenous growth of open-vent volcanoes.

Activity of Mt. Etna between 2004 and 2010
After ∼20 months of quiescence, and without precursory seismicity or ground deformation, a new eruption According to eq. 2, the thermal flux recorded on 17 November 2006 was produced by a magma output rate, Q out , equal to 6.5 ± 1.9 m 3 s −1 . The evident unbalance between Q in and Q out suggests the empting of a shallow magma chamber consistent with a stage of deflation of the volcano edifice, as measured by GPS network (cfr. Fig. 4 www.nature.com/scientificreports www.nature.com/scientificreports/ petrologic data suggest that the magma feeding the 2004-2005 activity was likely to have been stored in the shallow plumbing system during the 2000 and 2001 activity, where it experienced volatile loss and extensive crystallization 41,43,44 . Few months after the 2004 eruption, between November 2005 and January 2006, the amplitude of volcanic tremor increased abruptly, although neither effusive nor paroxysmal activity was observed at the summit vents 45 . This event was associated with mild inflation of the summit of the volcano and with increased SO 2 flux, allowing some authors to interpret it as a failed eruption, resulting from recharging of the volcanic feeder at depth (>3 km b.s.l.) 45 . Following this intrusive event, minor explosive activity occurred in the summit area (Bocca Nuova -BNand Voragine -VOR -craters; Fig. 1c) during 2006, whereas the Northeast Crater (NEC; Fig. 1c) was the site of deep-seated explosions, which increased in early July 2006. The heightening of NEC activity was the prelude to the 2006 eruption that started on 14 July and that was characterised by two main phases of activity 46 . The first phase (14-24 July) lasted 10 days during which Strombolian and effusive activity (emitting 3-6 × 10 6 m 3 of lavas) took place along a short fissure on the lower ESE flank of the SEC cone 47 . The second phase lasted from 31 August until 14 December and consisted of 20 eruptive episodes at or near the summit of the SEC cone producing ash plumes. From 12 October onward, an almost persistent effusive activity took place from a number of vents opened on the south of SEC 46 . (Fig. 1c). As a whole the 2006 eruptive activity produced two main lava fields emplaced on the south-west and on the south-east flank of the volcano, with a total lava volume of 20-37 × 10 6 m 3 42,46 . www.nature.com/scientificreports www.nature.com/scientificreports/ After the end of the 2006 eruption, a recharge phase was indicated by continuous inflation showed by GPS measurements 38,39 . The recharge phase was accompanied by seven short-lived lava-fountaining episodes that erupted a total of ∼12 × 10 6 m 3 of lava, and marked the shifting of activity from the SEC to the newly formed new SE crater (NSEC) [48][49][50] . The last powerful episode occurred on 10 May 2008, only 3 days before the onset of the 2008 eruption 51 . Unlike the previous two eruptions, the initial phase of the May 2008 eruption was preceded and accompanied by strong seismic release and by marked deformation of the volcano flanks 52 .
Moreover, in this case, the eruptive fissure propagated inside the Valle del Bove, producing the long-lived 2008-2009 eruption that lasted 419 days (the longest eruption since 1991-1993) and that emitted 68-74 × 10 6 m 3 of lava 42,53 (Fig. 1c). After the end of this eruption Mt. Etna was characterized only by continuous degassing and by sporadic minor explosive phenomena for 17 months 49

Results
The time series of Volcanic Radiative Power (VRP) recorded between 2004 and 2010 (methods) outlines a continuous thermal emissions (∼1.65 thermal detections per day) with VRP ranging from less than 1 MW to ∼8.5 GW (Fig. 2a).  Fig. 2a) are clearly visible and they separate from the low and continuous thermal emission that characterise the entire period investigated. The statistical analysis of the VRP (log-transformed), indicates a bimodal distribution, with two distinct regimes within the entire population of data (n obs = 3320; Fig. 2c). The low thermal regime, characterized by VRP < 40 MW ( < 7.6 log-transformed), was typical of inter-eruptive periods, when the thermal source was likely produced by weak strombolian activity, hot degassing cracks within the summit craters, as well as by the cooling of previously emplaced lava flows (during the post-eruptive phases; Fig. 2a). Conversely, the high thermal regime (VRP > 40 MW) corresponds to periods of major effusive activity, which also includes the 7 lava-fountaining episodes occurred in 2007 and 2008 (Fig. 2a). It is thus clear that a change from low to high thermal regimes (peaking at 4 MW and 252 MW, respectively) corresponds to a transition from open-vent to effusive activity, in manner similar to that already observed at Stromboli volcano 54 .
The SO 2 flux data ( Fig. 2b) also indicate a continuous degassing between 2004 and 2010 (∼0.37 SO 2 detections per day), but with an unimodal distribution of the data (n obs = 753; Fig. 2d). The average φSO 2 obtained from the www.nature.com/scientificreports www.nature.com/scientificreports/ OMI data (Methods) is equal to 3259 tonnes day −1 , in excellent agreement with the flux measured by ground measurements during 2005-2008 55 (see methods). It is worth nothing that during the inter-eruptive phases, the average flux was 3023 tonnes day −1 , while during the 3 effusive phases (grey fields in Fig. 2) the average φSO 2 raised to 4301 tonnes day −1 . This suggests that the effusive eruptions of Mt. Etna were characterized by a bulk increase of the degassing rate of only ∼42% with respect to the persistent, long term emission.
Notably, the failed eruption that occurred during the winter of 2005-2006 45 (yellow field in Fig. 2) was accompanied by a monthly-long pulse of φSO 2 , while there were no anomalous thermal emissions recorded in the VRP time series (the cooling trend remained essentially unperturbed; Fig. 2a). In contrast, all the lava-fountain episodes are evident in the VRP time series (Fig. 2a), but not particularly in the This may be due to the impulsive nature of the phenomenon (i.e. duration less than 1 day, on average) which coupled with the low sampling rate of OMI (1 image per day) makes challenging the detection of all short-lived SO 2 emissions associated to these episodes 36 .
In general, it can be observed that both thermal and gas emissions have been continuous during 2004-2010 (Fig. 2a,b), remarking the persistent open-vent behaviour of Mt. Etna. Nevertheless, the thermal flux shows wide and evident increases during major eruptions (almost two orders of magnitude), leading to a bimodal distribution of VRP (Fig. 2c). On the contrary, φSO 2 shows an unimodal distribution (Fig. 2d) and a more homogeneous trend, with smoothed, monthly-long oscillations accompanying the major eruptions and the intrusive episode (Fig. 2b).

Discussion
The magma supply rate and the magma output rate have been estimated using eqs 1 and 2, whose detailed parametrization is described in the method section. However, some aspects must be clarified regarding the interpretation of these parameters and their volcanological significance at open-vent volcanoes.
The magma supply rate calculated by eq. 1 (Q in ) is in fact a "magma degassing rate", since it provides an estimate of the quantity of magma that rises above the exsolution level of SO 2 and degasses 1,19,23 . Therefore, it should be interpreted in a broad sense, as the rate of magma entering the shallow portion of the magma system where a convective degassing cell is operative (about 3-4 km b.s.l. for Etna 4,23,33 ).  38 . Note the good correlation between the overall trends of the two datasets, as well as the excellent correspondence between the main phases of magma accumulation/withdrawal (a), leading to coeval and coherent inflation/deflation cycles (b). Coloured fields and bars as in www.nature.com/scientificreports www.nature.com/scientificreports/ On the contrary, the magma output rate (Q out ) calculated through eq. 2 provides an estimate of the lava discharge rate during the effusive phases 8 . However, during phases of open-vent activity, the magma output rate assumes an ambiguous meaning, and should be considered an "apparent output rate", since there is no net extrusion of lava. In this case, the application of eq. 2 can still be used to infer the rate at which the magma reaches the surficial levels in the conduit (i.e., few tens/hundreds meters at most), where it radiates and cools before being cycled back 8,36,56 . In this scenario, the threshold separating the low thermal regime (open-vent activity) from the high radiating regime (effusive activity) indicates that an ascending volumetric flux of 0.16 ± 0.05 m 3 s −1 represents a physical barrier for the convective magma transport within the conduit(s) of Mt. Etna. Above this threshold the magma can no longer be recycled in the conduit and must be partially or completely extruded.
Accordingly, the two magma fluxes (Q in and Q out ) sample the rate of magma transport at different depths 56 , being the SO 2 -derived flux relative to the magma circulation above the exsolution level of SO 2 , and the IR-derived flux relative to very shallow (i.e. during open-vent) or surficial (during effusive activity) levels.
The time series of Q in and Q out , as well as the respective cumulative volumes (V in and V out ) estimated for the 2004-2010 period, are shown in Fig. 3. The direct comparison of these rates clearly shows the overall unbalance due to the excessive degassing of Mt. Etna, with supply rate generally exceeding the output rate. This unbalance is particularly evident during inter-eruptive periods, while it is cancelled, or sometimes reversed, during the major eruptive phases (Fig. 3a).
Although the long-term magma budget indicates a general accumulation of magma, the trend shows midand short-term variations of ΔV m that can be compared with the deformation data reported for the same period. In this regard, we have used the published baseline distance between two CGPS stations (located on the west side of Etna; inset in Fig. 4b), which highlight the deformation pattern characterizing Mt. Etna activity between 2004-2010 38 (Fig. 4b).
This comparison shows an excellent correlation between the two datasets, suggesting that the magma budget retrieved from SO 2 and thermal data is closely linked to the volcano deformation. All the main stages of inflation (positive ΔV m ) are actually corresponding to periods of "excess degassing" produced by magma accumulation at depth, without net extrusion of magma. These stages essentially correspond to the open-vent activity, or to periods characterized by magmatic intrusion, like the one occurred in the winter 2005-2006 (yellow field in Fig. 4). On the other hand, all periods of magma drainage (negative ΔV m ) occurred during deflation stages, and in particular during specific periods of the main effusive eruptions (grey fields in Fig. 4a).
By using similar GPS analysis, other authors 39 identified 6 stages of inflation/deflation occurred between 2004 and 2007 (named from T5 to T10; Table 1). For each stage they estimated the depth (D) of the deformation source as well as the magma chamber volume changes (ΔV c ) that better explain the ground deformation pattern (Fig. 5a). These data indicate the presence of a magma chamber at 3.5-5.5 km b.s.l. where the pressurization/ depressurization occurs. The volumes changes (ΔV c ) inferred for each stage span from −8.4 × 10 6 m 3 (T5) to +6.4 × 10 6 m 3 (T7), and compare very well with the volumes calculated from the magma budget approach (ΔV m ; Table 1). This excellent correlation (R 2 = 0.989; Fig. 5b) indicates a ΔV m /ΔV c ratio spanning from 2.4 to 6.5, with a best-fit value equal to 3.6 (Fig. 5b). As described previously this ratio describes the compressibility of the system (eq. 3), which for bubble poor basaltic magmas stored at shallow crustal levels spans from 1.2 to 7.7 10,11 . These results are thus consistent with the occurrence of efficient magma convection in the upper magmatic systems of Mt. Etna allowing degassed magma to sinks back and accumulate within the deforming magma chamber.
These results have a series of important implications for the long-and short-term monitoring of volcanoes showing persistent activity.
Firstly, we demonstrated that the endogenous growth of Mt. Etna volcano can be quantified in great detail by means of satellite-derived thermal and SO 2 flux data, and that the resulting magma budget (ΔV m ) is compatible with the volume changes (ΔV c ) inferred by deformation measurements, and with the bulk compressibility of the magmatic system.
Starting from this evidence the following considerations can be drawn for any open-vent volcano:  www.nature.com/scientificreports www.nature.com/scientificreports/ • each open-vent volcano that degasses more magma than it emits, grows in an endogenous way; • this growth can be monitored by means of satellite-derived thermal and SO 2 measurements (allowing the calculation of ΔV m ; eqs 1 and 2), or by means of deformation measurements (allowing the calculation of ΔV c ; eq. 3); • The ratio between ΔV m and ΔV c is a function of the compressibility of the system and can be used to to infer depth of the accumulation zone as well as the magma and the crustal parameters; • The continuous monitoring of ΔV m and ΔV c allows to quantify the shallow magma budget of the system. Each temporal decorrelation between these two volumes is indicative of a change in (i) position of the magma accumulation zone (ii) initial volatile's content (iii) thermo-elastic properties of the system. • if open-vent activity is not accompanied by long-term detectable deformation (i.e. Merapi, Popocatépetl, Colima 14 ) the endogenous growth possibly occurs within a very compliant crust or at shallow depths, or in correspondence of a bubble-rich, highly compressible magma chamber, characterized by a high ΔV m /ΔV c ratio 13 .
We also demonstrated that thermal and gas emissions are not necessarily correlated throughout effusive activity 3 , nor during open-vent activity (Fig. 2). This implies that, even at persistently degassing volcanoes, the high-temperature thermal anomalies measured by satellite 57 (i.e. Te > 600 K), are not necessarily correlated to the emission of gas itself, as hypothesized by 58 . In fact, the magmatic intrusions occurred between December 2005 and January 2006, although it produced an evident pulsation of SO 2 flux from the summit craters (>20 ktonnes day −1 ), it was not accompanied by evident thermal anomalies, just because the magma remained at depth (>3 km b.s.l. 45 ). Accordingly, we here stress that also at persistently degassing volcanoes, the detection of MIROVA thermal anomalies, advise that the magma column has reached the surface, or is at very shallow depth.
A first recent comparison of thermal anomalies, degassing, and surface deformation, derived from satellite data at 47 of the most active volcanoes in Latin America, provides clear evidence of the potential of applying this approach at regional scale 59 . Equations 1, 2 and 3 can be applied in near real time if appropriate parametrization has been validated such as for Mt. Etna. We believe that in the future, the technological advancement represented www.nature.com/scientificreports www.nature.com/scientificreports/ by the new infrared, ultraviolet and microwave satellite's sensors will allow a sufficient degree of detail for monitoring on a daily basis the endogenous growth of open-vent volcanoes and, more generally, the efficiency of magma transport within the magmatic plumbing systems.

Methods
Volcanic Radiative power (VRp) via MoDIs-MIRoVA data. Thermal emissions of Mt. Etna, have been obtained by analysing the data provided by the Moderate Resolution imaging Spectroradiometer (MODIS). Two MODIS sensors, in orbit since February 2000 and May 2002 (on board of Terra and Aqua NASA's satellites, respectively) provide, on average, 4 images per day of the entire globe, with a spatial resolution of 1 km in the infrared. The data were processed through the MIROVA system 57 , an automatic volcanic hotspot detection system, based on the analysis of Middle Infrared (MIR) images. This system allows to detect, locate and quantify the thermal emissions sourced by volcanic activity by calculating the Volcanic Radiative Power (VRP). The VRP is a measurement of the heat flux radiated almost exclusively by hot lava surfaces (i.e. having effective temperature above 600 K) and is calculated via the MIR method 60 as: where R MIR,alert is the pixel integrated MIR radiance of the i th alerted pixel, R MIR,bk is the MIR radiance of the background, A pixel is the pixel size (1 km 2 for the resampled MODIS pixels), and 18.9 is a constant of proportionality 60 . The VRP has an uncertainty of ±30% 60 which cover the possible insulation conditions (integrated surface temperature from 600 K to 1500 K) of the observed lava flow (not shown in Fig. 2a for graphical clarity).
In this work, we analysed 8903 images acquired over Mt. Etna between 1 August 2004 and 31 January 2010. Within this dataset MIROVA detected a thermal anomaly in 3320 images (∼37%) allowing the construction of the VRP time series shown in Fig. 2a. The presence of clouds (or ash plumes), as well as the poor viewing geometry conditions may prevent or affect the thermal signal detected by the satellite 57 . In volcanoes that exhibit persistent thermal anomalies (i.e more than one a day; as for Etna) a practical way to overcome this attenuation is to considered the maximum daily VRP (Supplementary Dataset) as the most reliable value to calculate the daily magma output rate 57 . The latter has been calculated by using the eq. (2), where a single coefficient, called radiant density (c rad in in J m −3 ) describes the relationship between volumetric and radiant flux appropriate for the observed eruption 8 . An exhaustive review of the application of this method over basaltic lava flows is provided by 61 .
Here we used a radiant density of 2.5 × 10 8 J m −3 which has been already calibrated for Etnean lava and provides estimates of magma output rate with an uncertainty of ±30% 8,57 . According to 57 , this uncertainty embeds all the possible emplacement conditions (rheological, topographic and insulation) that characterize the Etnean lava flows. In order to obtain dense rock equivalent magma output rates (Q out ), we assumed ρ m = 2600 kg m −3 , and ρ lava = 2350 kg m −3 to give a bulk flow vescicularity of about 10% 62 . so 2 flux (φso 2 ) via ozone Monitoring Instrument (oMI). The OMI is one of the four instruments on board of AURA platform dedicated to monitor solar backscatter radiation at wavelengths spanning from 270 to 500 nm (visible and ultraviolet). It provides a daily global coverage in 14 orbits since the 1 st October 2004. Each image has a complete swath of 2600 km and a nominal pixel spatial resolution of 13 × 24 km at nadir. Here we used the OMSO2 Level 2G (0.125° × 0.125°) global product 63 that contains estimates of the total SO 2 vertical column density (VCD SO2 ; in D.U.), assuming different centres of mass altitudes (CMAs), at 0.9 km (Planetary Boundary Layer, PBL), 3 km (Lower tropospheric, TRL), 8 km (Middle tropospheric, TRM), and 17 km (Lower stratospheric, STL), respectively. Since 2010, a technical issue on the instrument's field of view makes certain rows of the OMI swath unusable for SO 2 retrieval. This led the occurrence of 2-3 days-long gaps between consecutive OMI observations that corrupt the quantification of the long-term degassing budget. Here, we processed OMI data acquired between 1 st October 2004 to 1 st February 2010, consisting of 1831 images.
To detect SO 2 emitted from Mt. Etna we adopted the same procedure described in 64 , composed by the following steps: (i) Extraction and cropping of mask of 20° × 20° centred on Mt. Etna, for any CMA product; (ii) Calculation of a contextual threshold at PBL level, defined by T so2 = μ + 3σ, where μ and σ represent the mean and the standard deviation of those pixels within the mask having VCD SO2 < 1 D.U (at PBL level); (iii) Recognition and definition of groups of adjacent pixels (clusters) that positively exceed this threshold; (iv) Visual inspection and manual selection of all the clusters attributed to the Etna's plume.
For any image where at least one cluster was selected, we calculated the SO 2 mass burden (M SO2 , in tonnes), at PBL, TRL, TRM and STL layers, by using the following equation, as proposed by 65  where npix, A pix and VCD SO2 represent, the number of alerted pixels, the pixel area (in km 2 ) and the vertical column density of SO 2 (in D.U.), of the selected plume. For Mt. Etna plume we used the SO 2 mass retrieved for PBL and TRL layers, because the information on the typical altitude reached by the volcanic plume related to the activity of the analysed volcanoes (less than 5 km). In particular, we used the PBL level for the most of the cases which provide good results to quantify passive ( 66 . However, during more powerful plume emissions, hereby defined as plumes producing more than 20 ktonnes in the PBL, we used the mass burden calculated for TRL layer. To convert the mass burden (M SO2 ) into SO 2 flux (φSO 2 ) we used the "delta-M" method 65 , according which: where Δt is the time interval (days) between two consecutive mass estimates M i and M i−1 , and τ is the SO 2 e-folding time, an SO 2 loss term that takes into account chemical, transport and dilution processes inside the plume 67 . We settled τ equal to 2 days, which has been proved to be an optimized value for a continuously degassing volcano in a stable wind regime 68 . The resulting time series consists of 754 discrete measurements of φSO 2 (Supplementary Dataset) which are compared with ground measurements (Fig. 6) obtained with the FLAME network and traverse method for the same period 55 . The uncertainty in the estimation of φSO 2 using eq. 6 derives from the complex combination of the uncertainties in the selection of the proper center of mass altitudes (hence M SO2 ) as well as on the assumed SO 2 loss term (τ). Although the exact error of each of these parameters remains difficult to quantify (since it may vary from image to image) the excellent agreement with ground measurements, in terms of average values (3259 tonnes day −1 from OMI; 3530 tonnes day −1 from Flame network; 3250 tonnes day −1 from traverse method) and trend (Fig. 6), provides a first order validation for φSO 2 derived from OMI whose uncertainty is considered comparable with ground measurements (i.e. ±30% 55 ). Finally, to convert the φSO 2 into magma supply rate (Q in ) we used the eq. (2) by assuming a total sulfur loss ΔS, equal to 3600 ppm. This value is derived from the maximum initial S content measured on melt inclusion of Mt Etna 69 and assuming a complete volatile loss during magma ascent. Note that the assumption of a lower initial sulphur content or an incomplete loss of the volatile during the rise of the magma, produces higher values of Q in . Consequently, the results from eq. 1 are to be considered minimum values on which the uncertainty of φSO 2 (±30%) is propagated.

Data Availability
The satellite datasets (VRP and φSO 2 time series) and the cumulative volumes (Vin and Vout) are available as Supplementary Datasets.