Changes in SO 2 Flux Regime at Mt. Etna Captured by Automatically Processed Ultraviolet Camera Data

: We used a one-year long SO 2 ﬂux record, which was obtained using a novel algorithm for real-time automatic processing of ultraviolet (UV) camera data, to characterize changes in degassing dynamics at the Mt. Etna volcano in 2016. These SO 2 ﬂux records, when combined with independent thermal and seismic evidence, allowed for capturing switches in activity from paroxysmal explosive eruptions to quiescent degassing. We found SO 2 ﬂuxes 1.5–2 times higher than the 2016 average (1588 tons / day) during the Etna’s May 16–25 eruptive paroxysmal activity, and mild but detectable SO 2 ﬂux increases more than one month before its onset. The SO 2 ﬂux typically peaked during a lava fountain. Here, the average SO 2 degassing rate was ~158 kg / s, the peak emission was ~260 kg / s, and the total released SO 2 mass was ~1700 tons (in 3 h on 18 May, 2016). Comparison between our data and prior (2014–2015) results revealed systematic SO 2 emission patterns prior to, during, and after an Etna’s paroxysmal phases, which allows us to tentatively identify thresholds between pre-eruptive, syn-eruptive, and post-eruptive degassing regimes.


Introduction
Active volcanoes are monitored by a variety of increasingly sophisticated volcanic gas sensing techniques [1][2][3], which are contributing to improved understanding and monitoring of volcanic activity. In particular, trends in volcanic gas composition and fluxes help detect subtle changes in the rates of magma ascent and degassing within shallow volcano plumbing systems, and allow to better confinepre-eruptive and syn-eruptive processes (gas in magmas being the main drivers of volcanic processes) [4]. Nevertheless, gas monitoring still lags behind more established seismic and geodetic techniques [5]. This is because the low temporal resolution of gas observations has traditionally hampered real-time analysis of fast-occurring volcanic processes, such as shallow intrusion of magma prior to eruption, and rapid gas ascent and release during explosive eruptions.
Volcanic SO 2 emissions provide key information on the rates of magma ascent in shallow (<3 km) magma plumbing systems, and are extensively monitored worldwide using instrumental networks of scanning ultraviolet spectrometers using the Differential Optical Adsorption Spectroscopy (DOAS)

UV Camera System
In the framework of the ERC-funded project BRIDGE (www.bridge.unipa.it), we designed a multi-instrument UV-absorption spectroscopy system for robust SO 2 flux measurements ( Figure S1). The system is composed of (i) an instrument module and (ii) an acquisition/processing module. The instrument module is equipped with two JAI CM-140GE-UV cameras sensible to UV-radiation, and one Ocean-Optics USB2000+ Spectrometer coupled to a telescope of rectangular, vertically-oriented Field Of View (FOV ≈ 0.3 • × 14 • ), and is spatially filtered to match the ≈12 • vertical width. Two different band-pass optical filters with Full Width at Half Maximum (FWHM) of 10 nm, and a central wavelength of 310 and 330 nm, respectively, are applied in front of the cameras to enhance differential Remote Sens. 2019, 11, 1201 3 of 21 UV absorption in the SO 2 bandwidth [45,46]. In addition, 520 × 676 pixel images are acquired at 10-bit resolution with a frame rate of 0.5 Hz.
To obtain a quantitative measure of SO 2 column density within the volcanic cloud, we calculate the proportionality ratio between absorbance and SO 2 concentration in a defined region of the image pointed by the Ocean-Optic USB2000+ Spectrometer [47]. The use of the UV spectrometer allows us to quantitatively measure the full UV spectrum, and then fit the theoretical SO 2 absorption cross-section [48] with the differential absorption between two consecutive spectra (acquired every 5 seconds). The Ocean-Optic USB2000+ Spectrometer in use has on-board a Sony ILX511B Linear Silicon CCD Array Detector at 2048 pixels, with a wavelength response of 200-1100 nm, a dynamic range of 8.5 × 108, and a SNR of 250:1 at full signal. Calibrated SO 2 column densities over the entire images are then obtained by integrating images achieved by the UV camera with information achieved by the spectrometer. The instrument module is powered with a 12 V power supply, and requires 15 W in a fully operational mode. A Fujitsu RX100 Workstation, connected to the instrument module, automatically acquires synchronous data from the instrument module, and processes data without the need of the operator. To do this, we designed algorithms to control acquisition and processing parameters, such as the automatic tuning of camera's and spectrometer's exposures, and automatic evaluation of optimal viewing conditions (see Section 2.3.2). The computer internal time drift is controlled by a specifically designed application that reads the time-stamp from an NMEA (National Marine Electronics Association) standard message coming from a GPS antenna. The instrument module communicated with the acquisition/processing module via wired or wireless TCP/IP connection. This UV camera system was installed at the Montagnola site (on Etna Volcano, Figure 1), and designed to stream real-time SO 2 flux results using a Wi-Fi data link. The objective is to capture SO 2 emissions associated with diverse volcanic processes and dynamics, including quiescent (passive) degassing, explosive eruptions (strombolian activity/lava fountaining), and effusive eruptions [49]. Montagnola is located at~3 km distance from the active summit vents and grants perfect views of the southern sector of the summit crater area (Figure 1).

