Mt. Etna Paroxysms of February–April 2021 Monitored and Quantified through a Multi-Platform Satellite Observing System

: On 16 February 2021, an eruptive paroxysm took place at Mt. Etna (Sicily, Italy), after continuous Strombolian activity recorded at summit craters, which intensiﬁed in December 2020. This was the ﬁrst of 17 short, but violent, eruptive events occurring during February–April 2021, mostly at a time interval of about 2–3 days between each other. The paroxysms produced lava fountains (up to 1000 m high), huge tephra columns (up to 10–11 km above sea level), lava and pyroclastic ﬂows, expanding 2–4 km towards East and South. The last event, which was characterised by about 3 days of almost continuous eruptive activity (30 March–1 April), generated the most lasting lava fountain (8–9 h). During some paroxysms, volcanic ash led to the temporary closure of the Vincenzo Bellini Catania International Airport. Heavy ash falls then affected the areas surrounding the volcano, in some cases reaching zones located hundreds of kilometres away from the eruptive vent. In this study, we investigate the Mt. Etna paroxysms mentioned above through a multi-platform satellite system. Results retrieved from Advanced Very High Resolution Radiometer (AVHRR), Moderate Resolution Imaging Spectroradiometer (MODIS), and Spinning Enhanced Visible and Infrared Imager (SEVIRI), starting from outputs of the Robust Satellite Techniques for Volcanoes (RST VOLC ), indicate that the 17th paroxysm (31 March–1 April) was the most powerful, with values of radiative power estimated around 14 GW. Moreover, by the analysis of SEVIRI data, we found that the 5th and 17th paroxysms were the most energetic. The Multispectral Instrument (MSI) and the Operational Land Imager (OLI), providing shortwave infrared (SWIR) data at 20/30 m spatial resolution, enabled an accurate localisation of active vents and the mapping of the areas inundated by lava ﬂows. In addition, according to the Normalized Hotspot Indices (NHI) tool, the 1st and 3rd paroxysm (18 and 28 February) generated the largest thermal anomaly at Mt. Etna after June 2013, when Landsat-8 OLI data became available. Despite the impact of clouds/plumes, pixel saturation, and other factors (e.g., satellite viewing geometry) on thermal anomaly identiﬁcation, the used multi-sensor approach allowed us to retrieve quantitative information about the 17 paroxysms occurring at Mt. Etna. This approach could support scientists in better interpreting changes in thermal activity, which could lead to future and more dangerous eruptions. was the most powerful, and probably radiated the largest amount of energy. These results conﬁrm the importance of using a multi-sensor approach in investigating and quantifying thermal volcanic activity. This study demonstrates the efﬁciency of the used multi-platform system in monitoring and quantifying Mt. Etna activity, despite the impact of a number of factors (e.g., clouds/plumes, pixel saturation, satellite viewing geometry, and atmospheric effects) on the identiﬁcation and mapping of volcanic thermal anomalies from space. This system may provide reliable information also about subtle phases of thermal unrest (e.g., those preceding intense eruptions), contributing to the volcanic risk mitigation also in other areas.


Introduction
The Advanced Very High Resolution Radiometer (AVHRR), the Moderate Resolution Imaging Spectroradiometer (MODIS), and the Spinning Enhanced Visible and Infrared
In this work, we assess the advantages and limitations of a multi-platform observing system, integrating data at different spatial, temporal, and spectral resolutions from low-(LEO) and geo-(GEO) orbiting satellite sensors, in monitoring and characterising the Mt. Etna (Sicily, Italy) paroxysms, i.e., intense short-lived eruptions occurring with Strombolian activity, lava fountaining, lava overflows, and tephra columns.
Mt. Etna ( Figure 1) is the highest and most active volcano in Europe. Its eruptive activity frequently occurs from summit craters, although some intense flank eruptions, which are usually more dangerous for both infrastructures and population in the presence of eruptive fissures opening at low altitudes, occurred in recent years [21][22][23]. In the last decade, several paroxysms took place at the summit craters [24]. During January 2011−December 2013, 46 eruptive episodes were recorded [25]. The paroxysms of January−July 2011 yielded the growth of the New South-East Crater (NSEC) at the eastern base of the 'old' South-East Crater (SEC, born in 1971 [26]), continuing in 2013 [27,28]. In 2015, other powerful paroxysms occurred at the Voragine crater (VOR), generating a high eruptive column (up to 15 km above sea level) [29]. In recent years, a brief lateral eruption took place on 24-27 December 2018. A huge movement of the eastern flank of the volcano and a seismic crisis (up to Mw 4.9) caused damages in a vast urbanised area of the eastern flank of Etna [12]. After this brief but intense flank eruption, the South-East Crater complex (SEC) awoke on 30 May 2019 and then reactivated on 18 and 27 July 2019, producing Strombolian activity. On 12 September 2019, the VOR crater came into activity with Strombolian explosions and intra-crater lava overflow, and then again, the SEC re- In the last decade, several paroxysms took place at the summit craters [24]. During January 2011-December 2013, 46 eruptive episodes were recorded [25]. The paroxysms of January-July 2011 yielded the growth of the New South-East Crater (NSEC) at the eastern base of the 'old' South-East Crater (SEC, born in 1971 [26]), continuing in 2013 [27,28]. In 2015, other powerful paroxysms occurred at the Voragine crater (VOR), generating a high eruptive column (up to 15 km above sea level) [29]. In recent years, a brief lateral eruption took place on 24-27 December 2018. A huge movement of the eastern flank of the volcano and a seismic crisis (up to Mw 4.9) caused damages in a vast urbanised area of the eastern flank of Etna [12]. After this brief but intense flank eruption, the South-East Crater complex (SEC) awoke on 30 May 2019 and then reactivated on 18 and 27 July 2019, producing Strombolian activity. On 12 September 2019, the VOR crater came into activity with Strombolian explosions and intra-crater lava overflow, and then again, the SEC reactivated in the spring of 2020, with frequent Strombolian eruptions culminating in four major episodes (December 2020), and finally in 17 eruptive paroxysms recorded during February-April 2021, which mostly occurred at intervals of around 2-2.5 days between Remote Sens. 2021, 13,3074 3 of 20 each other [30]. During the night of [18][19] May, a series of new eruptive paroxysms started and are still ongoing at the time of writing this article.
Here, we investigate the 17 Mt. Etna paroxysms of February-April 2021 through the outputs of the Robust Satellite Techniques for Volcanoes (RST VOLC ) [31], the RASTer (RST on ASTER [32]) system, and the Normalized Hotspots Indices (NHI) tool [33]. To our knowledge, this is one of the first studies focusing on those events, which were highly spectacular and particularly intense, receiving special attention from international, national, and social media.

