*3.1. Volcanic Radiative Power by Using MODIS-MIROVA Data*

MIROVA is a fully automatic volcano-dedicated hot-spot detection system ([49]) based on the images acquired by the MODIS sensor mounted on TERRA and AQUA NASA's satellites.

The combination of spatial and spectral principles to the middle infrared images acquired by MODIS makes possible the identification of the hot-spot contaminated pixel(s) (1 km resolution) and, consequently, the quantification of high temperature (>600 K) thermal anomalies, in terms of Volcanic Radiative Power (VRP, in Watt) (see [50]). The VRP (in Watt) is calculated through the MIR-method [50]:

$$\text{VRP}\_{\text{pIX}} = A\_{\text{pIX}} \times 18.7 \times (\text{L}\_{\text{4alert}} - \text{L}\_{\text{4bk}}) \tag{1}$$

where APix, L4alert and L4bk are the pixel area (106 m2 for MODIS), and the middle infrared radiance (W m−<sup>2</sup> sr−<sup>1</sup> μm<sup>−</sup>1) recorded by the alerted (L4alert) and background (L4bk) pixels, respectively. According to [50] for hot targets that have an integrated temperature comprised between 600 and 1500 K, the use of the constant of 18.7 sr·μm provides reliable estimates of radiant power with an uncertainty of ±30%. This makes VRP particularly appropriate for calculating the heat flux sourced by only the active portions of lava flows, being almost insensitive to cooling lava surfaces having temperature lower than ∼300 ◦C [49].

Overall, during the September 2018–January 2019 period, MIROVA elaborated more than 657 images (four overpasses per day) of which 184 (ca. 28%) detected a thermal anomaly. The VRP values range from less than 1 MW to a maximum of 2295 MW, this last being recorded on December 24 at 21:15 UTC.

#### *3.2. Time Average Discharge Rate and Erupted Volume via MODIS Data*

Satellite thermal data represent a useful parameter to estimate of Time Averaged lava Discharge Rate (TADR) and erupted volumes during effusive eruptions [51]. This thermal approach relies on the observed relationships that characterize the effusion rates, the active flow area, and the thermal flux, as documented by [51–54].

Here, we adopted an empirical approach, proposed by [55], that directly relates the volcanic radiant power (VRP), measured via MODIS, to the Time Average Discharge Rate (TADR) through a unique, best-fit parameter, called radiant density (crad; in J m<sup>−</sup>3);

$$\text{TADR} = \frac{\text{VRP}}{\text{crad}} \tag{2}$$

This empirical parameter represents the bulk efficiency of the lava body to radiate heat from its active surface, and is mainly controlled by the rheological, insulation, and topographic conditions at the time of emplacement [49,55]. Consequently, by integrating VRP measurements in time to obtain the total Radiant Energy (VRE) we are able, to estimate the erupted volume by assuming the appropriate crad, for the observed lava flow.

Based on previous calibration of the radiant density typical of Etnean lava flows [49,55], here we used a crad between 2 and 3.6 <sup>×</sup> <sup>10</sup><sup>8</sup> J m−<sup>3</sup> that makes it possible to constrain the erupted lava volumes between a minimum and maximum value.

It should be noted that the above method relies on the assumption that all the heat recorded by space is produced by the outpoured lava from the surface and the emplacing of lava flow or lava body. However, during open-vent activity, no net lava can be outpoured from the volcano, despite a thermal anomaly is still detected. In these cases the thermal approach can be used to infer the rate at which magma reaches the uppermost levels of the conduit (i.e., the bottom of a crater), before being cycled back in the convective magma column [56–59]. In this view, a variation in magnitude and persistence of heat radiation from summit craters could indicate a different level of magma column [60].

#### *3.3. SENTINEL 2 Images*

The Multispectral Instrument (MSI) is carried on-board SENTINEL-2A/2B platforms, two ESA satellites launched on polar, sun-syncronous orbit on June 2015 and March 2017, respectively. MSI provides multispectral data in 13 bands from VNIR to SWIR spectral region, useful for volcano monitoring applications. Indeed, the 20 m/pixel high-spatial resolution in the SWIR bands makes it

possible to detail morphometric thermal features of the heat volcanic sources and distinguish active volcanic craters, or sectors, during ongoing eruptions (with a decameter spatial detail; [14,42]).