Seismic and Thermal Data
We compared our SO 2 fluxes with other independent geophysical parameters such as tremor and thermal radiance primary to have a benchmark in the calibration of automatic SO 2 flux calculation algorithm. Last but not least, we combined these data to characterize volcanic activity.
We used seismic data recorded by the ETN station [42,50], located at Lapide Malerba, at 5 km from the summit area. ETN is equipped with a broad-band seismometer (Guralp CMG-40T, with a sensitivity of 800 V/(m/s) and eigen period of 30 s). The link between SO 2 flux and volcanic tremor at Mt. Etna [22,51,52] suggests that the tremor is generated by the same degassing dynamics. We then calculated volcanic tremor amplitude from raw traces recorded at ETN station, by averaging within a 1-minute length window the maximum RMS amplitude taken within a 1-s window.
Thermal remote sensing offers a great opportunity to follow volcanic unrest from ground and space to characterize volcanic activity in near-real time [53] and to estimate near vent large pyroclastic products and lava flow discharged during eruptions. We used satellite data from the MIROVA system [43,54] and ground-based thermal cameras [44] to constrain onset, duration, and intensity through the time of eruptive events occurring at Etna during 2016. In particular, MIROVA uses the data provided by the MODIS sensor, which acquires four images per day (two daytime and two nighttime) with a spatial resolution of 1 km in the infrared bands [43]. The heat flux retrieved from MIROVA data, called Volcanic Radiative Power (VRP), is a combined measurement of the area and the integrated temperature of the hot emitters (hot vents, lava flows, etc.) at the time of each image acquisition. These data are provided without correction for cloud cover and satellite view geometry [43]. These factors may introduce noise in the dataset, which has been proven not to affect the general trend associated with the activity of Mt. Etna [43]. A mask 5 × 5 km around the volcano summit is used to filter out thermal anomalies due to wildfires.

The Algorithm for Automatic Processing of UV Camera Data
Manual processing of UV camera results [33] has the advantage that images are checked and validated manually by an operator. However, manual procedures are time-consuming, especially when dealing with huge data flows from permanent monitoring stations. To overcome this limit, we designed a MATLAB-based algorithm to automatically process images, and thus obtain SO 2 fluxes in near real-time. Using a "low-cost" PC-workstation (with 8 GB RAM and Xeon E3 type CPU), we can successfully process data at~5× speed (i.e., a 5 minutes-long record is processed in~1 min), which allows nearly real-time monitoring. The automatic routine has been calibrated using the results of the manual procedure and includes: (1) automatic determination of background absorbance levels, (2) automatic determination of image goodness, using image quality indexes, (3) estimation of gas plume speed and its distribution across the crater area, and (4) calculation of SO 2 flux distribution throughout the summit crater area. These are described in more detail below.

