**3. Methods: RSTVOLC Algorithm**

The RSTVOLC algorithm [28] identifies volcanic hotspots by means of two local variation indices defined as:

$$\log\_{MIR}(\mathbf{x}, \mathbf{y}, t) = \frac{B T\_{MIR}(\mathbf{x}, \mathbf{y}, t) - \mu\_{MIR}(\mathbf{x}, \mathbf{y})}{\sigma\_{MIR}(\mathbf{x}, \mathbf{y})} \tag{1}$$

$$\otimes\_{MIR-TIR}(\mathbf{x}, \mathbf{y}, \mathbf{t}) = \frac{\Delta T(\mathbf{x}, \mathbf{y}, \mathbf{t}) - \mu\_{\Delta T}(\mathbf{x}, \mathbf{y})}{\sigma\_{MIR}(\mathbf{x}, \mathbf{y})} \tag{2}$$

In Equation (1), *BTMIR*(*x*, *y*, *t*) is the MIR (medium infrared) brightness temperature measured at the time *t* for each pixel (*x*, *y*), whereas *μMIR*(*x*, *y*) and *σMIR*(*x*, *y*) are the relative temporal mean and standard deviation. In Equation (2), Δ*T*(*x*, *y*, *t*) = *BTMIR*(*x*, *y*, *t*) − *BTTIR*(*x*, *y*, *t*), where *BTTIR*(*x*, *y*, *t*) is the brightness temperature measured in the TIR (thermal infrared) band at around 11 μm wavelength; *μ*Δ*T*(*x*, *y*) and *σ*Δ*T*(*x*, *y*) stand for the relative temporal mean and standard deviation. These terms are calculated after processing cloud-free satellite records selected according to specific homogeneity criteria (i.e., same spectral channel/s, calendar month, satellite overpass time). In particular, the OCA (One Channel Cloud-Detection Approach) RST-based method [38] is generally used to filter out cloudy pixels from the scenes. The iterative *kσ* clipping filter, which is also implemented within the RSTVOLC process, enables the removal of signal outliers (e.g., extremely hot pixels) [29].

The  *MIR*(*x*, *y*, *t*) index identifies anomalous signal variations in the MIR band of sensors like AVHRR (channel 3: 3.55–3.93 μm) and MODIS (channels 21/22: 3.929–3.989 μm), where hot magmatic surfaces reach the peak of thermal emissions [39,40]. The  *MIR*−*TIR*(*x*, *<sup>y</sup>*, *<sup>t</sup>*) index is used jointly with the previous one for minimizing spurious effects associated with non-volcanological signal fluctuations [28]. RSTVOLC, combining those indices, is capable of guaranteeing an efficient identification of volcanic thermal anomalies (a cloud-masking procedure is used in the daytime before running the algorithm) under different observational conditions (e.g., [12,18,41–43]).

#### **4. Results**

#### *4.1. Monitoring the Paroxysmal Events of May 2016*

To investigate the Mt. Etna eruptive events of May 2016, we show, in Figure 4, the curve of the volcanic tremor (top panel) and the time series of the radiant flux (*Qrad*) retrieved from AVHRR and MODIS data uncorrected for atmospheric effects (bottom panel). In more detail, in Figure 4b, we also show eruption chronology and overcast periods, for better assessing the impact of clouds on the achieved results. We estimated the radiant flux from infrared AVHRR data by outputs of a dual band-three components method [44]; the mean *Qrad* value calculated from two end-members was considered. The formulation proposed by [45] and amended by [46] was used to retrieve the same parameter from MODIS data.

**Figure 4.** (**a**) Root-mean-square (RMS) tremor amplitude at Etna Cratere del Piano (ECPN) seismic station, a.u. = arbitrary unit; (**b**) radiant flux retrieved from Advanced Very High Resolution Radiometer (AVHRR) (orange bars) and Moderate Resolution Imaging Spectroradiometer (MODIS) (blue bars) data from the 17 to 31 May. Eruption chronology and information about overcast periods (in light grey), provided by the One channel Cloud-Detection Approach (OCA) method, is also reported.

