2.3.4. Image Analysis and SO2 Flux Calculation

To enhance the contrast between volcanic emissions and atmospheric noise, and to limit dispersion effects and chemical conversions of SO2 in the atmosphere, image processing was conducted on a restricted image portion capturing an image sub-region above the crater area (Figure 2).

We also aimed at resolving gas contributions from different vents, and therefore, capturing changes in degassing dynamics and location [28]. To do this, we selected a rectangular sub-area over the crater terrace, along which we calculated the distribution of SO2 column densities and a plume velocity field (Figure 2). We calculated SO2 column densities and plume speed as close as possible to the vent, which minimizes the effects of wind and air entrainment within the plume that would produce dilution of SO2 concentrations farther downwind. This allowed us to detect changes in degassing dynamics across the crater terrace that were associated with changes in volcanic activity and regimes.

We then calculated velocity (mean, maximum, and associated standard deviation) and absorbance distribution along an ideal profile positioned in the middle of the sub-area, derived from averaging a series of parallel profiles within the area of analysis. From this, the SO2 density flux (in kgm−1s−1) was calculated by multiplying column densities associated with each pixel of the profile with the corresponding normal velocity component of motion (Figure 5, see also Reference [32]). Velocity profiles (Figure 5) were obtained by averaging the calculated two-dimensional velocity fields, and filtering out velocity points with low coherence. We also derived uncertainty in velocity determination along the profile by calculating the standard deviation associated with the velocity values used to determine the average velocity.

**Figure 5.** SO2 density flux calculated in the sub-area encompassing the summit craters. White arrows correspond to gas velocity vectors calculated on high coherence regions of the images. SO2 density flux distributions along the entire crater area and over the black dashed profile are calculated within the highlighted area, corresponding to vertical (Fy), western (FxW), and eastern (FxE) horizontal components, respectively. The SO2 total flux, for a given sector, was then calculated by integrating the density flux over the total length of the profile.

The SO2 flux was obtained by spatial 1D integration over the profiles. Discrimination of degassing contributions from the different craters, specifically the central (VOR+BN) and southeastern (SEC+NSEC) craters, was obtained by vectorial summation of the various SO2 flux components that border each crater. In doing so, we took into account if gas is moving away or toward the vent, which corresponds to a positive/negative flux contribution respectively.

Given the position of the station relative to the summit crater area (Figure 1), the gas contribution from the North-East crater (NEC) was not resolvable from our images, as hidden behind the SEC crater's ridge.

#### 2.3.5. Validation of the Automatic Method

Validity of the automatically processed SO2 fluxes was tested for a comparison with SO2 fluxes manually obtained using the Vulcamera software [33] (Figure 6). In particular, this was conducted for some selected days of acquisition characterized by a good weather condition and by clear and well-visible volcanic plume. The manual Vulcamera procedure involves calibrating images by individuation of a clear sky portion unaffected by degassing. Then, SO2 fluxes were calculated over an integration profile roughly perpendicular to wind direction, using plume speeds obtained from cross-correlation between consecutive frames along the selected profile [33].

**Figure 6.** Comparison between daily averaged SO2 fluxes calculated by a manual operator (using Vulcamera software, [33]) and those obtained by the automatic algorithm. The correlation coefficient between the two datasets is 0.75, with a 1:1 proportionality ratio. This correlation validated the automatic algorithm as an alternative, reliable method for SO2 flux calculation.

Comparison between manually and automatically calculated fluxes (Figure 6) demonstrated a correlation coefficient (R2) of ~0.75, with a best-fit regression line showing a ~1 proportionality factor. The main source of errors was associated with the algorithm for the automatic selection of the background sky area within an image, in particular during highly variable cloud conditions and strong winds. In such cases, background sky detection within images might not be optimal and may then result in a main underestimation of the real SO2 column densities, as shown by the position of outliers (Figure 6). Overall, this comparison validated the use of the automatic SO2 flux determination procedure, which paves the way to its full exploitation in real-time volcano monitoring, as already started on the Stromboli volcano (Italy) [32].

#### **3. Results: Application of the Automatic Real Time Algorithm: the Etna 2016 Case**