Image Processing and SO 2 Column Densities
Relative UV absorption by volcanic plume SO 2 is quantified by applying the Beer-Lambert law. Sets of synchronous images, taken by the two co-aligned cameras using different filters, are combined to obtain single absorbance images. This method, known as the "double filter method" [12,45], implies the use of two cameras with different filters, with one centered at 310 nm and the other centered on 330 nm, where UV radiation is/is not absorbed, respectively. The use of the two filters method allows compensating for aerosol attenuation/backscattering, to avoid any temporal mismatch associated with filter change while using a single camera, and to maintain the sampling rate at up to 0.5 Hz (two synchronous images every 2 s taken by two cameras [45]).
In our automatic routine, once a new raw image is acquired, a first quality check is made by the system, in order to keep the best exposure times to compensate any subtle changes in sunlight intensity. An automatic real-time tuning of exposure time is also applied to the UV spectrometer data, in order to obtain the best measurement dynamics within the UV bandwidth.
Synchronous images from the two cameras are then real-time corrected for vignetting effects associated with filters and optics, and normalized for relative exposure times. Residual intensities are then combined to obtain an un-calibrated absorbance image using the Lambert-Beer Law equation.
where A is the absorbance, I 310 and I 330 are pixel intensities associated with cameras mounting the 310 or 330 nm filter, while A 0 is the absorbance level associated with a clear background sky sub-area of the image (assumed to be unaffected by SO 2 absorption, I 0 310 and I 0 330 of Figure 2) and calculated as −log 10 (I 310 /I 330 ), following Kern [55].
In the automatic processing module, this background sky sub-area is automatically selected for each image by monitoring a distal sky horizontal section with respect to the vent position, and selecting the sector with the lowest absorbance intensity.
Residual absorbance is then converted into SO 2 column density integrating data from the co-located ultraviolet spectrometer, which is pointing to a known sub-area within the camera field of view. This procedure yields, in real-time, the proportionality ratio between absorbance and SO 2 column densities using the method described in McGonigle [47].  August. (f) Picture taken from Reference [56] showing BN crater collapse occurred on 10 October, and marking the end of enhanced degassing activity (see text).
Remote Sens. 2019, 11, x FOR PEER REVIEW 6 of 22 Figure 2.Calculation of the SO2 column densities using the dual camera method with 330 nm (a) and 310 nm (b) optical filters. Absorbance image (c) is obtained using the Lambert-Beer equation after image normalization with respect to background absorbance intensities (black circles in (a) and (b), automatically located). Black circles in (a) and (b) represent the sky areas where absorbance is assumed to be SO2 free (minimum absorbance level). Black rectangles (Ssky Sgrn) in (a) represent the areas used for calculation of the visibility index. The red circle in (c) shows the FOV of co-located UV-scanning spectrometer used to convert un-calibrated absorbance intensities into SO2 column densities.