Figure 4b shows that the radiant flux estimated by satellite ranged, during the period of interest, from a few MW up to a few GW. In particular, the paroxysmal event of May 18 was the most intense, leading to *Qrad* values, retrieved from daytime MODIS and AVHRR data acquired under comparable values of satellite zenith angle (SZA), in the range 2.5–4.0 GW. This paroxysm, occurring with lava fountaining from VOR and a lava overflow from BN, generated the major peak in the volcanic tremor, see Figure 4a. The following Strombolian activity lead to another sudden increase of this parameter but it was undetected by RSTVOLC because of clouds, see Figure 4b.

During 20–21 May, another significant increase of *Qrad*, up to about 2.0 GW, was recorded. This increase of radiant flux was determined by an eruptive activity similar to that of 18 May (see eruption chronology), leading to the third abrupt increment in the volcanic tremor. The latter increased, although in a less significant way than before, on 24 May when another lava fountaining activity occurred. Based on the retrieved *Qrad* values, this eruptive episode was slightly less intense than the previous one. Figure 4b shows that RSTVOLC also identified the Strombolian event of 22–23 May. Moreover, it detected a thermal anomaly in between the Strombolian activities of 21 May and 22 May, as indicated by values of the radiant flux retrieved in that period which were mostly lower than 1 GW. In addition, an abrupt increase of *Qrad* was recorded in short-time intervals (e.g., in between AVHRR and MODIS observations of 24 May at 20:23 UTC and 21:12 UTC), revealing some abrupt variations in the intensity of the volcanic thermal emissions. However, due to cloud coverage and satellite overpass times, some thermal activities (e.g., see 24 May) were undetected from space. On the other hand, even the discontinuous identification of a low-level thermal anomaly after 25 May, possibly associated to a weak degassing from VOR (leading to a minimum *Qrad* value of about 7.5 MW), is ascribable to clouds.

#### *4.2. Investigating the Thermal Activity of June–August 2016*

In Figure 5, we show the temporal trend of the total MIR radiance retrieved from nighttime AVHRR and MODIS records of June–August 2016. We investigated only data with relatively low values of satellite zenith angle (i.e., SZA < 40◦), unlike the eruptive events of May 2016, for minimizing the impact of the satellite viewing geometry on the thermal anomaly detected from space.

**Figure 5.** Total medium infrared (MIR) radiance retrieved from hotspot pixels identified by Robust Satellite Techniques–volcanoes (RSTVOLC) on nighttime AVHRR (blue bars) and MODIS (red bars) data of June–August 2016. Dotted green line indicates the start of hot degassing activity at VOR. Note that no MODIS data were acquired at Institute of Methodologies for Environmental Analysis (IMAA) during 1−3 July because of antenna problems; AVHRR data acquired in July were unprocessed due to geo-location issues.

During the first week of June, no MODIS data were acquired at IMAA because of some antenna problems. In July, AVHRR records were unprocessed due to some geolocation issues (the system was back to being fully operational in early August). Despite these limitations, reducing the number of available satellite scenes, RSTVOLC provided information about changes in Mt. Etna thermal activity.

Indeed, after the identification of a sporadic thermal anomaly in June RSTVOLC revealed, since early July, the occurrence of more continuous thermal emissions at the crater area. In more detail, Figure 5 shows that after 4 July, the total MIR radiances decreased, because of a less significant thermal anomaly (see Figure 6), showing small fluctuations until 3 August when an evident increase of this parameter was recorded. During 10–12 August, detected thermal anomaly was more extended in terms of hot spot pixels, as shown by RSTVOLC maps which are not reported here, leading to peak of total MIR radiance.

**Figure 6.** (**a**) MODIS channel 22 (MIR) images with an indication (see white arrow) of the volcanic thermal anomaly identified by RSTVOLC (brighter tones indicate higher brightness temperature values); 4 July 2016 at 21:06 UTC (BT22MAX = 297.65 K); (**b**) 6 July 2016 at 20:53 UTC (BT22MAX = 290.29 K).

(**a**) (**b**)