Etna is one of the volcanoes worldwide with the longest and most continuous SO2 flux record. SO2 flux measurements have become fundamental in volcano monitoring to define the rates of magma ascent and degassing within the shallow (<3 km) plumbing system. SO2 fluxes have been measured on Etna since the 1970s using the COSPEC (Correlation Spectrometer) [49,61–64] and, more recently, a network of Differential Optical absorption spectrometers (FLAME [65,66]). In view of its recurrent activity and robust past SO2 flux record, Mt. Etna is an ideal test site for validating our automatic processing method.

We reported below on the SO2 data automatically acquired and processed during 2016, which is a period characterized by substantial temporal changes in activity styles, including a phase of reduced degassing in the aftermaths of the December 2015 eruption [30,67]. This was followed by gradual activity escalation culminating into the May 2016 eruptive phase, and in a new vent opening episode on the eastern VOR crater rim in August [68–70].

Our aim was to test if different SO2 degassing regimes, related to such diverse activity styles, could be resolved and characterized in automatic (and in nearly real-time) using our permanent SO2 camera system.

#### *3.1. Etna's Activity in 2016, SO2 Flux Records, and Comparison with Seismic, Thermal Dataset, and Field Observations*

The SO2 flux time-series of 2016, which were generated by using the automated processing algorithm proposed in this study, is illustrated in Figure 7. The figure highlights that the significant variability in volcanic activity style in 2016 reflected a highly dynamic SO2 flux behavior (Figure 7). As illustrated in Figure 7, our 2016 temporal record shows daily averaged SO2 fluxes ranging between a few hundred to ~6000 tons per day (t/d). The associated standard deviations range from 100 to 4000 t/d. To assist interpretation of SO2 flux variations, we also reported seismic tremor and thermal radiance time-series (Figure 7), where the latter is expressed as Volcanic Radiative Power (VRP) [71].

**Figure 7.** (**a**) Cumulative daily averaged SO2 fluxes are compared with thermal cumulative activity detected by MIROVA. (**b**) SO2 flux daily average time-series are compared with seismic tremor amplitude (**c**), and thermal radiance (**d**). Three volcanic degassing levels are distinguished by different colors: green, low average degassing, orange, enhanced degassing, pink: intense degassing. The May 2016 Eruptive phase is preceded by more than one month long period of enhanced degassing activity when no significant increase in tremor and thermal is identified. Enhanced degassing activity is also detected between July and September, according to relatively higher tremor and thermal signals.

Combined analysis of field observations reported in Reference [56] including thermal, seismic, and SO2 fluxes time-series allowed us to distinguish three main periods of activity: (1) a pre-eruptive period (January to 16 May 2016), in which volcanic activity remained low (January to March) or gradually resuming (April to 16 May), (2) an eruptive period between 16 May and 25 May, characterized by intense strombolian activity, lava flows, and three short-lived lava fountaining episodes that

occurred in a brief time lapse between 18 May and 21 May, (3) a post-eruptive period (26 May to 31 December), that included a brief period of reduced activity until the end of June, which is followed by gradual (re)intensification of volcanic activity culminating with opening of new, strongly degassing incandescent vent at VOR on 7 August. This strong degassing declined during the subsidence of the Bocca Nuova (BN) crater's floor that occurred on October 10. These periods are described in more detail in the following sections.

#### 3.1.1. Volcanic Activity from 1 January 2016 to 15 May 2016 (Pre-Eruptive Period)

Volcanic activity from January to March 2016 was mainly characterized by sporadic ash emissions from the NSEC, and by passive degassing mainly from NEC. The daily average of SO2 fluxes fluctuated at low levels around ~1500 tons/day. The seismic tremor was stable at around low levels for Etna, and thermal activity occurred at reduced levels (<4 MW) (Figure 7).

Starting from the beginning of April, SO2 emissions gradually intensified and reached daily averaged fluxes of ~2000 t/d, and seismic tremor fluctuated within a subtle increasing trend. No significant thermal anomaly was still observed (Figure 7), and no significant volcanic activity change was reported from field observations, with NEC still passively degassing and NSEC producing a little more frequent ash emissions with occasional blocks ejected.

#### 3.1.2. The May Eruptive Period (16–25 May)