Mt. Etna Paroxysms of February-April 2021
On 16 February 2021 (Figure 2), after an almost continuous Strombolian activity recorded in the previous months at Mt. Etna summit craters, which culminated with the major Strombolian episodes of 13-14, 21 and 22 December 2020, a sequence of powerful paroxysms took place at SEC. Table A1 (see Appendix A) details the 17 paroxysms, which, until 1 April, followed one another with intervals ranging from 32 h to over 6 days, providing information about eruption features.
Remote Sens. 2021, 13, x FOR PEER REVIEW 3 of 20 activated in the spring of 2020, with frequent Strombolian eruptions culminating in four major episodes (December 2020), and finally in 17 eruptive paroxysms recorded during February-April 2021, which mostly occurred at intervals of around 2-2.5 days between each other [30]. During the night of 18-19 May, a series of new eruptive paroxysms  started and are still ongoing at the time of writing this article. Here, we investigate the 17 Mt. Etna paroxysms of February-April 2021 through the outputs of the Robust Satellite Techniques for Volcanoes (RSTVOLC) [31], the RASTer (RST on ASTER [32]) system, and the Normalized Hotspots Indices (NHI) tool [33]. To our knowledge, this is one of the first studies focusing on those events, which were highly spectacular and particularly intense, receiving special attention from international, national, and social media.

Mt. Etna paroxysms of February-April 2021
On 16 February 2021 (Figure 2), after an almost continuous Strombolian activity recorded in the previous months at Mt. Etna summit craters, which culminated with the major Strombolian episodes of 13-14, 21 and 22 December 2020, a sequence of powerful paroxysms took place at SEC. Table A1 (see Appendix A) details the 17 paroxysms, which, until 1 April, followed one another with intervals ranging from 32 h to over 6 days, providing information about eruption features. The 17th paroxysm was the longest of the time series (with lava fountains of 540 min duration), the eruptive column was usually high, up to 11 km above sea level (a.s.l.) (e.g., the 6th paroxysm), and lava flows moved towards Valle del Bove and the southern flank (see Table A1).
During the eruptive events, primitive magma, which is compositionally similar to the source magma that forms in the Earth's mantle at a depth of over 40 km [34], very rich in gas, was expelled (e.g., [35]). The coarser and heavier material fell by gravity around the eruptive vent and therefore on the flanks of the SEC, which through this continuous contribution has systematically changed shape and size in each of the 17 eruptions. The lighter portion of the fragmented magma, however, continued to rise in the atmosphere for kilometres, until it reached the 'floating' altitude, around 8-10 km above sea level, where ashes and lapilli were taken over by the winds of the high troposphere and then transported to areas very far from the volcano, up to many hundreds of kilometres. Therefore, most of the Etna urban centres were repeatedly covered by a heavy mantle of ashes and lapilli. A phenomenon that is only apparently harmless. Indeed, the fallout of the coarser volcanic material, with bombs larger than 10-15 cm falling over 7 km away The 17th paroxysm was the longest of the time series (with lava fountains of 540 min duration), the eruptive column was usually high, up to 11 km above sea level (a.s.l.) (e.g., the 6th paroxysm), and lava flows moved towards Valle del Bove and the southern flank (see Table A1).
During the eruptive events, primitive magma, which is compositionally similar to the source magma that forms in the Earth's mantle at a depth of over 40 km [34], very rich in gas, was expelled (e.g., [35]). The coarser and heavier material fell by gravity around the eruptive vent and therefore on the flanks of the SEC, which through this continuous contribution has systematically changed shape and size in each of the 17 eruptions. The lighter portion of the fragmented magma, however, continued to rise in the atmosphere for kilometres, until it reached the 'floating' altitude, around 8-10 km above sea level, where ashes and lapilli were taken over by the winds of the high troposphere and then transported to areas very far from the volcano, up to many hundreds of kilometres. Therefore, most of the Etna urban centres were repeatedly covered by a heavy mantle of ashes and lapilli. A phenomenon that is only apparently harmless. Indeed, the fallout of the coarser volcanic material, with bombs larger than 10-15 cm falling over 7 km away from the eruptive vent, can injure people and animals, damage and burn vegetation, break crystals, and damage car bodyworks. Volcanic ash makes roads slippery and dangerous for road traffic. The passage of cars on the roads covered with lapilli further fragments the volcanic material, making it microscopic in size (up to PM10 size), easily inhaled, and therefore dangerous for human health. Finally, when the material falls on Catania, Sigonella, and Reggio Calabria airports, they must temporarily close, to ensure the safety of flights and clean the runways. Multiplying these problems 17 times in 44 days, we can understand the hardship faced by the population as well as the difficulties that the authorities have had to face.