To assess RSTVOLC detections, as well as changes of thermal activity revealed by Figure 5, we analysed the cloud-free Sentinel 2A-MSI (Multispectral Instrument) data of June–August 2016 made freely available online by the Sentinel Hub.

The MSI has 13 spectral channels centered in the VNIR and SWIR bands having a different spatial resolution [47]. In Figure 7, we show a number of false color composite images, nine in total from 2 June to 28 August, magnified over the Mt. Etna crater area. These RGB (red-green-blue) products, at a 20-m spatial resolution, were generated by using bands 12 (2.19 μm), 8A (0.865 μm) and 4 (0.665 μm). The figure shows that some crater pixels assumed the red color owing to the dominance of the SWIR component.

Since high-temperature magmatic surfaces are highly radiant even in the SWIR region (e.g., [48]), and considering that other features (e.g., clouds) are clearly recognizable in Figure 7, the presence of a thermal anomaly at the Mt. Etna crater area, before and after the start of hot degassing activity at VOR, was confirmed. Besides, it can be noted that some VOR pixels were brighter, especially, on satellite scenes of 2 July (i.e., two days before the evident increase of total MIR radiance revealed by Figure 5) and 21 August (see Figure 7).

The false color composite imagery, at a 30-m spatial resolution, of Figure 8 generated from Landsat 8-OLI (Landsat 8-Operational Land Imager) data of 5 July and 6 August provided by UGSG (United States Geological Survey), further corroborate information provided by RSTVOLC. The figure confirms that a thermal anomaly affected the Mt. Etna crater area before the opening of the new degassing vent; see pixels, in red, highly radiant in OLI band 7 (2.1–2.3 μm).

**Figure 7.** RGB (Red= Band 12, 2.19 μm; Green = Band 8A, 0.865 μm; Blue = Band 4, 0.665 μm) Sentinel-2 products, at 20 m spatial resolution, of June–August 2016 generated from Level 1C data excluding images completely overcast over target area. The white arrow indicates pixels more radiant in the shortwave infrared (SWIR) band at VOR.

**Figure 8.** RGB (Red = Band 7, 2.1–2.3 μm; Green = Band 5, 0.845–0.885 μm; Blue = Band 2, 0.450–0.515 μm) Landsat 8-Operational Land Imager (OLI) images, at a 30-m spatial resolution, showing hotspot pixels (in red), magnified at the top-right side of each panel, affecting the Mt. Etna area. (**a**) 5 July 2016; (**b**) 6 August 2016.

To localize, in a more accurate way, the source of the thermal anomaly identified by the satellite, we analyzed the temporal trend of the SWIR radiance. In more detail, we used the dark object subtraction (DOS) method [49] to correct SWIR radiance for effects of atmospheric scattering by subtracting the radiance value of the darkest object from each pixel of the image (e.g., [50]). We focused our analyses on the VOR area since the thermal anomaly detected by RSTVOLC generally did not include the NSEC; see AVHRR pixel polygon in green overlapped with the OLI sub-scene in the inset

of Figure 9. The latter shows the similar behavior of SWIR curves at nine VOR pixels, apart from 5 July when the SWIR signal decreased only at VOR5 and VOR7.

**Figure 9.** Temporal plot of the OLI-SWIR radiance retrieved, after applying the dark object subtraction (DOS) method, over nine different VOR pixels. In the inset, the OLI sub-image of 6 August, with indication of most radiant pixel (VOR2), and the AVHRR pixel polygon (in dotted green) associated to thermal anomaly flagged by RSTVOLC before the opening of the new degassing vent.

In addition, the plot shows that the strongest increase in the SWIR signal was recorded at VOR2 on 6 August, see the orange curve. By the temporal trend of the total SWIR radiance retrieved, pixel by pixel, along the A-B transect region of Figure 10, we found that VOR2, whose geographic coordinates of the center are reported in Figure 9, was the most radiant pixel in OLI band 7 (see value of the analyzed parameter retrieved in correspondence to the dotted red line). Hence, it is reasonable to suppose that the volcanic thermal emissions from the VOR2 area affected, in a more significant way, the thermal anomaly detected by RSTVOLC after eruptions of May and before the start of high-temperature degassing activity on 7 August.