Automatic Determination of the Optimal Viewing Condition
The optimal plume viewing conditions, and the presence of a clear sky, are required for reliable SO2 density measurements. However, weather conditions are extremely variable on Etna's summit, and often prevent optimal SO2 observation. To minimize uncertainties due to poor weather conditions, we set-up a sub-routine for real-time calculation of two visibility indexes. The first visibility index (Fog index) is calculated as the unsigned ratio between the mean pixel intensity associated with the camera FOV's portions capturing sky and ground, respectively. Tests we conducted on real and synthetic images show that the higher this ratio is, the better the visibility condition is. SO2 measurements are then selected by setting a threshold on the visibility index, and discarding measurements below the threshold (e.g., biased by a poor visibility condition).
Detecting a "sky" signal well above the "ground" signal ( Figure 2) is a required but not sufficient condition for reliable SO2 measurements. This is especially true in the presence of a highly condensed plume, where the SO2 absorbance signal can be masked. Thus, a second automatic procedure was developed and run in real-time, which allows us to select only images with a clear SO2 signal above atmospheric noise. This latter procedure is based on the principle of combining absorbance and 310 nm images associated with the plume. A well detected and measurable SO2 signal requires that lower intensities in the 310 nm image are measured in the plume relative to its surroundings (because SO2is absorbing solar radiation), and that higher intensities are consistently obtained in the absorbance image (see the Lambert Beer equation). This condition is only verified if the SO2 signal is high enough to emerge above the atmospheric noise. To discriminate this condition in real-time, we defined a correlation index ( Figure 3) as the correlation coefficient between absorbance and 310-nm pixel intensities over a cross-section intersecting the plume (Figures 3 and 4). The correlation coefficient is defined as C(i,j)/((C(i,i)*C(j,j))^(1/2), where C is the covariance matrix, i and j are pixel intensities over the cross-section of absorbance and 310-nm images, respectively. In such plots, the closer the correlation coefficient is to the value -1, the more absorbance can be related to gas. Images that do not satisfy this condition (e.g., that have a Correlation Index < -0.5) are disregarded by the automatic computation ( Figure 4b).

Figure 2.
Calculation of the SO 2 column densities using the dual camera method with 330 nm (a) and 310 nm (b) optical filters. Absorbance image (c) is obtained using the Lambert-Beer equation after image normalization with respect to background absorbance intensities (black circles in (a) and (b), automatically located). Black circles in (a) and (b) represent the sky areas where absorbance is assumed to be SO 2 free (minimum absorbance level). Black rectangles (Ssky Sgrn) in (a) represent the areas used for calculation of the visibility index. The red circle in (c) shows the FOV of co-located UV-scanning spectrometer used to convert un-calibrated absorbance intensities into SO 2 column densities.

Automatic Determination of the Optimal Viewing Condition
The optimal plume viewing conditions, and the presence of a clear sky, are required for reliable SO 2 density measurements. However, weather conditions are extremely variable on Etna's summit, and often prevent optimal SO 2 observation. To minimize uncertainties due to poor weather conditions, we set-up a sub-routine for real-time calculation of two visibility indexes. The first visibility index (Fog index) is calculated as the unsigned ratio between the mean pixel intensity associated with the camera FOV's portions capturing sky and ground, respectively. Tests we conducted on real and synthetic images show that the higher this ratio is, the better the visibility condition is. SO 2 measurements are then selected by setting a threshold on the visibility index, and discarding measurements below the threshold (e.g., biased by a poor visibility condition).
Detecting a "sky" signal well above the "ground" signal ( Figure 2) is a required but not sufficient condition for reliable SO 2 measurements. This is especially true in the presence of a highly condensed plume, where the SO 2 absorbance signal can be masked. Thus, a second automatic procedure was developed and run in real-time, which allows us to select only images with a clear SO 2 signal above atmospheric noise.
This latter procedure is based on the principle of combining absorbance and 310 nm images associated with the plume. A well detected and measurable SO 2 signal requires that lower intensities in the 310 nm image are measured in the plume relative to its surroundings (because SO 2 is absorbing solar radiation), and that higher intensities are consistently obtained in the absorbance image (see the Lambert Beer equation). This condition is only verified if the SO 2 signal is high enough to emerge above the atmospheric noise. To discriminate this condition in real-time, we defined a correlation index ( Figure 3) as the correlation coefficient between absorbance and 310-nm pixel intensities over a cross-section intersecting the plume (Figures 3 and 4). The correlation coefficient is defined as where C is the covariance matrix, i and j are pixel intensities over the cross-section of absorbance and 310-nm images, respectively. In such plots, the closer the correlation coefficient is to the value −1, the more absorbance can be related to gas. Images that do not satisfy this condition (e.g., that have a Correlation Index < −0.5) are disregarded by the automatic computation ( Figure 4).   Gas is visible only when the fog index is greater than 4 and the correlation index is less than -0.5.

Plume Velocity Field
A robust plume velocity field is mandatory for reliable SO2 flux measurements. Errors in plume velocity have shown to contribute to 40% or more of the overall error in the determined fluxes [13,57].
The UV camera approach offers the unique opportunity to track the gas while dispersing right after atmospheric emission, which minimizes errors in plume speed determination of yet more established DOAS and COSPEC methods. These methods indirectly infer plume speed from either on-site measurement of wind velocity or from meteorological models [11,15].
The UV camera approach allows us to derive the velocity profile over the summit craters by applying an optical flow algorithm that tracks gas fronts in consecutive frames [58].   Gas is visible only when the fog index is greater than 4 and the correlation index is less than -0.5.

Plume Velocity Field
A robust plume velocity field is mandatory for reliable SO2 flux measurements. Errors in plume velocity have shown to contribute to 40% or more of the overall error in the determined fluxes [13,57].
The UV camera approach offers the unique opportunity to track the gas while dispersing right after atmospheric emission, which minimizes errors in plume speed determination of yet more established DOAS and COSPEC methods. These methods indirectly infer plume speed from either on-site measurement of wind velocity or from meteorological models [11,15].
The UV camera approach allows us to derive the velocity profile over the summit craters by applying an optical flow algorithm that tracks gas fronts in consecutive frames [58]. . Gas is visible only when the fog index is greater than 4 and the correlation index is less than −0.5.

Plume Velocity Field
A robust plume velocity field is mandatory for reliable SO 2 flux measurements. Errors in plume velocity have shown to contribute to 40% or more of the overall error in the determined fluxes [13,57].
The UV camera approach offers the unique opportunity to track the gas while dispersing right after atmospheric emission, which minimizes errors in plume speed determination of yet more established DOAS and COSPEC methods. These methods indirectly infer plume speed from either on-site measurement of wind velocity or from meteorological models [11,15].
The UV camera approach allows us to derive the velocity profile over the summit craters by applying an optical flow algorithm that tracks gas fronts in consecutive frames [58].
Optical flow consists ofthe apparent motion pattern of image objects between two consecutive frames, caused by the movement of either the object or the camera, and is valid under the assumptions that the pixel intensities of an object do not change significantly between consecutive frames, and that the neighboring pixels have similar motion.
If a pixel, with intensity I(x,y,t), where (x,y) are the pixel coordinates and t is the time in first frame, moves by distance (dx,dy) in the next frame taken after time dt, it can be assumed that: Then, from the Taylor series approximation of the right-hand side, removing common terms and dividing by dt, one gets the following equations.
where ∂I ∂x and ∂I ∂y are the image gradients, ∂I ∂t is the gradient along time, and u and v are horizontal and vertical velocities that are unknown. Lucas & Kanade [59] provide a method (LK) to derive these unknown velocities, by solving the basic optical flow Equations (3) and (4) for all the pixels using the least squares criterion, and by combining information from nearby pixels. We then applied the LK algorithm, included within the Open-CV toolbox, to our dataset [60]. We tested the performance of this method by applying it to artificial images with known particle velocities. The method has successfully determined velocity field with an error of <5%.
Absorbance images, obtained using the two filters method, contain gas-rich and ash-free portions of the plume with a higher absorbance relative to the background, and/or to ash-rich or particle-rich plume portions. We exploit this feature to track only gas moving fronts in consecutive frames by filtering them from other moving features such as lapilli and ash and by applying the LK method to the absorbance images rather than to the raw images directly acquired by our dual camera system. Velocities are then calculated by selecting the best features to track within the image, and which correspond to the areas with the highest pixel intensities (i.e., high SO 2 column densities) and high spatial coherence in consecutive frames (taken every 2 s).

Image Analysis and SO 2 Flux Calculation
To enhance the contrast between volcanic emissions and atmospheric noise, and to limit dispersion effects and chemical conversions of SO 2 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 SO 2 column densities and a plume velocity field ( Figure 2). We calculated SO 2 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 SO 2 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 SO 2 density flux (in kgm −1 s −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. 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. 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.

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].
Comparison between manually and automatically calculated fluxes ( Figure 6) demonstrated a correlation coefficient (R 2 ) 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 The SO 2 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 SO 2 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.

Validation of the Automatic Method
Validity of the automatically processed SO 2 fluxes was tested for a comparison with SO 2 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, SO 2 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]. 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].