Data
The AVHRR sensor, aboard NOAA (National Oceanic and Atmospheric Administration) satellites, provides data in five spectral channels, from visible to thermal infrared. The AVHRR3, flying on the MetOp (Meteorological Operational)-A, B, and C satellites, scans the Earth's surface in six spectral bands, including the channel 3A (1.6 µm) operating in daytime [36]. MODIS, on board the Aqua and Terra satellites, acquires VIS and TIR data four times per day in 36 spectral bands (e.g., [3]). SEVIRI, aboard Meteosat Second Generation (MSG) geostationary satellites, provides data in four VNIR (0.4-1.6 µm) and eight IR channels (3.9-13.4 µm), with a temporal sampling up to 5 min in the Rapid Scanner mode [37]. ASTER, aboard the Terra platform, operates in 14 spectral bands from VNIR to TIR (15 m to 90 m spatial resolution), with a repeat cycle of 16 days. OLI, aboard the Landsat-8 satellite (16-day revisit time), acquires VNIR to SWIR data in nine spectral channels [38]. The MSI instrument, aboard the Sentinel-2 constellation (5 day revisit time at the equator in cloud-free conditions), acquires data in 13 spectral bands from VNIR to SWIR, with a spatial resolution ranging from 10 m to 60 m [39].
Satellite data from those sensors ( Table 1) were used to monitor and quantify the Mt. Etna paroxysms, using the hotspot detection methods described in detail in the next section. The RST VOLC algorithm is an optimised configuration of the Robust Satellite Techniques (RST) for volcanological applications (see [31,40] for details).
RST VOLC uses two indices in combination, based on the Absolutely Local Index of Change of the Environment (ALICE [41]), to detect thermal anomalies [41]. The ⊗ MIR (x,y,t) index, which is computed using the AVHRR channel 3B, the MODIS channel 21/22, and the SEVIRI channel 4 (see Table 1), is highly sensitive to hot magmatic features. The ⊗ MIR-TIR (x,y,t) index, by exploiting also TIR data (i.e., AVHRR channel 4, MODIS channel 31, and SEVIRI channel 9), is used to increase the confidence level of detection [40]. Since RST VOLC is a tuneable technique, different cutting levels of those indices (e.g., ⊗ MIR (x,y,t) > 3 AND ⊗ MIR-TIR (x,y,t) > 3) are used to identify volcanic hotspots [31].
Those indices are calculated by processing pluriannual time series of homogeneous (i.e., same spectral channel/s, same month, and overpass times) cloud-free satellite records. The One-channel Cloudy-radiance-detection Approach (OCA [42]), and a standard cloudmasking procedure, from the European Organization for the Exploitation of Meteorological Satellites (EUMETSAT), are used to filter out cloudy pixels from SEVIRI data [43].
The RST detection scheme applied to ASTER TIR data, under Google Earth Engine (GEE) platform, is widely described in [32].