**Figure 10.** Spatial profile of the total SWIR radiance retrieved along the A–B transect region intersecting the VOR area (see OLI data sub-scene in the inset of the figure). The most radiant pixel, whose location is indicated by the dotted red line, having geographic coordinates of the center 37◦45 6.08"N, 14◦59 45.08"E, corresponds to VOR2 in Figure 9.

#### **5. Discussion**

In this paper, we have investigated the Mt. Etna activity of May–August 2016, exploiting the information provided by the satellite-based monitoring system operating at IMAA (whose products are currently used only for research purposes) implementing the RSTVOLC algorithm. The latter runs on both AVHRR and MODIS data, which currently guarantee more than 10 observations per day over Italy, increasing the probability of processing cloud-free satellite scenes [51].

The results of this study confirm that RSTVOLC performed in a similar way when detecting thermal anomalies regardless of the satellite data used; i.e., despite the different features of AVHRR and MODIS instruments. Indeed, by combining values of the radiant flux, it was possible to investigate eruptions occurring in May in a more continuous way than using data from a single satellite sensor; although, a different retrieval method was used. Indeed, despite the possible inaccuracy in estimating the radiant flux due to the analysis of satellite records which were uncorrected for atmospheric effects, we observed a good agreement between temporal changes of *Qrad* and information provided by field observations. High values of radiant flux characterized periods of lava effusion/fountaining, whereas low values of the same parameter were retrieved in the presence of Strombolian eruptions, as shown in Figure 4.

Regarding the few discrepancies with ground-based observations, they were mostly related to clouds. Although recent literature studies performed by means of ground-based hyperspectral imagers have suggested that these features could be accounted for under certain circumstances [45], clouds generally represent a common issue for satellite-based methods developed for monitoring thermal volcanic activity. In this work, clouds had a not negligible impact on the results of the thermal anomaly detection because, on some days, they partially or completely obscured the underlying hot surfaces. Consequently, some Strombolian activities were identified by satellite some hours after the eruption onset; RSTVOLC detected, for instance, the first thermal anomaly over the Mt. Etna area on 18 May, on AVHRR overpass of 01:17 UTC (universal time coordinated); although, a Strombolian activity was already in progress the day before, see eruption chronology in Figure 4. Moreover, a few daytime MODIS scenes strongly affected by clouds, which were not completely removed by the OCA method, showed artifacts. These features were recognized and filtered out from the analyzed time series.

Despite the impact of clouds and satellite viewing geometry on the thermal anomaly identification (e.g., as for AVHRR data acquired at 17:05 UTC), this study provides additional information concerning eruptive events occurring in May. Specifically, Figure 2 showed that, in spite of possible lava cooling effects, the eruptive activity of 21 May probably lasted longer than was indicated by field observations. Although, the occurrence of a short-lived eruptive event in the late afternoon of the same day cannot be excluded. The unfavorable meteorological conditions possibly affected the direct/visual observation of these phenomena. Concerning the short-term variations of *Qrad* also revealed in Figure 2, they require further investigation to be assessed, analyzing, for instance, infrared SEVIRI data.

RSTVOLC detections performed after the eruptive events of May revealed that the intensity of fumarolic emissions from the fracture fields cutting the entire summit area changed before 7 August, possibly due to a preparatory phase of hot degassing activity occurring at VOR. This is even more evident from looking at the curve of the radiant flux (in blue), retrieved from nighttime MODIS data, in Figure 11. The plot shows that *Qrad* was around 30 MW on 4 July, then decreased in the following weeks, when its average value was ~7.5 MW, and increased up to 18 MW on 3 August (see right axis); i.e., four days before the opening of the new degassing vent marked by the dotted black line. After a slight reduction, the radiant flux once again increased on 7 August owing to the high-temperature degassing activity in progress at VOR, reaching its peak of around 33 MW three days later, see Figure 11.