In the early morning of 16 May, strombolian activity resumed at NSEC and NEC, and became very strong at the latter crater on 17 May (also reported by References [68–70]). On 18 May at 10:50 UTC, a lava fountain started at VOR, which had been quiescent since 3 December, 2015. This event marked the beginning of a paroxysmal sequence, lasting until 25 May. The sequence included two additional short-lived lava fountaining episodes at VOR, on 19 May and 21 May, and ended with an intense strombolian and lava flow activity that lasted several hours (Figures 7 and 8). Lava effusion accompanied all the strongest explosive episodes, issuing from both fissures on the summit cone and overflow from its crater rim. In particular, overflowing from the BN crater rim formed a large lava field that extended downslope up to 3 km. The cumulative volume of lava flows and pyroclastic fall deposits was preliminary estimated at 7 to 10 Mm<sup>3</sup> [56], which is similar to what erupted during the December 2015 paroxysmal sequence [30]. Lastly, the summit crater's area was affected by intense deformation with fracturing, subsidence, and a formation of a ~1 km-long and nearly NS oriented graben-like structure (Figure 1b).

**Figure 8.** SO2 fluxes measured from 11 May 2016 to 31 May 2016 encompassed the 16–25 May eruptive period (blue bars show the maximum and minimum values while white dot indicates the daily average), which show a good agreement with seismic tremor amplitude, with relative maxima fitting the phases of enhanced seismic activity linked with paroxysmal events during the eruptive period.

Thermal radiance and the seismic tremor showed coherent increases of three and two order of magnitude, respectively (Figures 7 and 8). The three lava fountains were clearly marked by short-lived peaks in seismic tremor amplitude on 18, 19, and 21 May (Figure 8). The 24–25 May episode of intense strombolian activity/lava fountaining was associated with a wider, longer-lived phase of seismic tremor increase (Figure 8).

The daily averaged SO2 fluxes increased up to 6000 tons/day (Figure 7), pointing to heightened degassing, nearly tripled with respect to the pre-eruptive phase. A detail of the SO2 fluxes recorded during the May 2016 paroxysmal sequence is illustrated in Figure 8, where alternation of eruptive and repose periods were evident in the degassing record, with peaks in daily averaged SO2 flux in three (18, 21 and 25 May) out of four days of paroxysmal activity.

It is also worth noting that only the first (18 May) lava fountaining episode occurred during the UV-camera acquisition temporal interval (see Section 3.2). The high daily averaged SO2 fluxes obtained on this specific day (Figure 8), thus, reflected the "explosive" SO2 contribution during the paroxysmal event. In contrast, the SO2 flux peaks on 22 and 25 May could not be explained by syn-explosive SO2 release (the lava fountains occurred outside the camera acquisition hours), but rather reflected heightened passive degassing, and/or milder (strombolian) explosive activity, prior to/after the paroxysmal episode itself.

#### 3.1.3. Volcanic Activity after the Eruptive Phase (26 May–31 Dec)

Reduced degassing emissions (<2000 t/d) were measured after the eruptive phase, from May 26th until the end of June (~1 month), which they could interpret as the aftermath of voluminous gas/magma release during the previous December 2015 [30,39,67,72] and May 2016 eruptive sequences. Field observations indicated that the summit craters were weakly fuming and occluded by lavas and pyroclasts. A progressive subsidence occurred on the VOR crater's floor, where lunar cracks formed on the crust of the spatter deposits that had filled the crater. Since early July, intensification of SO2 degassing (average daily fluxes increased from ~900 to ~3000 tons/day, Figure 7) was accompanied by crack widening near the eastern VOR's crater rim, along the graben-like structure (Figure 1), which culminated, on 7 August, by opening a new 20-m large pit vent, characterized by vigorous high temperature degassing and glowing at night (Figure 1c) [69]. Consistently, thermal activity, as detected by MIROVA, increased in early July from <5 to ~10 MW, which was also reported by Reference [69], and a notable increase was also observed on the seismic tremor (Figure 7). However, repeated inspections on the summit area did not reveal any evidence of explosive activity such as fallout material deposited around the pit vent (as reported in Reference [56]).

From 10 October, after a small explosion, a large subsidence of the BN north-western inner floor was observed. This affected lavas and pyroclastic materials that had filled the central craters during the 2015–2016 paroxysmal sequences. This inner crater subsidence, characterized by episodic collapses, was nicely paralleled by declining SO2 fluxes down to low values (~1000 tons/day on average). RMS seismic amplitude and thermal radiance slightly decreased as well.