The SENTINEL-2 data are made available by the ESA Copernicus Service Data Hub (https: //cophub.copernicus.eu/dhus/#/home), were downloaded through the cloud storage service of Amazon Web Service S3 (AWS-S3, https://registry.opendata.aws/SENTINEL-2/). Over the Etna geographical area, SENTINEL-2 satellites have a revisit time of about 2–3 days considering the two inspecting platforms 2A and 2B, enabling us to acquire about 10 high-resolution images per month, on average.

Following Massimetti et al. [61] we used the TOA (Top of the Atmosphere) reflectance data with a band combination 12–11–8a in the SWIR spectral region (R: 2190 nm, ρ12; G: 1610 nm, ρ11; B: 865 nm, ρ8a).

To analyze the thermal signature and to enlighten the presence of hot-spots, we applied a simple approach based on spectral ratios analysis of 12–11–8a bands, detecting and locating the number of "hot" pixels for each image in which a thermal volcanic emission is ongoing [62]. The algorithm to detect hot-spot contaminated pixel applied here was already successfully tested on several worldwide volcanic targets in comparison with MODIS-MIROVA heat flux [61]. For each pixel detected as "hot", we calculated a Thermal Index (T.I.) as the sum of the reflectances in the three SWIR bands analyzed (e.g., T.I. = ρ <sup>12</sup>+ ρ <sup>11</sup>+ ρ 8A) which is here considered as a proxy of the heat source temperature. We use the Thermal Index map to create stacked thermal profiles over the summit crater area of the Etna volcano (Figure 4), enabling us to distinguish the different craters in terms of persistence and to quantify the thermal contribute of each defined sector (see section Results).

To summarize, the SENTINEL-2 data are used to:


#### *3.4. Infrasound Arrays and Seismic Tremor*

Infrasound monitoring at the Etna volcano is performed with two 4 elements small aperture infrasound arrays [22,63,64], named ETN and MVT and installed, respectively, on September 2007 and in 2015. The ETN array was installed at 2100 m a.s.l., at a distance of 5500 m from the summit craters on the southern rim of Valle del Bove, while the MVT was installed on the southern flank of volcano at 1800 m a.s.l. and at a distance of about 6500 m from the craters. Both arrays are equipped with a FIBRA (www.item-geophysics.it) with four elements deployed following a triangular geometry and with an aperture (maximum distance between two array elements) of ~250 m and ~150 m, respectively ([22]; see Figure 1). Each array element is equipped with a differential pressure transducer with a sensitivity of 25 mV/Pa, maximum pressure range of ± 100 Pa, and flat frequency response between 0.01 and 100 Hz. Infrasound pressure data are converted to digital at each array element, transmitted with fiber optic cable to the central array element, where data is collected and GPS time stamped. The use of fiber optic made it possible to increase the signal-to-noise ratio and reducing the damages related to lightning and electrical discharges, with the ETN array being operational continuously since 2007. At ETN the array is co-located with a Guralp CMG-6T broadband seismometer, with eigen-period of 10 s and sensitivity of 2000 V/m s<sup>−</sup>1.

Infrasound array data are processed in the time domain by applying a grid search procedure for a source within the crater area. The analysis is applied on 5-s-long time window. The sound speed velocity is fixed within the 330–360 m s−<sup>1</sup> apparent velocity range, which accounts for a temperature effect on the acoustic propagation velocity and for the source-receiver elevation differences (1300 ± 100 m for the ETN array and 1600 ± 100 m for the MVT array). These two constraints allow the automatic removal of all the sources different from the volcano. The algorithm ([65,66] for details) searches for coherent signals recorded across the array, and we define a detection when the semblance exceeds 0.5. For each detection, the corresponding back-azimuth and the mean acoustic pressure are thus calculated according to the delay times observed among the different channels. Considering the aperture of the two arrays (<250 m) and the typical frequency of 1 Hz, the expected azimuth resolution is ~2◦ [67], which corresponds to ~190–230 m horizontal resolution at a source-to-receiver distance of 5500–6500 m. This processing is applied to the infrasonic data recorded at both the ETN and MVT arrays.

Here, we analyze the MVT array data because its capability, due to the location (Figure 1), to discriminate between Central Craters (BN/VOR and NE) and SE (NSE and SE) explosive activity.