Figure 6.
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.

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][62][63][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][69][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. Figure 6. Comparison between daily averaged SO 2 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 SO 2 flux calculation.

Etna's Activity in 2016, SO2 Flux Records, and Comparison with Seismic, Thermal Dataset, and Field Observations
Comparison between manually and automatically calculated fluxes ( Figure 6) demonstrated a correlation coefficient (R 2 ) 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 SO 2 column densities, as shown by the position of outliers ( Figure 6). Overall, this comparison validated the use of the automatic SO 2 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].

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 SO 2 flux record. SO 2 flux measurements have become fundamental in volcano monitoring to define the rates of magma ascent and degassing within the shallow (<3 km) plumbing system. SO 2 fluxes have been measured on Etna since the 1970s using the COSPEC (Correlation Spectrometer) [49,[61][62][63][64] and, more recently, a network of Differential Optical absorption spectrometers (FLAME [65,66]). In view of its recurrent activity and robust past SO 2 flux record, Mt. Etna is an ideal test site for validating our automatic processing method.
We reported below on the SO 2 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][69][70].
Our aim was to test if different SO 2 degassing regimes, related to such diverse activity styles, could be resolved and characterized in automatic (and in nearly real-time) using our permanent SO 2 camera system.

Etna's Activity in 2016, SO 2 Flux Records, and Comparison with Seismic, Thermal Dataset, and Field Observations
The SO 2 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 SO 2 flux behavior ( Figure 7). As illustrated in Figure 7, our 2016 temporal record shows daily averaged SO 2 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 SO 2 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]. 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. 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 SO 2 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. 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 SO 2 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, SO 2 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.

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][69][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 3 [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). 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. 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 SO 2 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 SO 2 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 SO 2 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 SO 2 fluxes obtained on this specific day (Figure 8), thus, reflected the "explosive" SO 2 contribution during the paroxysmal event. In contrast, the SO 2 flux peaks on 22 and 25 May could not be explained by syn-explosive SO 2 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.

Volcanic Activity after the Eruptive Phase (26 May-31 Dec)
Reduced degassing emissions (<2000 t/d) were measured after the eruptive phase, from May 26 th 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 SO 2 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 SO 2 fluxes down to low values (~1000 tons/day on average). RMS seismic amplitude and thermal radiance slightly decreased as well.

Syn-Explosive SO 2 Emissions during the Lava Fountaining Event
The 18 May lava fountaining event, entirely captured by the SO 2 camera (Figure 9), allowed us to explore to what extent SO 2 cameras could resolve the degassing dynamics associated with a lava fountaining event. Lava fountains are of special concern at Etna since the volcanic ash they inject into the atmosphere is a potential threat to aviation and population living in the surroundings [73]. These events, while very well monitored and understood [34,[36][37][38][39][40][41][42]50,[74][75][76][77], are poorly characterized in terms of their associated gas emission rates and volumes. Our 18 May results (Figure 9), therefore, represent one of the first syn-explosive gas records on the volcano [28,30].
indicated that negative peaks in the SO2 flux time-series were systematically associated with the presence of ash. The latter severely impacts SO2 detection via UV cameras [21], particularly in near-vent measurements where plumes can be very ash-rich and, thus, optically opaque. Thermal and visual observations showed that ash caused high-frequency fluctuations (short-lived negative peaks) in the SO2 flux record particularly during the paroxysm climax (~ 11.30-12.00), but also prior to the lava fountain onset, e.g., after 09.00 when the visual camera captured the first ash emission with no thermal anomaly yet detected (Figure 9). Onset of the lava fountain at 10:57 UTC, as constrained by co-acquired images of the INGV-OE thermal and visual cameras (Figure 9), was clearly detected as a visible SO 2 flux increase up to~260 kg/s (22,000 t/d), relative to a pre-eruptive level of~12 kg/s (~1000 t/d). In the following hour, while activity escalated to peak at~11.30-12.00 (see thermal records), a fluctuating and irregular SO 2 flux trend is registered ( Figure 9). Co-acquired thermal and visual images (see panels in Figure 9), clearly indicated that negative peaks in the SO 2 flux time-series were systematically associated with the presence of ash. The latter severely impacts SO 2 detection via UV cameras [21], particularly in near-vent measurements where plumes can be very ash-rich and, thus, optically opaque. Thermal and visual observations showed that ash caused high-frequency fluctuations (short-lived negative peaks) in the SO 2 flux record particularly during the paroxysm climax (~11.30-12.00), but also prior to the lava fountain onset, e.g., after 09.00 when the visual camera captured the first ash emission with no thermal anomaly yet detected (Figure 9).

Discussion
Our results demonstrated the ability of the novel automatic processing routine to capture fluctuations in the SO 2 regime, which responded to changes in Etna's activity style (Figure 7). The automatically processed SO 2 data were consistent with those manually obtained ( Figure 6). However, requiring no operator time and being obtained/delivered in nearly real-time, they represented a clear advantage for monitoring purposes.
Our 2016 UV camera-based dataset also provided novel insights into the relationship between rates/modes of SO 2 release and eruptive/degassing styles. These latter were quite diverse in 2016 since they ranged from non-eruptive quiescent degassing to intense paroxysmal activity in May. Nicely, our automatically processed SO 2 fluxes peaked during heightened activity (Figures 7 and 8) during the May 2016 eruptive sequence, and during the July-August degassing unrest that led to opening of the VOR incandescent vent on 7 August.
One important observation is the mild but significant SO 2 flux increase in April 2016, as demonstrated by a clear change in the slope of the cumulative SO 2 mass (see orange-coloured area in Figure 7a). The SO 2 flux increase started more than one month before onset of the May eruptive sequence, which is a period when the seismic tremor was at average levels and no significant thermal anomaly was detected (Figure 7). We interpreted these escalating SO 2 fluxes as reflecting the slow but systematic increase in the magma supply rate to the shallow (<3 km [78]) Etna's magma feeding system that triggered the May 2016 paroxysmal sequence [68]. However, while the relatively subtle changes in SO 2 passive degassing in April 2016 may have represented a precursory sign for the imminent (May 2016) eruption, we still noticed that a similar SO 2 increase was observed from July to September 2016 (orange-coloured sub-interval in Figure 7a). This did not culminate into an eruption yet (if not for the opening of the 7 August summit vent). Thus, the volcano's feedback to increased shallow magma emplacement (as indicated by increasing SO 2 fluxes) may be different from time to time, perhaps in response to distinct stress regimes prevailing on the upper part of the edifice, and/or temporally varying feedbacks between magma ascent and rates of gravitational spreading of the mobile eastern flank [79,80].
We also characterized SO 2 emissions associated with a lava fountaining event at a high spatial and temporal resolution [28,30]. Even though ash is a serious issue for ground-based SO 2 remote sensing during explosive eruptions, we still noted that the entirety of 18 May eruptive episode was well marked by elevated SO 2 fluxes well above background emissions in the pre-event and post-event phases. In fact, SO 2 emissions manifestly dropped down at only 14:00 UTC, when declining thermal emissions consistently marked lava fountaining termination. From such, we found it useful to tentatively estimate the cumulative SO 2 mass released by the 18 May lava-fountain, by integrating the signal over the eruption duration. We obtained a total SO 2 mass released in the event of~1700 tons, and an average flux of 158 kg/s (13,600 t/d) for a total duration of~3 h. We caution this inferred mass corresponds to a lower range estimate, due to the presence of volcanic ash that severely depressed the measured SO 2 signal during the eruption climax. It is well established that SO 2 fluxes are directly linked to the rate of magma ascent and degassing [78]. Thus, temporal variations in SO 2 fluxes do reflect changes in magma feeding to the volcano's shallow plumbing system, and, as such, may help track transition in activity style, from quiet passive degassing to eruptive periods [81]. On Etna, we now have 3 years of UV camera observations available ( [28,30], this study), during which transition from quiescence to eruption has frequently been observed. We, thus, examined our dataset in the attempt to tentatively identify any possible systematic SO 2 flux threshold/trend corresponding to such an activity switch.
For this purpose, in Figure 10, we compared the cumulative SO 2 flux trends (time-normalized) for three different periods encompassing three eruptive paroxysmal episodes, which occurred in August 2014, December 2015, and May 2016, respectively. For each of three events, we calculated the averaged SO 2 fluxes (corresponding to the slopes of the cumulative curve) in the periods before, during, and after eruption, and we found significant similarities between the three events. In each of the three 2014-2016 events, the pre-eruptive fluxes fell in a relatively narrow range, between 1900 and 2500 t/d. The syn-eruptive (during the paroxysmal sequence) fluxes were typically higher, and spanned between 3000 and 5200 t/d (Figure 10), and were the highest during the December 2015 paroxysmal sequence that was consistently the most energetic in the past few years [67]. Lastly, each of the three post-paroxysmal phases was characterized by reduced SO 2 emissions, ranging between 600 and 900 t/d, which impliesa reduced magma supply and degassing of a volatile depleted (residual) magma after each eruptive episode. These ranges were fully in agreement with those indicated by Reference [49] and by the analysis of a 13-year-long dataset. and 2500 t/d. The syn-eruptive (during the paroxysmal sequence) fluxes were typically higher, and spanned between 3000 and 5200 t/d (Figure 10), and were the highest during the December 2015 paroxysmal sequence that was consistently the most energetic in the past few years [67]. Lastly, each of the three post-paroxysmal phases was characterized by reduced SO2 emissions, ranging between 600 and 900 t/d, which impliesa reduced magma supply and degassing of a volatile depleted (residual) magma after each eruptive episode. These ranges were fully in agreement with those indicated by Reference [49] and by the analysis of a 13-year-long dataset.
Our preliminary results are suggestive of the existence of a systematic pattern in SO2 emissions that, if confirmed, would imply a recurrent degassing process/mechanism prior to, during, and after the Etna's eruptive periods. Clearly, additional data are required to corroborate this initial hypothesis. Pre-eruptive, syn-eruptive, and post-eruptive phases show similar SO2 fluxes for all these three events. In each of these three events, pre-eruptive fluxes range between 1900 and 2500 t/d, syn-eruptive fluxes are the highest (3000-5200 t/d), while post-eruptive fluxes are systematically the lowest (<900 t/d). Our preliminary results are suggestive of the existence of a systematic pattern in SO 2 emissions that, if confirmed, would imply a recurrent degassing process/mechanism prior to, during, and after the Etna's eruptive periods. Clearly, additional data are required to corroborate this initial hypothesis.

Conclusions
SO 2 imaging at Mt Etna during 2016 revealed different styles of gas emissions, which reflected changes in volcanic activity, from quietly passive degassing to eruptive activity (lava fountaining, intense strombolian activity, and lava flowing).
To real-time characterize and monitor this very dynamic activity period, we designed a novel routine to automatically calculate SO 2 fluxes including computer-based detection of image quality, and calculation of plume speed time-series using computer vision libraries. This automatic processing routine allowed us to obtain real-time information on volcano degassing dynamics at high spatial and temporal resolution, which results in a further step in instrumental volcanic gas monitoring. We validated the methodology through a comparison with manually processed results and integration with independent thermal and seismic observations. All these independent datasets showed coherent temporal variations that validated the use of UV cameras for detecting subtle changes in volcanic and degassing activity. Our novel method, thus, promises a step ahead in instrumental volcanic gas monitoring.
Our automatically derived SO 2 flux time-series were used to constrain degassing regimes on Mt. Etna in 2016. We have shown that our automatically processed SO 2 fluxes peaked during heightened activity, such as during the May 2016 eruptive sequence, and during a phase of elevated degassing in July-August 2016, which culminated with the opening of a new summit incandescent vent. The May 2016 sequence was preceded by a circa one-month-long phase of a mild but detectable SO 2 flux increase, and was followed by an abrupt drop in degassing. Comparison between our SO 2 flux time-series in May 2016 and those associated with two other eruptive paroxysmal sequences (occurred in 2014 and 2015) highlighted strong similarities in SO 2 flux dynamics, which implies the possible existence of thresholds that distinguish between degassing regimes, prior to, during, and after eruptions. Pre-paroxysm SO 2 fluxes were found to have consistent values (of~2000 t/d) during the three episodes. Similarly, the highest SO 2 fluxes (from 3000 t/d up to 5200 t/d on a daily average basis) were identified during the three eruptive sequences, while post-eruptive were systematically characterized by reduced degassing (<1000 t/d). This result, if confirmed by future observations, may bring implications for identifying switches in volcanic activity regime. We believe this methodology can be successfully exported on other open-vent active volcanoes, after a relatively brief period of calibration.
We also tested the ability of UV camera records to characterizing SO 2 emissions during a lava fountain episode. This is, to the best of our knowledge, one of the first SO 2 camera record at high spatial and temporal resolution of an ongoing lava fountaining [28,30]. We have reported on high levels of degassing during the lava fountain (of up to 260 kg/s), from which we assessed a lower limit (due to the ash presence within the plume) for the cumulative SO 2 mass of~1700 tons emitted during a lava fountain episode.