**Figure 11.** SWIR radiance (left axis), from OLI band 7 (red points) and MSI band 12 (green triangles) data of June–August 2016, measured over VOR2 area and corrected using the DOS method (see Section 4.2). In blue, the temporal trend of radiant flux (right axis) retrieved after filtering MODIS data for values of SZA < 40◦. The dotted black line indicates the start of hot degassing activity at VOR.

Those variations of radiant flux were consistent with information independently provided by OLI and MSI data at the VOR2 location, corresponding to the area where the new degassing vent would later open, see the previous section. Figure 11 shows that SWIR radiance measured at the VOR2 crater area was higher in early July compared to mid-end July, and significantly increased on 6 August (see left axis), strengthening the hypothesis formulated above.

The information retrieved by satellite is perfectly compatible with the actual evolution of the system of degassing fractures, which progressively opened in the summit area of Mt. Etna. It should be pointed out that, at the end of May 2016, eruptions from all summit craters had conduits almost totally obstructed by the erupted volcanic products, and, therefore, they were no longer capable of degassing efficiently. Under these conditions, it is very likely that the pressure of the volcanic gases that rise inside the obstructed conduits have exerted a considerable pressure, also causing the heating and the thermal alteration of the rocks in the immediate vicinity of the conduits themselves, looking for alternative outputs. Therefore, the fracture system opening in the summit area could have been an alternative way to allow the volcanic gas to reach the surface. At first, this gradual fracturing process triggered the thermal anomalies observed from space during the first days of July. Finally, the opening of the degassing vent on August 7 close to the VOR's rim could be the culmination of this fracturing process, which eventually drained most of the gas from the obstructed NEC conduit.

This evolution has already been observed several times during the last twenty years at Etna. Starting from January 1998, a system of N-S fractures opened, and then progressively propagated, affecting exactly the same portion of the summit area, generating dry fractures, fumarolic activities, and eruptive fissures [33], similar to what happened in May–August 2016. Thus, the opening of the vent on 7 August was possibly triggered by the intersection of this developing fracture system with the North-East crater (NEC) plumbing system, leading to the draining of gas that had remained blocked inside its conduits after the May eruption. Finally, it should be emphasized that the opening of this fracture system seems to be triggered by an acceleration of the deformations affecting the eastern flank of the volcano [26,27,35–37]); therefore, it is possible that the opening of new degassing or erupting vents will continue to take place in this portion of the summit crater area.

### **6. Conclusions**

This study demonstrates that multi-platform observing systems may also provide a relevant contribution for monitoring active volcanoes in areas where efficient ground-based surveillance systems exist, enabling the identification of low-temperature fumarole fields whose intensity variations may precede new and more significant phases of thermal unrest.

In this context, the performance of the RSTVOLC system may be further increased using VIIRS (Visible Infrared Imaging Radiometer Suite), flying onboard Suomi NPP (Suomi National Polar-orbiting Partnership) and JPPS-1 (Joint Polar Satellite System) satellites. VIIRS collects data in 22 different spectral bands, ranging from 0.412 μm (VIS) to 12.01 μm (TIR) and including 17 spectral channels at a 750-m spatial resolution, i.e., the moderate resolution bands (M-bands) and the Day/Night panchromatic band (DNB), and five imaging resolution bands (I-bands), having a spatial resolution of 375 m at the nadir. Among the channels, the I4 band (3.55–3.93 μm) should guarantee further improvements in the identification of weak thermal anomalies (i.e., those of low temperature and/or spatial extent), as indicated by some investigations currently in progress, making RSTVOLC even more effective in detecting subtle changes of thermal volcanic activity at Mt. Etna.

**Author Contributions:** F.M., M.N. and N.P. conceived the work and wrote the paper. All the other authors contributed in analyzing and interpreting results.

**Funding:** This research received no external funding.

**Acknowledgments:** Landsat 8-OLI imagery used in this work were provided by the USGS (United States Geological Survey) through the Earth Explorer portal (https://earthexplorer.usgs.gov/). Sentinel-2 data are made free available online by the SENTINEL Hub (http://www.sentinel-hub.com/).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