NHI Algorithm
The NHI multichannel algorithm uses two normalised indices to identify and map high-temperature features through satellite data at mid-high spatial resolution [44]. Those indices analyse the TOA (Top of the Atmosphere) radiance (W m -2 sr -1 µm -1 ) measured in the SWIR (at 1.6 µm and 2.2 µm) and NIR (Near Infrared) (0.8 µm) bands of MSI and OLI. The algorithm, whose wide description can be found in [44], flags as 'hot' the pixels with positive values of one or both the above-mentioned indices. NHI was tested with success over different volcanic areas using data also from previous sensors such as TM (Thematic Mapper), ETM + (Enhanced Thematic Mapper, Plus), and ASTER [45], whose SWIR observations were available until April 2008 due to the failure of the SWIR instrument [46].
The NHI tool (https://sites.google.com/view/nhi-tool (accessed on 23 June 2021)), by exploiting the high computational resources and the extended data archive of the GEE platform, enables the analysis of volcanic thermal anomalies at the global scale. The tool includes some additional spectral tests to increase the accuracy in mapping hot targets minimising false detections. The latter are mainly associated with the multispectral misregistration of Sentinel-2 MSI imagery, affecting also other recent hotspot detection methods (e.g., [47]). A detailed description of the NHI tool performance can be found in [33].

AVHRR Estimations
To estimate the radiative power from infrared AVHRR data [48], starting from the outputs of the dual band three-component method, we used the following formulation: where Q rad is the radiative power [W], n is the number of detected hot pixels, A i is the lava flow area for the pixel i having a crusted and molten portion p c i and p h i with temperatures T c i and T h i , and σ is the Stefan-Boltzmann constant (5.67 × 10 -8 Wm -2 K -4 ) [48]. In this study, we assumed ε = 0.98 (e.g., [49,50]) as representative of the Etnean lava flows to calculate the radiative power (mean value from lower and upper limit; T c = 100 • C; Remote Sens. 2021, 13, 3074 6 of 20 500 • C), without correcting data for atmospheric effects to cross-compare results with the MODIS ones.

SEVIRI and MODIS Estimations
To retrieve the radiative power from SEVIRI and MODIS data, starting from thermal anomalies flagged by RST VOLC , we used the formulation proposed in previous studies to characterise thermal anomalies with temperatures ranging in between 600 and 1500 K [51,52], and largely used in literature (e.g., [53,54]): In Equation (2), A pixel is the pixel area (km 2 ) at relative scan angle (A pixel from MODIS data was calculated as ∆S × ∆T, where ∆S and ∆T are the along-scan and along-track pixel size [52,55], τ MIR is the atmospheric transmittance in the MIR band, σ is the Stefan-Boltzmann constant, ε f is the hotspot emissivity over all wavelengths, ε f,MIR is the hotspot emissivity in the MIR band, L MIR (W m -2 sr -1 µm -1 ) is the MIR spectral radiance from the hot target, L b,MIR is the background radiance, and a is a sensor-dependent constant ( Table 2). Table 2. Sensor-dependent constants [51][52][53].
In this work, L b,MIR from MODIS data refers to the temporal mean of the spatial average calculated over a 7 × 7 pixel box located at the Mt. Etna southern area (i.e., at a distance of about 10 km). The OCA method was used to discard cloudy pixels from this computation. To retrieve L b,MIR from SEVIRI data, we first calculated the temporal mean of the MIR radiance, i.e., µ LMIR (x, y, t), for each pixel of the scene at the observation time t (96 time slots) by analysing a long-time series of L MIR data [41]. For each hotspot pixel detected by RST VOLC , the background radiance represents the spatial average of µ LMIR (x, y, t) retrieved over a 7 × 7 pixel window (about 91 × 91 km 2 ), around the pixel (x,y). To compute the spatial average, we excluded both sea and no-data pixels.

Thermal Anomaly Identification and Mapping
In Figure 3, we show some thermal anomaly maps, referring to the period 16 February-10 March 2021, generated through the NHI tool and the RASTer system. Maps display the hotspot pixels in different colours. Regarding the NHI detections, red and yellow pixels indicate the areas where volcanic thermal emissions were strong or less intense, respectively. Remote Sens. 2021, 13, x FOR PEER REVIEW 7 of Figure 3. Thermal anomaly maps generated using the NHI tool (red/yellow pixels) and the RASTer (orange pixels) system. Maps show the thermal anomaly at the summit craters, and the detected lava flows associated with paroxysms at SEC.
After a gradual increase in Strombolian activity and volcanic tremor recorded by th Istituto Nazionale di Geofisica e Vulcanologia (INGV) surveillance network, a lava ove flow started on 16 February at around 17:05 UTC. The lava overflow caused the parti collapse of the cone at the eastern eruptive vents of SEC, generating a pyroclastic flo along the western wall of the Valle del Bove (see Table A1). Soon after, an intense lav fountaining (500-600 m high, 45 min lasting) occurred, generating a tall (≈10 km abov sea level) eruptive column [30]. A few hours before the start of the paroxysm, a therm anomaly, which was spatially more extended than previous days, affected the Mt. Etn crater area, as indicated by hot pixels flagged by NHI over VOR-BN and SEC (see th S2-MSI map of 16 February). The RASTer map, from ASTER TIR data of 17 Februar shows that the lava flow emitted during the first paroxysmal event affected the Valle d Bove, moving towards the E-SE direction (see orange pixels).
The second paroxysm took place during the night of 17-18 February, i.e., 33 h afte the previous eruptive episode. This event was similar to the previous one and generate a lava flow travelling towards the Valle del Bove, in the NE and SE, and on the hig southern flank of Etna, in the SW direction [57]. The S2-MSI and L8-OLI maps of 1 February provided the same information about the shape and spatial extent of lava flow revealing the occurrence of a thermal activity also at NEC. Figure 3 shows that the fir flow moved in the SW direction reaching an altitude of about 2840 m, while the secon flow extended up to about 1680 m in the Valle del Bove. The maps of 21 and 23 Februar display the lava flows associated with the 4th and 5th paroxysm, which mostly affecte the same areas. The fourth paroxysm started in the late afternoon of 20 February and wa preceded by a weak Strombolian activity. At around 21:30 UTC, a small lava overflo occurred; explosive activity then increased, and a lava fountain activity began at 22:0 UTC, becoming more intense in the following hours. This activity, generating an eruptiv column extending up to about 10 km a.s.l., decreased starting from 01:10 UTC [30]. Th fifth paroxysm took place during 22-23 February with lava fountains lasting about 5 min, and the emission of a high eruptive column. It had a time duration of about 14 including also the Strombolian phase preceding the fountain (Table 1A).
The L8-OLI and Terra-ASTER maps of 25 and 26 February display the therm anomaly associated with the sixth paroxysm, which took place on 24 February, lastin  Table A1). Soon after, an intense lava fountaining (500-600 m high, 45 min lasting) occurred, generating a tall (≈10 km above sea level) eruptive column [30]. A few hours before the start of the paroxysm, a thermal anomaly, which was spatially more extended than previous days, affected the Mt. Etna crater area, as indicated by hot pixels flagged by NHI over VOR-BN and SEC (see the S2-MSI map of 16 February). The RASTer map, from ASTER TIR data of 17 February, shows that the lava flow emitted during the first paroxysmal event affected the Valle del Bove, moving towards the E-SE direction (see orange pixels).
The second paroxysm took place during the night of 17-18 February, i.e., 33 h after the previous eruptive episode. This event was similar to the previous one and generated a lava flow travelling towards the Valle del Bove, in the NE and SE, and on the high southern flank of Etna, in the SW direction [57]. The S2-MSI and L8-OLI maps of 18 February provided the same information about the shape and spatial extent of lava flows, revealing the occurrence of a thermal activity also at NEC. Figure 3 shows that the first flow moved in the SW direction reaching an altitude of about 2840 m, while the second flow extended up to about 1680 m in the Valle del Bove. The maps of 21 and 23 February display the lava flows associated with the 4th and 5th paroxysm, which mostly affected the same areas. The fourth paroxysm started in the late afternoon of 20 February and was preceded by a weak Strombolian activity. At around 21:30 UTC, a small lava overflow occurred; explosive activity then increased, and a lava fountain activity began at 22:00 UTC, becoming more intense in the following hours. This activity, generating an eruptive column extending up to about 10 km a.s.l., decreased starting from 01:10 UTC [30]. The fifth paroxysm took place during 22-23 February with lava fountains lasting about 50 min, and the emission of a high eruptive column. It had a time duration of about 14 h, including also the Strombolian phase preceding the fountain (Table 1A).
The L8-OLI and Terra-ASTER maps of 25 and 26 February display the thermal anomaly associated with the sixth paroxysm, which took place on 24 February, lasting about six hours. Lava fountain activity fed a new lava flow expanding towards SW [58], which was well detected by satellite, as indicated by the NHI and RASTer maps.
Lava flow associated with the 7th paroxysm was identified almost in concomitance with the end of the eruptive event (i.e., from the S2-MSI scene of 28 February at 09:44 UTC). Consequently, the lava flow, extending up to an altitude of about 2000 m in the Valle del Bove, appeared particularly intense (see the high number of red pixels). The seventh paroxysm occurred after a three-and-a-half-day break and was preceded by a weak Strombolian, having started on 28 February at around 08:10 UTC, rapidly evolving to lava fountains, lasting 54 min. The paroxysm also produced a lava flow moving towards East, and a tall eruptive column (11 km a.s.l.) causing the ash fall-out in the East direction.
The thermal anomaly map from the S2-MSI scene of 3 March at 09:51 UTC displays the lava flow associated with the 8th paroxysm, which had occurred the day before with lava fountains, lasting about two hours, and an eruptive column dispersing in the South direction. Thermal anomaly associated with the paroxysm was partially detected by satellite (as indicated by the manual inspection of satellite imagery), due to possible lava cooling effects. In general, those effects such as clouds/plumes had a significant impact on thermal anomaly identification, as for the maps of 6 and 10 March. On 4 March, the 9 th paroxysm occurred. It was preceded by a Strombolian activity increasing 22 min after midnight, reaching the climax phase about two hours later. The lava fountain lasted just over two hours. The 10th paroxysm started on 7 March between 00:00 and 01:00 UTC and lasted about 10 h, emitting a small lava flow from a vent opening at the base of SEC. Activity then evolved to lava fountains at 06:00 UTC lasting 60 min, generating a high eruptive column, dispersing in the East direction, and causing the ash fall-out over the neighbour villages [59]. The map of 6 March, referring to the period in between the 9th and 10th paroxysm, was strongly affected by cloud coverage, while the map of 10 March shows the lava flow emitted the day before, during the 11th paroxysm, from a vent located at the S base of the SEC [60], which was masked by clouds mainly in the distal part.
The map from the S2-MSI scene of 13 March shows a small thermal anomaly detected after the end of the 12th paroxysm. This event, which was preceded by a weak Strombolian activity occurring on 12 March since 03:20 UTC, generated lava fountaining lasting 175 min, and an eruptive column extending up about 9-10 km a.s.l. (Table A1) On the other hand, the analysis of the false colour composite imagery revealed that this feature was spatially more extended than that indicated by NHI. The latter did not detect the lava flow in the distal part (where it was less radiant in the SWIR bands), due to possible effects of lava cooling.
Finally, the map of 30 March shows that a low-to-moderate thermal activity was in progress at the summit craters, some hours before the occurrence of the 17th paroxysm. This event was preceded by an intense Strombolian activity occurring since 14:30 UTC [64]. A small lava emission was noticed through INGV thermal surveillance cameras, emitted from the SEC's eastern base at 13:00 UTC on 31 March. Lava fountains started a few hours later and lasted 540 min. The eruptive column reached 9-10 km a.s.l., dispersing in the SSW direction. Lava flows expanded towards E, SE, and SSE (in Valle del Bove) and towards S and SW (on the southern flank) down to 1820 m a.s.l.
Hence, on the one hand, S2-MSI, L8-OLI, and Terra-ASTER data provided detailed information about active vents and emitted lava flows, particularly when they were less affected by cloud/volcanic plumes. On the other hand, they did not enable the continuous monitoring of Mt. Etna activity, despite the advantages of data integration, which guarantees a temporal sampling of 2-3 days on average over the Mt. Etna area (see Figure 4). Therefore, thermal activity associated with the paroxysms of February-April 2021 required the use of high temporal resolution satellite data to be monitored and quantified in a more efficient way from space, as shown and discussed in detail in the next sections.
Remote Sens. 2021, 13, x FOR PEER REVIEW 9 of 20 affected by cloud/volcanic plumes. On the other hand, they did not enable the continuous monitoring of Mt. Etna activity, despite the advantages of data integration, which guarantees a temporal sampling of 2-3 days on average over the Mt. Etna area (see Figure 4). Therefore, thermal activity associated with the paroxysms of February−April 2021 required the use of high temporal resolution satellite data to be monitored and quantified in a more efficient way from space, as shown and discussed in detail in the next sections.

AVHRR and MODIS Observations
The RSTVOLC system performs the near-real-time monitoring of Italian volcanoes through infrared AVHRR and MODIS data, directly acquired and processed at the Institute of Methodologies for Environmental Analysis (IMAA) of the National Research Council (CNR). Data integration currently guarantees up to nine satellite observations per day (AVHRR-3 data are not analysed in the daytime because of channel 3A). The system generates thermal anomaly products (i.e., maps and ASCII files), available within a few minutes after the sensing time [31], which may be delivered to the users upon request. Figure 5 displays the temporal trend of radiative power retrieved over the period 16 February−1 April 2021, integrating results from the above-mentioned sensors, starting from the RSTVOLC detections.

AVHRR and MODIS Observations
The RST VOLC system performs the near-real-time monitoring of Italian volcanoes through infrared AVHRR and MODIS data, directly acquired and processed at the Institute of Methodologies for Environmental Analysis (IMAA) of the National Research Council (CNR). Data integration currently guarantees up to nine satellite observations per day (AVHRR-3 data are not analysed in the daytime because of channel 3A). The system generates thermal anomaly products (i.e., maps and ASCII files), available within a few minutes after the sensing time [31], which may be delivered to the users upon request. Figure 5 displays the temporal trend of radiative power retrieved over the period 16 February-1 April 2021, integrating results from the above-mentioned sensors, starting from the RST VOLC detections.
The figure shows that AVHRR provided relevant information especially about the 1st, 4th, and 9th event, due to satellite overpass times fitting with the time of paroxysms. On the other hand, the other eruptive episodes (e.g., the 3rd, 8th, 12th, 15th, and 17th paroxysms), mostly occurring in the daytime, were better identified using MODIS. The integration of AVHRR and MODIS observations indicates that some eruptive events were particularly intense (see the 1st, 3rd, and 17th paroxysms). Among them, the 17th paroxysm appeared as the most powerful, with a value of radiative power reaching about 14 GW (from the daytime MODIS scene of 1 April at 09:33 UTC; see the thermal anomaly in Figure 6b). Figure 5 then shows the progressive increase/decrease in thermal volcanic activity recorded before/after some eruptive events (e.g., see the 1st and 12th paroxysms), and the identification of a subtle thermal anomaly also between the 16th and 17th paroxysms. The figure shows that AVHRR provided relevant information especially about the 1st, 4th, and 9th event, due to satellite overpass times fitting with the time of paroxysms. On the other hand, the other eruptive episodes (e.g., the 3rd, 8th, 12th, 15th, and 17th paroxysms), mostly occurring in the daytime, were better identified using MODIS. The integration of AVHRR and MODIS observations indicates that some eruptive events were particularly intense (see the 1st, 3rd, and 17th paroxysms). Among them, the 17th paroxysm appeared as the most powerful, with a value of radiative power reaching about 14 GW (from the daytime MODIS scene of 1 April at 09:33 UTC; see the thermal anomaly in Figure 6b). Figure 5 then shows the progressive increase/decrease in thermal volcanic activity recorded before/after some eruptive events (e.g., see the 1st and 12th paroxysms), and the identification of a subtle thermal anomaly also between the 16th and 17th paroxysms.
These results corroborate the advantages of AVHRR and MODIS data integration in monitoring Mt. Etna activity (see also [17]). Nonetheless, Figure 5 shows that some eruptive events (e.g., the 2nd, 6th, 10th, and 11th paroxysms) were monitored in a less efficient way, because of their time duration. SEVIRI data processed through the RSTVOLC algorithm enabled the promptly identification of thermal anomalies and the continuous monitoring of volcanic activity (see next section). These results corroborate the advantages of AVHRR and MODIS data integration in monitoring Mt. Etna activity (see also [17]). Nonetheless, Figure 5 shows that some eruptive events (e.g., the 2nd, 6th, 10th, and 11th paroxysms) were monitored in a less efficient way, because of their time duration. SEVIRI data processed through the RST VOLC algorithm enabled the promptly identification of thermal anomalies and the continuous monitoring of volcanic activity (see next section).  Figure 7 displays the temporal trend of radiative power retrieved from infrared SEVIRI data (at 15 min time interval; acquired at the School of Engineering of University of Basilicata) of 16 February-2 April 2021. The plot also shows the time interval of each paroxysm (in pink), and the cloudiness (grey line) retrieved using the cloud-masking procedure described in the 'methods' section.  Figure 7 displays the temporal trend of radiative power retrieved from infrared SEVIRI data (at 15 min time interval; acquired at the School of Engineering of University of Basilicata) of 16 February-2 April 2021. The plot also shows the time interval of each paroxysm (in pink), and the cloudiness (grey line) retrieved using the cloud-masking procedure described in the 'methods' section.

SEVIRI Observations
(a) (b) Figure 6. (a) Thermal anomaly map automatically generated at IMAA-CNR from the night-time AVHRR overpass of 16 February at 20:51 UTC (i.e., after the end of the 1st paroxysm), red and yellow pixels, respectively, indicate the areas where thermal emissions were stronger or less intense; (b) The thermal anomaly (in red) detected by RSTVOLC over the Mt. Etna area on the daytime MODIS scene of 1 April at 09:33 UTC (channel 22 in the background) leading to the peak of radiative power in Figure 5. Figure 7 displays the temporal trend of radiative power retrieved from infrared SEVIRI data (at 15 min time interval; acquired at the School of Engineering of University of Basilicata) of 16 February-2 April 2021. The plot also shows the time interval of each paroxysm (in pink), and the cloudiness (grey line) retrieved using the cloud-masking procedure described in the 'methods' section. On that day, the 3rd paroxysm occurred (between 08:45 and 08:50 UTC), with a lava fountaining activity that was recorded after a lava overflow occurring since 07:55 UTC. This paroxysm, ending at 10:10 UTC, generated a dense ash plume dispersing towards SE direction, and causing the ash-fall over Zafferana Etnea and Acireale villages [30]. Based on RST VOLC detections, thermal anomaly of 19 February became more intense between 08:30 and 08:45 UTC, marking the start of the paroxysm. At 10:30 UTC, i.e., soon after the end of the eruptive event, volcanic radiative power was estimated around 3.8 GW. After the discontinuous identification of weak thermal activity in the morning of 20 February (between 06:00 and 08:30 UTC), a new phase of thermal unrest occurred at 21:00 UTC preceding the fourth paroxysm. Thermal anomaly, identified two days later, at 22:45 UTC, marked the start of the fifth paroxysm. The latter, based on information inferred from SEVIRI data, started at 23:30 UTC and reached the peak (≈2.8 GW) a few hours later (i.e., on 23 February at 00:45 UTC). On 24 February, RST VOLC detected a subtle thermal anomaly between 06:00 and 07:45 UTC and at 16:30 UTC, preceding the sixth paroxysm. The latter reached the peak at 20:45 UTC. While the identification of the thermal anomaly associated with the 8th paroxysm was strongly affected by clouds (the radiative power was estimated after the end of the paroxysm), during the period 28 February at 18:15 UTC-2 March at 16:00 UTC, a short quiescence period occurred, in agreement with field observations (see Table A1).

SEVIRI Observations
On the other hand, another phase of thermal unrest was identified on 2 March, starting from 16:15 UTC, i.e., before the start of the ninth paroxysm. The following two events were less intense in terms of VRP (with values below 3 GW), whereas the 12th paroxysm was among the most powerful. Clouds then strongly affected the crater area, impacting on the identification of thermal anomaly associated with the 13th, 14th, and 15th paroxysms (see high cloudiness values in Figure 7). The 13th paroxysmal event started on 14 March at 20:10 UTC, with the increase in Strombolian activity at SEC. The lava fountaining occurred a few minutes after midnight and ceased by 02:43 UTC [65]. The following event (the 14th of the time series) began on 17 March when the Strombolian activity recorded at 00:30 UTC increased in intensity, changing into lava fountaining at around 02:15 UTC. The paroxysm ended at 06:00 UTC [66,67]. Explosions at SEC and the lava flow in the Valle del Bove persisted until the evening of 18 March, when weather conditions improved. The 15th paroxysm occurred on 19 March with lava fountains recorded since 08:35 UTC, lasting about 120 min. This event followed the same pattern of previous events, with the increase in volcanic tremor, a Strombolian activity evolving to lava fountains, and the lava flow emission [66]. Figure 7 shows that a subtle thermal anomaly was then identified by RST VOLC before the 16th paroxysm. From the night of 25 March, and over more than 150 h, the system did not flag any hotspot over the target area. This quiescence period was the longest recorded between two consecutive paroxysms. Thermal anomaly detected from 31 March at 18:45 UTC marked the start of a new lava fountaining. The 17th paroxysm appeared, however, much less intense than that indicated by the analysis of infrared MODIS data (e.g., VRP ≈ 2.6 GW on 1 April at 08:45 UTC). The intensity of thermal emissions then progressively decreased, and a subtle hotspot affected the Mt. Etna area until 2 April at 07:30 UTC. This information agrees with field observations, indicating that in the late morning of 1 April, explosive activity at SEC was almost exhausted, unlike the effusive activity, which gradually ran out during the night of 1-2 April.

Discussion
We investigated the Mt. Etna paroxysms of February-April 2021 through a multiplatform observing system, retrieving information about location, size, shape, and relative intensity level of thermal anomalies detected by satellite. Figure 8 displays the plot of the total hotspot area retrieved from S2-MSI and L8-OLI scenes of June 2013-March 2021 using the NHI tool. The figure shows that, although other significant thermal anomalies (e.g., 26 October 2013, 18 March 2017) affected the target area before the eruptive events investigated in this work, those associated with the paroxysms of 16-18 February 2021 were the largest in size. It should be pointed out, however, that the low temporal sampling of the analysed satellite data, such as other factors (e.g., clouds/eruptive plumes), probably affected the results of Figure 8.
Despite those limitations, L8-OLI, S2-MSI, and Terra-ASTER data allowed us to investigate also changes in thermal volcanic activity occurring at Mt. Etna. This is indicated, for instance, by the small thermal anomaly identified soon after the lava effusion of 15 December 2020 [68], becoming spatially more extended after mid-January 2021. This detection reveals a general increase in thermal volcanic activity at the summit craters about one month before the start of the first paroxysmal event of 2021 (see Figure 8).
The latter, such as the following paroxysms, was better monitored and quantified using satellite data at high temporal resolution (see Section 5.2). Indeed, especially SEVIRI data enabled the prompt identification and the continuous monitoring of thermal volcanic activity (e.g., [5,7,69]), except when clouds/plumes partially or completely obscured the target area (e.g., the 8th paroxysm). Figure 7 shows that the percentage of cloudy pixels was generally high (i.e., above 50%). Because of clouds/plumes, some eruptive events (e.g., the 8th, 11th, and 16th paroxysms) were identified later than field observations, while other paroxysms (e.g., the 13th, 14th, and 15th paroxysms) were partially identified by satellite. Clouds affected less than 30% of the analysed sub-scenes (a circle area of 12 km in radius centred over the Mt. Etna craters) when SEVIRI data were fully exploited (e.g., the 3rd and 12th paroxysm).
scenes of June 2013−March 2021 using the NHI tool. The figure shows that, although other significant thermal anomalies (e.g., 26 October 2013, 18 March 2017) affected the target area before the eruptive events investigated in this work, those associated with the paroxysms of 16-18 February 2021 were the largest in size. It should be pointed out, however, that the low temporal sampling of the analysed satellite data, such as other factors (e.g., clouds/eruptive plumes), probably affected the results of Figure 8.
Despite those limitations, L8-OLI, S2-MSI, and Terra-ASTER data allowed us to investigate also changes in thermal volcanic activity occurring at Mt. Etna. This is indicated, for instance, by the small thermal anomaly identified soon after the lava effusion of 15 December 2020 [68], becoming spatially more extended after mid-January 2021. This detection reveals a general increase in thermal volcanic activity at the summit craters about one month before the start of the first paroxysmal event of 2021 (see Figure 8). The latter, such as the following paroxysms, was better monitored and quantified using satellite data at high temporal resolution (see Section 5.2). Indeed, especially SEVIRI data enabled the prompt identification and the continuous monitoring of thermal volcanic activity (e.g., [5,7,69]), except when clouds/plumes partially or completely obscured the target area (e.g., the 8th paroxysm). Figure 7 shows that the percentage of cloudy pixels was generally high (i.e., above 50%). Because of clouds/plumes, some eruptive events (e.g., the 8th, 11th, and 16th paroxysms) were identified later than field observations, while other paroxysms (e.g., the 13th, 14th, and 15th paroxysms) were partially identified by satellite. Clouds affected less than 30% of the analysed sub-scenes (a circle area of 12 km in radius centred over the Mt. Etna craters) when SEVIRI data were fully exploited (e.g., the 3rd and 12th paroxysm).
Regarding the radiative power, values retrieved from SEVIRI data appeared underestimated when compared to temporally coincident MODIS observations. In addition, even with filtering MODIS data for low values of the satellite zenith angle, to take the poor viewing conditions (e.g., [11,20,70,71]) into account, differences in terms of VRP recorded especially during the 17th paroxysm were particularly high (≈12.5 GW). Those differences were less significant when weak thermal activities occurred. The MIR channel saturation, having a lower impact on results from MODIS data (we used channel 21, offering a high dynamic range when channel 22 saturated) and the overestimation of Regarding the radiative power, values retrieved from SEVIRI data appeared underestimated when compared to temporally coincident MODIS observations. In addition, even with filtering MODIS data for low values of the satellite zenith angle, to take the poor viewing conditions (e.g., [11,20,70,71]) into account, differences in terms of VRP recorded especially during the 17th paroxysm were particularly high (≈12.5 GW). Those differences were less significant when weak thermal activities occurred. The MIR channel saturation, having a lower impact on results from MODIS data (we used channel 21, offering a high dynamic range when channel 22 saturated) and the overestimation of cloudy pixels observed on a number of SEVIRI scenes, due to the used cloud-masking procedure, were the main factors affecting the estimates of VRP. It should be mentioned, however, that the radiative power may be underestimated also because of the digital filter and the geometric resampling applied to SEVIRI raw data before their delivery via EUMETCast [72]. An overestimation of this parameter may instead occur due to the Point Spread Function (PSP), which is so large that even a small lava flow within a single pixel may generate significant effects on the measured radiance [73]. Finally, the impact of atmospheric effects on estimates of radiative power should be also considered (e.g., we recorded an increase in VRP up to~1 GW from SEVIRI data by setting a mean τ MIR value of 0.82, for a mid-latitude winter atmosphere).
Despite those factors, SEVIRI allowed us to characterise the Mt. Etna paroxysms also in terms of radiated energy. Figure 9 displays the plot of this parameter, which was calculated by time integrating the values of radiative power for each paroxysmal event, and the mean cloudiness values (dotted line). Based on the information provided by SEVIRI, the eruptive events of 22-23 February (the 5th paroxysm) and 31 March-1 April (the 17th paroxysm) were the most energetic (the 8th paroxysm was not quantified in terms of radiated energy also because of cloudy pixel overestimation). Therefore, by combining information from Figure 9 to AVHRR and MODIS observations, we can assert that the 17th paroxysm having the longest time duration was the most powerful, and probably radiated the largest amount of energy. These results confirm the importance of using a multi-sensor approach in investigating and quantifying thermal volcanic activity. the eruptive events of 22-23 February (the 5th paroxysm) and 31 March-1 April (the 17th paroxysm) were the most energetic (the 8th paroxysm was not quantified in terms of radiated energy also because of cloudy pixel overestimation). Therefore, by combining information from Figure 9 to AVHRR and MODIS observations, we can assert that the 17th paroxysm having the longest time duration was the most powerful, and probably radiated the largest amount of energy. These results confirm the importance of using a multi-sensor approach in investigating and quantifying thermal volcanic activity. From a volcanological point of view, between mid-February and the first days of April 2021, Mt. Etna produced 17 paroxysmal eruptions in just 44 days. Short but violent eruptions (capable of emitting a total of 40-50 million cubic meters of magma) occurred, which, on the whole, can be assimilated almost to a single eruptive event of considerable size: huge pyroclastic columns elevated over ten kilometres to the upper limit of the troposphere; lava flows 3-4 km long distributed on the east and south high flanks. Those eruptions were produced by the SEC, which doubled in size between 2011 and 2015. However, although summit eruptions are the most frequent at Mt. Etna, as shown in Figure 10, displaying the map of the areas more frequently affected by thermal anomalies retrieved from L8-OLI and S2-MSI data of 2013-2021, they do not represent the most dangerous eruptive events for this volcano. From a volcanological point of view, between mid-February and the first days of April 2021, Mt. Etna produced 17 paroxysmal eruptions in just 44 days. Short but violent eruptions (capable of emitting a total of 40-50 million cubic meters of magma) occurred, which, on the whole, can be assimilated almost to a single eruptive event of considerable size: huge pyroclastic columns elevated over ten kilometres to the upper limit of the troposphere; lava flows 3-4 km long distributed on the east and south high flanks. Those eruptions were produced by the SEC, which doubled in size between 2011 and 2015. However, although summit eruptions are the most frequent at Mt. Etna, as shown in Figure 10, displaying the map of the areas more frequently affected by thermal anomalies retrieved from L8-OLI and S2-MSI data of 2013-2021, they do not represent the most dangerous eruptive events for this volcano. Indeed, lateral eruptions are potentially much more dangerous, because they open vents at low altitudes, and therefore near or inside cities, generating lava flows that can destroy large urban areas in a short time [22,23]. But summit eruptions can be the prelude to lateral eruptions and must be analysed with extreme care.
After the brief but dynamically violent lateral eruption of 24-27 December 2018 [12,24,74], Etna first produced a modest and discontinuous summit activity, which increased in the second half of 2020, culminating first in almost continuous Strombolian Indeed, lateral eruptions are potentially much more dangerous, because they open vents at low altitudes, and therefore near or inside cities, generating lava flows that can destroy large urban areas in a short time [22,23]. But summit eruptions can be the prelude to lateral eruptions and must be analysed with extreme care.
After the brief but dynamically violent lateral eruption of 24-27 December 2018 [12,24,74], Etna first produced a modest and discontinuous summit activity, which increased in the second half of 2020, culminating first in almost continuous Strombolian activity at the SEC (January 2021), and then in the 17 eruptive paroxysms of February-April 2021. Numerous other paroxysmal events started from the night of 18-19 May 2021 and not yet finished at the time of writing. This is a succession of events similar to that recorded during 2000-2001 when over 70 eruptive paroxysms were the prelude to two violent and dangerous flank eruptions [75,76].
The increase in the intensity of the summit eruptions mentioned above was perfectly recorded by the set of satellite observations, at different temporal, spatial, and spectral resolutions, analysed in this work. They alerted possibly in advance of imminent phenomena and provided information also about the time duration and the energy released during the different paroxysms. Both those aspects are of great relevance for imagining the possible eruptive scenarios.

Conclusions
The paroxysms of 16 February-1 April 2021 were among the most intense eruptive events occurring at Mt. Etna during the last decade, as indicated by the results of this work.
The latter indicates that the 17th paroxysm was the most powerful, with values of radiative power up to about 14 GW (e.g., we retrieved VRP values in the range 1-2 GW during the paroxysms of May 2016 [17]), and probably the most energetic of the time series. Additionally, based on the information retrieved from L8-OLI and S2-MSI data, the 1st and 3rd paroxysm generated the largest thermal anomaly at Mt. Etna after April 2013. However, this information needs to be better assessed, also because of the low temporal sampling of analysed data especially before June 2015, when S2-MSI data were not available.
This study demonstrates the efficiency of the used multi-platform system in monitoring and quantifying Mt. Etna activity, despite the impact of a number of factors (e.g., clouds/plumes, pixel saturation, satellite viewing geometry, and atmospheric effects) on the identification and mapping of volcanic thermal anomalies from space. This system may provide reliable information also about subtle phases of thermal unrest (e.g., those preceding intense eruptions), contributing to the volcanic risk mitigation also in other areas.

Acknowledgments:
The authors wish to thank EUMETSAT and the Aeronautica Militare Italiana for its support to have access to MSG-SEVIRI data used in this work.

Conflicts of Interest:
The authors declare no conflict of interest.    The paroxysm was preceded by an intense Strombolian activity occurring on March 30th. Lava fountains, (time duration 540 min), eruptive column (9-10 km a.s.l.), ash impact on the SSW sector, lava flows towards E, SE, SSE (Valle de Bove) and S, SW (southern flank) up to 3.3 km long down to 1820 m a.s.l.