Volcanic Clouds Characterization of the 2020–2022 Sequence of Mt. Etna Lava Fountains Using MSG-SEVIRI and Products’ Cross-Comparison

: From December 2020 to February 2022, 66 lava fountains (LF) occurred at Etna volcano (Italy). Despite their short duration (an average of about two hours), they produced a strong impact on human life, environment, and air trafﬁc. In this work, the measurements collected from the Spinning Enhanced Visible and InfraRed Imager (SEVIRI) instrument, on board Meteosat Second Generation (MSG) geostationary satellite, are processed every 15 min to characterize the volcanic clouds produced during the activities. In particular, a quantitative estimation of volcanic cloud top height (VCTH) and ash/ice/SO 2 masses’ time series are obtained. VCTHs are computed by integrating three different retrieval approaches based on coldest pixel detection, plume tracking, and HYSPLIT models, while particles and gas retrievals are realized simultaneously by exploiting the Volcanic Plume Retrieval (VPR) real-time procedure. The discrimination between ashy and icy pixels is carried out by applying the Brightness Temperature Difference (BTD) method with thresholds obtained by making speciﬁc Radiative Transfer Model simulations. Results indicate a VCTH variation during the entire period between 4 and 13 km, while the SO 2 , ash, and ice total masses reach maximum values of about 50, 100, and 300 Gg, respectively. The cumulative ash, ice, and SO 2 emitted from all the 2020–2022 LFs in the atmosphere are about 750, 2300, and 670 Gg, respectively. All the retrievals indicate that the overall activity can be grouped into 3 main periods in which it passes from high (December 2020 to March 2021), low (March to June 2021), and medium/high (June 2021 to February 2022). The different products have been validated by using TROPOspheric Monitoring Instrument (TROPOMI) polar satellite sensor, Volcano Observatory Notices for Aviation (VONA) bulletins, and by processing the SEVIRI data considering a different and more accurate retrieval approach. The products’ cross-comparison shows a generally good agreement, except for the SO 2 total mass in case of high ash/ice content in the volcanic cloud.


Introduction
By releasing large quantities of particles and gases into the atmosphere, volcanic eruptions can have a significant impact on human health [1,2], the environment [3][4][5][6], and climate [7][8][9][10][11] and pose a severe threat to aviation safety [12].The residence time in the atmosphere of the emitted particles depends on their sizes and the height at which they are ejected.Typically, particles with a radius of less than about 20 µm can remain in the atmosphere for weeks and travel thousands of kilometers downwind.In addition to the particles, the most abundant gases are H 2 O, CO 2 , and SO 2 [13,14].The water vapor in combination with the ash particles, which act as condensation nuclei, can form water particles that, under particular pressure and temperature conditions, turn into ice.Volcanic ice clouds are spectrally indistinguishable from ice weather clouds but can be even more dangerous for aircraft as they can hide possible ash layers.
Etna (Sicily, Italy, 37.7 • N; 15.0 • E, Figure 1a), one of the most active volcanoes in the world, has exhibited since 2011 intense explosive activities ranging from Strombolian episodes to short lava fountain (LF) events [15][16][17].In 2016, this eruptive regime gradually transitioned to mild Strombolian explosive activities of long duration associated with isolated episodes of lava flows from the volcano's summit craters [18].The Christmas 2018 Etna eruption was preceded by a period of moderate explosive activity and small lava flows at the summit craters [19,20].
Remote Sens. 2023, 15, x FOR PEER REVIEW 2 of 22 particles that, under particular pressure and temperature conditions, turn into ice.Volcanic ice clouds are spectrally indistinguishable from ice weather clouds but can be even more dangerous for aircraft as they can hide possible ash layers.Etna (Sicily, Italy, 37.7°N; 15.0°E, Figure 1a), one of the most active volcanoes in the world, has exhibited since 2011 intense explosive activities ranging from Strombolian episodes to short lava fountain (LF) events [15][16][17].In 2016, this eruptive regime gradually transitioned to mild Strombolian explosive activities of long duration associated with isolated episodes of lava flows from the volcano's summit craters [18].The Christmas 2018 Etna eruption was preceded by a period of moderate explosive activity and small lava flows at the summit craters [19,20].On 13 December 2020, after about 18 months of eruptive pause, the volcano entered a new phase characterized by several LF episodes, all of them from the New Southeast Crater (NSEC) [23].After three other events (21 and 22 December 2020 and 18 January 2021), from 16 February 2021, Etna produced a series of intense and frequent LFs (every ~1-2 days) characterized by several kilometers' high volcanic clouds and ash fallout in the surrounding areas of the volcano.This had a strong impact on viability, stability of the roofs, air traffic (the Catania airport is a major international hub), agriculture, water contamination, and on the health of the local population [23].This first strong phase, with eight episodes in February and eleven episodes in March, finished with the 31 March-1 April 2021 LF (Figure 1b).After that, Mt.Etna took a break of about 50 days, and then a new event occurred on 19 May 2021.From this moment, the activity regained vigor with numerous fountains (eleven events in May, seventeen in June, seven in July) quite close in time but usually less energetic than the previous ones.Finally, the 2021 sequence came to an end with two episodes in August (9 and 29), one in September (21), and the last one on 23 October 2021.In 2022 only two LFs were observed on 10 and 21 February.
In this work, the data collected by the Spinning Enhanced Visible and InfraRed Imager (SEVIRI) instrument on board Meteosat Second Generation (MSG-4, Meteosat-11) geostationary satellite were used to fully characterize the volcanic clouds produced by all the different LFs by deriving their geometry (height and extension) and content (ash, ice, and SO2).The complete list of LFs from December 2020 to February 2022 was derived from the Volcano Observatory Notices for Aviation (VONA) bulletin emanated from Istituto Nazionale di Geofisica e Vulcanologia-Osservatorio Etneo (INGV-OE) when the wording "lava fountain" was present.Overall, 59 lava fountains (see Table A1) were processed; regarding the 66 episodes reported in Calvari and Nunnari [24], we considered as a single LF the three close episodes of 13-14 December 2020 and the two close episodes of 22-23 On 13 December 2020, after about 18 months of eruptive pause, the volcano entered a new phase characterized by several LF episodes, all of them from the New Southeast Crater (NSEC) [23].After three other events (21 and 22 December 2020 and 18 January 2021), from 16 February 2021, Etna produced a series of intense and frequent LFs (every ~1-2 days) characterized by several kilometers' high volcanic clouds and ash fallout in the surrounding areas of the volcano.This had a strong impact on viability, stability of the roofs, air traffic (the Catania airport is a major international hub), agriculture, water contamination, and on the health of the local population [23].This first strong phase, with eight episodes in February and eleven episodes in March, finished with the 31 March-1 April 2021 LF (Figure 1b).After that, Mt.Etna took a break of about 50 days, and then a new event occurred on 19 May 2021.From this moment, the activity regained vigor with numerous fountains (eleven events in May, seventeen in June, seven in July) quite close in time but usually less energetic than the previous ones.Finally, the 2021 sequence came to an end with two episodes in August (9 and 29), one in September (21), and the last one on 23 October 2021.In 2022 only two LFs were observed on 10 and 21 February.
In this work, the data collected by the Spinning Enhanced Visible and InfraRed Imager (SEVIRI) instrument on board Meteosat Second Generation (MSG-4, Meteosat-11) geostationary satellite were used to fully characterize the volcanic clouds produced by all the different LFs by deriving their geometry (height and extension) and content (ash, ice, and SO 2 ).The complete list of LFs from December 2020 to February 2022 was derived from the Volcano Observatory Notices for Aviation (VONA) bulletin emanated from Istituto Nazionale di Geofisica e Vulcanologia-Osservatorio Etneo (INGV-OE) when the wording "lava fountain" was present.Overall, 59 lava fountains (see Table A1) were processed; regarding the 66 episodes reported in Calvari and Nunnari [24], we considered as a single LF the three close episodes of 13-14 December 2020 and the two close episodes of [22][23] February 2021, hardly distinguishable from SEVIRI.Furthermore, the four episodes of 27-28 May 2021 were excluded due to the strong presence of meteorological clouds in the area which did not allow volcanic cloud discrimination and detection.
The paper is organized as follows: in Section 2, the satellite, the atmospheric data, and the methods used for the volcanic height and particles/gas retrievals computation are described.In Section 3, all the results of the 2020-2022 Etna lava fountains sequence are reported: volcanic cloud heights, Brightness Temperature Differences (BTD) thresholds used for ash/ice particle discrimination, and total masses' time series of ash/ice/SO 2 .In Section 3.4, the volcanic cloud heights and total masses are compared with the results obtained by other sensors, and a different retrieval procedure is applied to SEVIRI data.Finally, the conclusions are outlined in Section 4. In the Appendix are reported: a table summarizing all the LFs, and for each, the main results presented in the paper, the comparison with the LF heights and ground deposits, and a list of the acronyms used.As Supplementary Materials, the complete time series plots of VPR total masses of ash/ice/SO 2 obtained from the processing of more than two thousand SEVIRI images are displayed.

Materials and Methods
In this work, the data collected by SEVIRI were considered.SEVIRI has 12 spectral channels from visible (VIS) to Thermal InfraRed (TIR), a nadir spatial resolution of 3 km at sub-satellite point, and a temporal resolution of 5 or 15 min (Rapid Scan or Earth Full Disk, respectively).A satellite acquisition system (Multimission Acquisition SysTem, MAST) has been developed at INGV to collect and pre-process the SEVIRI images in real time [25].All the SEVIRI images used in this work have been resized for a "Central Mediterranean Area" and resampled on a regular grid of 3 × 3 km 2 .For each LF, the image processing starts when the eruption begins until the volcanic cloud is no longer detectable (by dilution or by leaving the considered area).

Volcanic Cloud Top Height
The Volcanic Cloud Top Height (VCTH) is one of the most important parameters for aircraft security [25,[27][28][29]].However, it also represents a key input for volcanic cloud retrieval procedures [30,31] and for a correct initialization of the volcanic ash dispersion and deposition models [32][33][34][35].In this work, we compare the VCTH obtained from 3 different methods, all applied to SEVIRI images: darkest pixel, cloud tracking, and HYSPLIT.

Darkest Pixel Method (DP)
The consolidated "Darkest Pixel" (DP) procedure is based on the comparison between the minimum SEVIRI 10.8 µm brightness temperature (BT 10.8 ) of the pixels contained in a fixed area over the summit craters and the atmospheric temperature profile measured in the same area and at the same time of satellite acquisition [30,36,37].In this work, an area of 17 × 17 pixels around Etna (about 50 × 50 km 2 ) was chosen and the maximum value of the top altitude time series (obtained from several consecutive SEVIRI images) was taken as VCTH.The method is very simple and usually permits obtaining the VCTH with good accuracy.However, the DP method also has several weaknesses.The main one happens when the plume is not completely opaque.In this case, the radiance from the surface and the lower part of the atmosphere beneath the plume increases the top-of-atmosphere (TOA) radiance measured by satellite, leading to underestimated plume heights.To consider the possible non-complete opacity of the pixel, the SEVIRI 10.8 µm brightness temperatures were decreased by 2 K as suggested by [36].Another drawback of the method is due to the atmospheric temperature vertical behavior: the uncertainty on VCTH retrieval increases significantly around the tropopause because of the low temperature gradient and, due to the rise of stratospheric temperature, this often leads to double solutions (see Figure 2a).In this work, the lowest height has always been chosen, because injections of the plume into the stratosphere are quite unusual for the characteristics of the 2020-2022 Etna activity.The last drawback of the DP method is related to the possible thermal disequilibrium with the surrounding atmosphere.In this study, nine cases with a cloud top temperature lower than the minimum air temperature were found, resulting in no point of intersection with the atmospheric temperature profiles.For these specific cases, the nearest value (tropopause height) was considered as VCTH.
main one happens when the plume is not completely opaque.In this case, the radiance from the surface and the lower part of the atmosphere beneath the plume increases the top-of-atmosphere (TOA) radiance measured by satellite, leading to underestimated plume heights.To consider the possible non-complete opacity of the pixel, the SEVIRI 10.8 μm brightness temperatures were decreased by 2 K as suggested by [36].Another drawback of the method is due to the atmospheric temperature vertical behavior: the uncertainty on VCTH retrieval increases significantly around the tropopause because of the low temperature gradient and, due to the rise of stratospheric temperature, this often leads to double solutions (see Figure 2a).In this work, the lowest height has always been chosen, because injections of the plume into the stratosphere are quite unusual for the characteristics of the 2020-2022 Etna activity.The last drawback of the DP method is related to the possible thermal disequilibrium with the surrounding atmosphere.In this study, nine cases with a cloud top temperature lower than the minimum air temperature were found, resulting in no point of intersection with the atmospheric temperature profiles.For these specific cases, the nearest value (tropopause height) was considered as VCTH.
The uncertainties of the DP method were estimated by considering +/−2 K of the value obtained from the coldest pixel temperature minus 2 K [36].The mean uncertainty for all the LFs considered is +/−0.3km. Figure 2 shows an example of the VCTH estimation from the DP procedure applied on the 12 March 2021 LF.In the left plot, the volcanic cloud's coldest pixel temperature value with its uncertainties (blue vertical solid and dashed lines, respectively), obtained from the SEVIRI image collected at 09:25 UTC, is compared with the NCEP temperature profile collected in the same region at the closest time (12 March 2021 12 UTC at 37.5°N, 15.0°E, red solid line).In the right plot, the VCTH time series obtained by processing the different SEVIRI images of the same day every 15 min is displayed.The NCEP temperature profiles of 12 March 2021 at 06 and 12 UTC were used for SEVIRI images before and after 9 UTC, respectively.

Cloud Tracking Method (CT)
The high data frequency of SEVIRI images can be exploited to retrieve wind speed and direction of the volcanic clouds for each event [37,38].First of all, it is necessary to perform the detection of the volcanic cloud (see Section 2.2).In case of an intense but shortlived eruption (such as LFs), it is quite easy to find and track, in the series of images, the pixels with the minimum radiance at 10.8 μm, which corresponds to the maximum ash/ice concentration.Computing the distance from the top of the volcano for at least 2-3 h from the start of the eruption and using a linear fit (blue line), the speed of the volcanic cloud is obtained (see Figure 3a).Finally, by comparing the retrieved peak speed with the wind The uncertainties of the DP method were estimated by considering +/−2 K of the value obtained from the coldest pixel temperature minus 2 K [36].The mean uncertainty for all the LFs considered is +/−0.3km. Figure 2 shows an example of the VCTH estimation from the DP procedure applied on the 12 March 2021 LF.In the left plot, the volcanic cloud's coldest pixel temperature value with its uncertainties (blue vertical solid and dashed lines, respectively), obtained from the SEVIRI image collected at 09:25 UTC, is compared with the NCEP temperature profile collected in the same region at the closest time (12 March 2021 12 UTC at 37.5 • N, 15.0 • E, red solid line).In the right plot, the VCTH time series obtained by processing the different SEVIRI images of the same day every 15 min is displayed.The NCEP temperature profiles of 12 March 2021 at 06 and 12 UTC were used for SEVIRI images before and after 9 UTC, respectively.

Cloud Tracking Method (CT)
The high data frequency of SEVIRI images can be exploited to retrieve wind speed and direction of the volcanic clouds for each event [37,38].First of all, it is necessary to perform the detection of the volcanic cloud (see Section 2.2).In case of an intense but short-lived eruption (such as LFs), it is quite easy to find and track, in the series of images, the pixels with the minimum radiance at 10.8 µm, which corresponds to the maximum ash/ice concentration.Computing the distance from the top of the volcano for at least 2-3 h from the start of the eruption and using a linear fit (blue line), the speed of the volcanic cloud is obtained (see Figure 3a).Finally, by comparing the retrieved peak speed with the wind speed of an atmospheric profile collected at the same time and position (red line, Figure 3b), the volcanic cloud altitude can be derived.It is important to note that, unlike the DP method, this procedure does not necessarily produce the maximum height (VCTH), but depending on the location of the center of mass, a lower altitude may result.Unfortunately, due to the characteristics of atmospheric profiles, much more than a single intersection point can often be found.When this happens, the use of wind direction can provide help, even though sometimes only a range of height can be obtained; in these cases, the mean value was considered.Other sources of uncertainty are the non-linearity of the plume trajectory due to the non-uniformity of the wind field, which can produce large errors with this method.
speed of an atmospheric profile collected at the same time and position (red line, Figure 3b), the volcanic cloud altitude can be derived.It is important to note that, unlike the DP method, this procedure does not necessarily produce the maximum height (VCTH), but depending on the location of the center of mass, a lower altitude may result.Unfortunately, due to the characteristics of atmospheric profiles, much more than a single intersection point can often be found.When this happens, the use of wind direction can provide help, even though sometimes only a range of height can be obtained; in these cases, the mean value was considered.Other sources of uncertainty are the non-linearity of the plume trajectory due to the non-uniformity of the wind field, which can produce large errors with this method.
The uncertainties of the CT method were computed considering the standard error of the estimate of the linear fit (cyan lines of Figure 3).The mean uncertainty for all the LFs considered is +/−0.7 km.

HYSPLIT Method (HY)
Finally, a similar procedure based on the Hybrid Single-Particle Lagrangian Integrated Trajectory (HYSPLIT) model has been also considered [39][40][41].By plotting several forward trajectories (obtained from the GFS-0.25°archive) at different altitudes starting from the top of Etna at the time of the eruption start, it is possible to visualize which one better corresponds to the volcanic cloud position for a SEVIRI image collected some hours later.In its basic form, the HYSPLIT model does not consider the particulate transport, so formally it would only be suitable for deriving the height of the SO2 cloud.Regardless, the HYSPLIT trajectories tool is a model widely used for the estimation of the volcanic cloud altitudes under the main hypothesis that ash, ice and SO2 are totally collocated [37,42,43].In this case, too, the characteristics of wind speed and direction can produce large uncertainties if there is not a marked vertical wind shear.On the other hand, the dispersion of the volcanic cloud can often highlight the different heights of the various parts of the plume.In both these cases it was assumed the maximum height was VCTH.The forward trajectories were chosen at 1 km altitude steps, so the uncertainty of the HY method was assumed +/−0.5 km. Figure 4 shows an example of the HY method applied on the 9 August 2021 Etna LF.The 12 h forward trajectories start from 02:00 UTC (beginning of the eruption) at 8 km (red), 9 km (blue), and 10 km (green) above sea level.The corresponding SEVIRI image (at 14:00 UTC) highlights the presence of the volcanic cloud (contoured in yellow) in the Mediterranean Sea, between 35°N-36°N and 18°E-20°E.From the comparison between the trajectories and the volcanic cloud position, the VCTH retrieved is about 9 km.The uncertainties of the CT method were computed considering the standard error of the estimate of the linear fit (cyan lines of Figure 3).The mean uncertainty for all the LFs considered is +/−0.7 km.

HYSPLIT Method (HY)
Finally, a similar procedure based on the Hybrid Single-Particle Lagrangian Integrated Trajectory (HYSPLIT) model has been also considered [39][40][41].By plotting several forward trajectories (obtained from the GFS-0.25 • archive) at different altitudes starting from the top of Etna at the time of the eruption start, it is possible to visualize which one better corresponds to the volcanic cloud position for a SEVIRI image collected some hours later.In its basic form, the HYSPLIT model does not consider the particulate transport, so formally it would only be suitable for deriving the height of the SO 2 cloud.Regardless, the HYSPLIT trajectories tool is a model widely used for the estimation of the volcanic cloud altitudes under the main hypothesis that ash, ice and SO 2 are totally collocated [37,42,43].In this case, too, the characteristics of wind speed and direction can produce large uncertainties if there is not a marked vertical wind shear.On the other hand, the dispersion of the volcanic cloud can often highlight the different heights of the various parts of the plume.In both these cases it was assumed the maximum height was VCTH.The forward trajectories were chosen at 1 km altitude steps, so the uncertainty of the HY method was assumed +/−0.5 km. Figure 4 shows an example of the HY method applied on the 9 August 2021 Etna LF.The 12 h forward trajectories start from 02:00 UTC (beginning of the eruption) at 8 km (red), 9 km (blue), and 10 km (green) above sea level.The corresponding SEVIRI image (at 14:00 UTC) highlights the presence of the volcanic cloud (contoured in yellow) in the Mediterranean Sea, between 35 • N-36 • N and 18 • E-20 • E. From the comparison between the trajectories and the volcanic cloud position, the VCTH retrieved is about 9 km.

Volcanic Cloud Detection and Discrimination
In this work, the volcanic cloud detection was performed by visual inspection of an RGB composition obtained using a combination of the BTs of the channels centered at 8.7, 10.8, and 12 μm (BT8.7,BT10.8,BT12) that allows identifying both the volcanic cloud particles and SO2.For each SEVIRI image, a "plume mask" was obtained drawing a region of interest (ROI) around the volcanic cloud.Once the cloud has been identified, it is also necessary to discriminate the types of particles from which it is mainly composed.For this task, the BTD approach was used [44,45].This method allows discriminating between ash and ice/water particles by exploiting the different spectral absorption in the TIR spectral range.In this interval, the absorption of ash particles at 10.8 μm is larger than that at 12 μm.The opposite happens for ice/water clouds, which absorb more significantly at longer TIR wavelengths.Therefore, the BTD, defined as the difference between the brightness temperature computed at 10.8 and 12 μm (BT10.8−BT12),turns out to be generally negative (BTD < 0) for the region affected by ash and positive (BTD > 0) for the region containing ice/water clouds.In almost all the LFs considered, the formation of high quantities of water/ice particles (strongly positive BTD values) has been observed.Due to the high altitudes generally reached by the LF volcanic clouds, we considered the formation of ice much more likely than the presence of water droplets in liquid form.Figure 5 shows an example of RGB composition and BTD.In the RGB image (panel a), the presence of volcanic ash, ice, and SO2 is identified by the red, dark, and light green colors, respectively, while the b and w BTD image (panel b) indicates the existence of both ash (dark part, negative values) and ice (white part, positive values) particles.

Volcanic Cloud Detection and Discrimination
In this work, the volcanic cloud detection was performed by visual inspection of an RGB composition obtained using a combination of the BTs of the channels centered at 8.7, 10.8, and 12 µm (BT 8.7 , BT 10.8 , BT 12 ) that allows identifying both the volcanic cloud particles and SO 2 .For each SEVIRI image, a "plume mask" was obtained drawing a region of interest (ROI) around the volcanic cloud.Once the cloud has been identified, it is also necessary to discriminate the types of particles from which it is mainly composed.For this task, the BTD approach was used [44,45].This method allows discriminating between ash and ice/water particles by exploiting the different spectral absorption in the TIR spectral range.In this interval, the absorption of ash particles at 10.8 µm is larger than that at 12 µm.The opposite happens for ice/water clouds, which absorb more significantly at longer TIR wavelengths.Therefore, the BTD, defined as the difference between the brightness temperature computed at 10.8 and 12 µm (BT 10.8 −BT12), turns out to be generally negative (BTD < 0) for the region affected by ash and positive (BTD > 0) for the region containing ice/water clouds.In almost all the LFs considered, the formation of high quantities of water/ice particles (strongly positive BTD values) has been observed.Due to the high altitudes generally reached by the LF volcanic clouds, we considered the formation of ice much more likely than the presence of water droplets in liquid form.Figure 5 shows an example of RGB composition and BTD.In the RGB image (panel a), the presence of volcanic ash, ice, and SO 2 is identified by the red, dark, and light green colors, respectively, while the b and w BTD image (panel b) indicates the existence of both ash (dark part, negative values) and ice (white part, positive values) particles.
The BTD approach is effective and simple to apply, even if it can lead to false positive ash detections (pixels wrongly detected as ashy) in some cases such as on clear surfaces during the night, on deserts, on very cold or ice surfaces, etc. [46].Otherwise, the BTD can lead to false negative ash detection (pixels wrongly detected as non-ashy) that may arise in case of high-water vapor content.This condition, frequent on Etna volcano, can hide and then cancel out the ash particles' effects on the BTD, thus revealing fewer ashy pixels than those that exist [46,47].For this reason, a correction, based on Radiative Transfer Models (RTM) computations, has to be applied [36,47].Figure 6a shows an example of the inverted arches curves obtained from ice (cyan lines) and ash pumice type (black lines, [48]) using MODTRAN 5.3 [49].In this plot, the Aerosol Optical Depth at 550 nm (AOD 0.55 ) varies from 0 (right, bigger BT 10.8 values) to 10 (left, lower BT 10.8 values) for both ash and ice, while Effective Radius (R e ) varies from 0.55 to 10 µm (down to up) for ash and from 1.39 to 50 µm (up to down) for ice.To obtain a BTD threshold to use for discriminating between ash and ice, a quite good approximation is to consider a horizontal linear fit (red line) of the points with maximum R e (10 µm for ash and 50 µm for ice).Figure 6b shows the trend of the BTD thresholds (obtained as just described) as a function of the months of the year, VCTH The BTD approach is effective and simple to apply, even if it can lead to false positive ash detections (pixels wrongly detected as ashy) in some cases such as on clear surfaces during the night, on deserts, on very cold or ice surfaces, etc. [46].Otherwise, the BTD can lead to false negative ash detection (pixels wrongly detected as non-ashy) that may arise in case of high-water vapor content.This condition, frequent on Etna volcano, can hide and then cancel out the ash particles' effects on the BTD, thus revealing fewer ashy pixels than those that exist [46,47].For this reason, a correction, based on Radiative Transfer Models (RTM) computations, has to be applied [36,47].Figure 6a shows an example of the inverted arches curves obtained from ice (cyan lines) and ash pumice type (black lines, [48]) using MODTRAN 5.3 [49].In this plot, the Aerosol Optical Depth at 550 nm (AOD0.55)varies from 0 (right, bigger BT10.8 values) to 10 (left, lower BT10.8values) for both ash and ice, while Effective Radius (Re) varies from 0.55 to 10 μm (down to up) for ash and from 1.39 to 50 μm (up to down) for ice.To obtain a BTD threshold to use for discriminating between ash and ice, a quite good approximation is to consider a horizontal linear fit (red line) of the points with maximum Re (10 μm for ash and 50 μm for ice).Figure 6b shows the trend of the BTD thresholds (obtained as just described) as a function of the months of the year, VCTH   The BTD approach is effective and simple to apply, even if it can lead to false positive ash detections (pixels wrongly detected as ashy) in some cases such as on clear surfaces during the night, on deserts, on very cold or ice surfaces, etc. [46].Otherwise, the BTD can lead to false negative ash detection (pixels wrongly detected as non-ashy) that may arise in case of high-water vapor content.This condition, frequent on Etna volcano, can hide and then cancel out the ash particles' effects on the BTD, thus revealing fewer ashy pixels than those that exist [46,47].For this reason, a correction, based on Radiative Transfer Models (RTM) computations, has to be applied [36,47].Figure 6a shows an example of the inverted arches curves obtained from ice (cyan lines) and ash pumice type (black lines, [48]) using MODTRAN 5.3 [49].In this plot, the Aerosol Optical Depth at 550 nm (AOD0.55)varies from 0 (right, bigger BT10.8 values) to 10 (left, lower BT10.8values) for both ash and ice, while Effective Radius (Re) varies from 0.55 to 10 μm (down to up) for ash and from 1.39 to 50 μm (up to down) for ice.To obtain a BTD threshold to use for discriminating between ash and ice, a quite good approximation is to consider a horizontal linear fit (red line) of the points with maximum Re (10 μm for ash and 50 μm for ice).Figure 6b   As Figure 6b clearly shows, the combination of the various parameters (PTH, SST, VCTH, and VZA) can contribute significantly to the variation of the BTD threshold (from 0.6 to 1.7).For this reason, the BTD values used in this work were obtained from RTM simulations specific for SEVIRI (MSG4) for each 2020-2022 LF episode (see Section 3.2).

Particles and SO 2 Retrieval Method
The Volcanic Plume Retrieval (VPR) procedure is a linearization of the radiative transfer equation capable to retrieve from multispectral satellite images the AOD 0.55 , the R e , the columnar abundance (m), and the total mass (M) of particles (ash, ice, water droplets, etc.) contained in a volcanic cloud [31,38,50].The channels used for the volcanic particles' retrievals are those centered at 10.8 µm and 12 µm, while for the SO 2 estimation, the channel centered at 8.7 µm is considered.Because of the ash/ice particles' absorption in the whole TIR spectral range, their quantitative estimations are taken into account to correct for the SO 2 amounts.The lack of this correction would otherwise lead to a significant overestimation of the amount of SO 2 present in the volcanic cloud [51].
The only VPR input required at run time is the plume temperature which can be easily obtained from a relevant atmospheric profile knowing the plume altitude.Before running, VPR also needs a "plume mask" so volcanic cloud detection and discrimination (ash or ice) is necessary.The main advantage of the VPR procedure is that it is easy to use and very fast.Coupled with SEVIRI's high temporal resolution, it allows for quick and reliable volcanic cloud retrievals during both day and night.
Several comparisons have been made in the past between the VPR and other procedures based on the calculation of atmospheric corrections, showing good agreement [22,52].Among the particle's properties, R e is that with the lower accuracy, while a better performance is related to AOD 0.55 and M. When no particles are present in the volcanic cloud or if the aerosol transmittance would be perfectly known, the SO 2 accuracy is very high; unfortunately, this excellent result is reduced in most real cases since the SO 2 retrieval is highly dependent on that of the ash/ice particles.In particular, the presence of high quantities of ice can produce greater SO 2 uncertainties [50].In this work, we considered an overall error of ±40% for particles (ash/ice) and ±50% for SO 2 total mass, due to the presence of large amounts of ice.The column height reported in the INGV VONA bulletin (www.ct.ingv.it/index.php/monitoraggio-e-sorveglianza/prodotti-del-monitoraggio/comunicati-vona, accessed on 16 March 2022) was considered [37,53].This value is predominantly obtained from visible surveillance cameras placed around Etna (during daytime and with good visibility conditions) or, more occasionally, from satellites (darkest pixel procedure).Unfortunately, in some cases, the VCTH is not given, while when more than one value was present (due to LF progress), the maximum was considered.The uncertainty of the VONA column height was set to ±0.5 km, according to Scollo et al. [53].

TROPOMI
The TROPOspheric Monitoring Instrument (TROPOMI) is a spectrometer on board Sentinel-5 Precursor (S5P) polar orbit satellite that covers a spectral range from ultraviolet (UV) to short wave infrared (SWIR), with a spatial resolution of 5.5 × 3.5 km 2 and a revisit time of about 1 day [54].The UV channels were used for the validation of volcanic cloud altitude and SO 2 .Particularly significant for TROPOMI is the ability to detect SO 2 columnar abundance of about 0.02 g/m 2 (0.7 DU), i.e., about thirty times better than the sensitivity of the multispectral satellite sensors [22].The plume height of SO 2 was retrieved from TROPOMI measurements using the algorithm of [55].In brief, the radiance measurements are analyzed using a Covariance-Based Retrieval Algorithm (COBRA) in which the SO 2 vertical column and plume height are jointly retrieved based on a look-up table approach.The SO 2 total mass is obtained by summing all contributions from the pixels of the plume.Note that by doing so, some pixels are rejected (those for which the height retrieval did not converge, essentially for vertical columns < 5 DU), and therefore the SO 2 total mass is probably a lower estimate.The uncertainty of the TROPOMI height and SO 2 total mass was set to +/−2 km [55] and 35% [22], respectively, an error estimate which might be too optimistic for conditions with a lot of volcanic ash.

LUT p
The Look-Up Tables procedure (LUT p ) is a well-known particles and gases retrieval method [30,45,47,51,56].It uses the same three TIR SEVIRI channels as VPR (8.7, 10.8, and 12 µm) and, through linear and bilinear interpolation with specific RTM simulations, permits to obtain particle AOD 0.55 , R e , mass, and SO 2 mass of a volcanic cloud.Generally, the LUT p is considered more accurate than the VPR [52], but this comes at the cost of a greater number of input parameters and a much longer calculation.For validating the VPR results, the LUT p was applied to the same SEVIRI images using the same RTM simulations used for the BTD thresholds estimation (see Section 3.2) plus 11 SO 2 values (from 0 to 10 g/m 2 , step 1 g/m 2 ).An overall error of ±40% for particles and SO 2 total mass was considered [21,30].

Volcanic Cloud Top Height
Figure 7 shows the VCTH for all the 2020-2022 Etna LFs using the three different procedures described in Section 2.1 (DP, CT, and HY).Each method has strengths and weaknesses depending on the conditions, so to get the "best estimate", a weighted average is presented (see Equations ( 1)-( 3)).The choice of the weights derives from the intrinsic characteristics of the three methods and from our experience: DP is considered the most reliable for high altitudes (strong eruptions with completely opaque pixels), while HY turns out better for low altitudes (weak eruptions).Because the CT method and the HYSPLIT model use a single atmospheric profile and a complete set of meteorological data, respectively, the latter is considered generally more reliable.The overall VCTH uncertainty was obtained from the combination of the uncertainties (according to Equations ( 1 1)-( 3)) used as input in the VPR procedure.The "Date" labels (x-axis) are in the format year-month-day, with letters a and b to distinguish the LFs that occurred on the same day.
Generally, the three methods agree well with each other: the mean absolute difference is 1.1, 1.2, and 1.4 km for DP vs. CT, DP vs. HY, and CT vs. HY, respectively.The VCTH ranges from 4 km (18 January 2021) to 13 km (23 October 2021).As the figure shows, the results emphasize three main periods (highlighted by the background colors of the plot): after an initial phase with low altitudes volcanic cloud (except for 21 December 2020), the period from 16 February to 19 March 2021 was characterized by high volcanic plumes with 10 episodes with VCTH greater than 10 km.Then, a long phase (3 months) occurred with 15 episodes with lower VCTH that ranges between 5-7 km, and from 19 June 2021, a new VCTH increase with several values greater than 10 km.Finally, this long eruptive phase ended with four strong (but spaced over a time range of 5 months) paroxysms from 21 September 2021 to 21 February 2022 with VCTH greater than 11 km.The complete list of VCTH values is given in Table A1.
It is also interesting to note that the retrieved SEVIRI-VCTH trend is in good agreement with the LF maximum heights trend shown by Calvari and Nunnari [24], obtained by exploiting the measurements collected from the ground-based TIR cameras (see Figure A1), with a Pearson's correlation coefficient of 0.65.The average difference between VCTH and LF maximum heights is 3.5 km (considering that NSEC is about 3.3 km asl high).

BTD Thresholds
The BTD threshold values used in this work to discriminate ashy and icy pixels were obtained from specific RTM simulations related to each LF considered as described in Section 2.2.The 2020-2022 atmospheric vertical profiles (PTH) at (37.5 • N, 15.0 • E) from NCEP/NCAR dataset were used; depending on eruption timing (nighttime or daytime), the 00 UTC or 12 UTC was chosen.From each of them, the total precipitable water (PW) was computed, too.SST values were obtained from "NOAA Daily Optimum Interpolation Sea Surface Temperature" (https://www.psl.noaa.gov/data/gridded/data.noaa.oisst.v2.highres.html,accessed on 19 July 2022 [57][58][59]).These data are provided at a 0.25 • × 0.25 • spatial resolution and were averaged over an area between 35 • -40 • N and 12.5 • -17.5 • E. Sea surface emissivity was fixed to 0.98.The VCTH weighted mean values (see Section 3.1, Table A1) were considered, with a constant plume thickness of 1 km.Pumice ash [48] and ice optical properties (extinction coefficient, single scattering albedo, and asymmetry parameter) derive directly from the application of the Mie theory considering a specific particle's refractive index and size distribution.The refractive indexes derive from the Aerosol Refractive Index Archive-ARIA database, compiled from the EODG group of Oxford University (http://eodg.atm.ox.ac.uk/ARIA, accessed on 27 October 2022), while the size distribution is considered log-normal.Regarding the VZA, a fixed value was considered for a geostationary SEVIRI sensor (VZA = 45 • , typical of the Etna area).In this way, for each LF considered, only one BTD threshold value was used.
In Figure 8c, all the SEVIRI-BTD values used are reported.They range between 0.48 and 1.83 with a mean (standard deviation) value of 1.04 (0.33).The major contribution to BTD variations is due to VCTH (Figure 8b), even if PW and SST (see Figure 8a) are significant, too.Sometimes similar VCTH values do not correspond to the same BTD values.For example, the 21 December 2020 and 10 February 2022 LFs have quite the same VCTH (11.4 and 11.5 km, respectively) and very different BTD (1.7 and 1.0, respectively).For SEVIRI, the Pearson's correlation coefficient between BTD thresholds and VCTH, PW, and SST is 0.73, 0.57, and 0.27, respectively.Given the wide temporal extension (more than 1 year) and the large range of volcanic cloud heights considered (from 4 to 13 km), the average value BTD = 1 can be considered a first good approximation for Etna, instead of the standard BTD = 0.
SST is 0.73, 0.57, and 0.27, respectively.Given the wide temporal extension (more than year) and the large range of volcanic cloud heights considered (from 4 to 13 km), the av erage value BTD = 1 can be considered a first good approximation for Etna, instead of th standard BTD = 0.

Total Mass of Ash, Ice, and SO2
Once the VCTH has been computed and the volcanic cloud detection and discrimi nation performed, the VPR procedure can be applied.From each satellite image, maps o ash/ice Re [μm], AOD0.55 and ma/mi/mSO2 [g/m 2 ] were computed.Finally, Ma/Mi/MSO2 [Gg total masses were easily calculated by adding the contribution of all the pixels multiplie by their area.Figure 9 shows an example of columnar abundance maps and time series o Ma/Mi/MSO2 for 23 October 2021 LF.The time trend highlights an almost simultaneou growth of both ash and ice from about 9 UTC until 10 UTC.After that, a decrease in as and an increase in ice quantities are observed, compatible with the progressive ice gener ation due to the ash particles that act as condensation nuclei.The SO2 presents a stron increase until about 10:15 UTC (when the LF ends, represented by the dashed red vertica line) and then a slight increase until 13:30 UTC.This latter effect is probably not a rea

Total Mass of Ash, Ice, and SO 2
Once the VCTH has been computed and the volcanic cloud detection and discrimination performed, the VPR procedure can be applied.From each satellite image, maps of ash/ice R e [µm], AOD 0.55 and m a /m i /m SO2 [g/m 2 ] were computed.Finally, M a /M i /M SO2 [Gg] total masses were easily calculated by adding the contribution of all the pixels multiplied by their area.Figure 9 shows an example of columnar abundance maps and time series of M a /M i /M SO2 for 23 October 2021 LF.The time trend highlights an almost simultaneous growth of both ash and ice from about 9 UTC until 10 UTC.After that, a decrease in ash and an increase in ice quantities are observed, compatible with the progressive ice generation due to the ash particles that act as condensation nuclei.The SO 2 presents a strong increase until about 10:15 UTC (when the LF ends, represented by the dashed red vertical line) and then a slight increase until 13:30 UTC.This latter effect is probably not a real increase in SO 2 but a consequence of the huge growth of the ice content that causes some criticalities in the retrieval correction procedure.The following decrease in SO 2 mass is due to the volcanic cloud dilution and chemical recombination.It should be noted that the dilution of the ice/ash particles (decrease in the total mass) usually begins earlier and is faster (the descent is steeper) than that of SO 2 , which remains in the atmosphere for a longer time.The complete VPR total masses' time series for all the LFs considered can be found in the Supplementary Materials (Figure S1).
In Figure 10, the VPR M a /M i /M SO2 for all the LFs were reported (in logarithmic scale).For each LF, the maximum value obtained from the time series was considered.To avoid the anomalous SO 2 increase after the end of LF, the SO 2 maximum value in the time range between the start and the end (+1 h) of the eruption was considered, in which the LF end time was obtained from Calvari and Nunnari [24].In most cases, the ice content is greater than ash but there are still some episodes where the opposite occurs, both with low and high plume heights.Together with M a /M i /M SO2 , Figure 10 also reports the VCTH to emphasize the general good correlation between the height and total masses as documented in Mastin et al. [60] where different studies on this topic are analyzed.The Pearson's correlation coefficient between VCTH and the logarithm of M x is 0.75 (ash), 0.64 (ice), and 0.57 (SO 2 ).Finally, note how the generation of ice is related not only to the volcanic cloud height but also to the season: during the summer (yellow background color), with almost the same plume height, the amount of ice is lower than that observed in winter (green background color).increase in SO2 but a consequence of the huge growth of the ice content that causes some criticalities in the retrieval correction procedure.The following decrease in SO2 mass is due to the volcanic cloud dilution and chemical recombination.It should be noted that the dilution of the ice/ash particles (decrease in the total mass) usually begins earlier and is faster (the descent is steeper) than that of SO2, which remains in the atmosphere for a longer time.The complete VPR total masses' time series for all the LFs considered can be found in the Supplementary Materials (Figure S1).In Figure 10, the VPR Ma/Mi/MSO2 for all the LFs were reported (in logarithmic scale).For each LF, the maximum value obtained from the time series was considered.To avoid the anomalous SO2 increase after the end of LF, the SO2 maximum value in the time range between the start and the end (+1 h) of the eruption was considered, in which the LF end time was obtained from Calvari and Nunnari [24].In most cases, the ice content is greater than ash but there are still some episodes where the opposite occurs, both with low and high plume heights.Together with Ma/Mi/MSO2, Figure 10 also reports the VCTH to emphasize the general good correlation between the height and total masses as documented in Mastin et al. [60] where different studies on this topic are analyzed.The Pearson's correlation coefficient between VCTH and the logarithm of Mx is 0.75 (ash), 0.64 (ice), and 0.57 (SO2).Finally, note how the generation of ice is related not only to the volcanic cloud height but also to the season: during the summer (yellow background color), with almost the same plume height, the amount of ice is lower than that observed in winter (green background color).The VPR total mass results have also been compared with the lava deposits on the ground obtained by Ganci et al. [61] exploiting the same SEVIRI images (Figure A2).The two different products show quite good agreement, with a Pearson's correlation coefficient between the lava volume and the logarithm of (Ma + Mi + MSO2) of 0.57.The VPR total mass results have also been compared with the lava deposits on the ground obtained by Ganci et al. [61] exploiting the same SEVIRI images (Figure A2).The two different products show quite good agreement, with a Pearson's correlation coefficient between the lava volume and the logarithm of (M a + M i + M SO2 ) of 0.57.

Products Cross-Comparison
Figure 11 shows the cross-comparison between SEVIRI, INGV-VONA bulletin, and TROPOMI VCTHs.Generally, TROPOMI heights are lower than those obtained from SEVIRI and surveillance cameras.However, there is a fundamental difference: the VCTH-DP and VONA estimates are related to aerosols emissions near the volcano (at the beginning of the eruption), while TROPOMI height refers to SO 2 which corresponds to a rather welldispersed plume at quite some distance from Etna.Dilution of the cloud is certainly important in this case.Moreover, the effect of ash on the measurement sensitivity of TROPOMI can lead to biases in the retrieved SO 2 heights.Overall and considering all these different characteristics, the agreement of the various trends is quite good.Quantifying the differences between SEVIRI-VCTH and TROPOMI height, we find large discrepancies (>3 km) only for a few LFs (8 out of 40): 16 February, 24 and 27 June, 6 and 14 July, 29 August, 23 October in 2021, and 21 February 2022.The dilution effect (S5P passage is more than 12 h after the LF climax time) is probably the explanation for the 16 February, 6 July, and 29 August 2021 cases, while the high presence of ash/ice could be the cause of the TROPOMI underestimation for the other cases.
Figure 10.VPR ash/ice/SO2 total masses and VCTH time series for all the 2020-2022 Etna LFs.The "Date" labels (x-axis) are in the format year-month-day, with letters a and b to distinguish the LFs that occurred on the same day.
The VPR total mass results have also been compared with the lava deposits on the ground obtained by Ganci et al. [61] exploiting the same SEVIRI images (Figure A2).The two different products show quite good agreement, with a Pearson's correlation coefficient between the lava volume and the logarithm of (Ma + Mi + MSO2) of 0.57.

Products Cross-Comparison
Figure 11 shows the cross-comparison between SEVIRI, INGV-VONA bulletin, and TROPOMI VCTHs.Generally, TROPOMI heights are lower than those obtained from SEVIRI and surveillance cameras.However, there is a fundamental difference: the VCTH-DP and VONA estimates are related to aerosols emissions near the volcano (at the beginning of the eruption), while TROPOMI height refers to SO2 which corresponds to a rather well-dispersed plume at quite some distance from Etna.Dilution of the cloud is certainly important in this case.Moreover, the effect of ash on the measurement sensitivity of TRO-POMI can lead to biases in the retrieved SO2 heights.Overall and considering all these different characteristics, the agreement of the various trends is quite good.Quantifying the differences between SEVIRI-VCTH and TROPOMI height, we find large discrepancies (>3 km) only for a few LFs (8 out of 40): 16 February, 24 and 27 June, 6 and 14 July, 29 August, 23 October in 2021, and 21 February 2022.The dilution effect (S5P passage is more than 12 h after the LF climax time) is probably the explanation for the 16 February, 6 July and 29 August 2021 cases, while the high presence of ash/ice could be the cause of the TROPOMI underestimation for the other cases.Figure 12 shows the scatter plots between VPR and LUT p total masses for all the 59 LFs considered.The two procedures appear in good agreement for ash and ice retrievals, while for SO 2 the differences are significant in some cases.SO 2 LUT p results are generally lower than VPR, especially for the bigger activities, when high quantities of ash/ice are present in the volcanic cloud.In these cases, VPR SO 2 retrievals are characterized by bigger uncertainties, as described in Section 2.3, thus probably give out an overestimated value for the SO 2 total mass.However, all the VPR and LUT p ash/ice values agree (within the errors), while regarding SO 2 only one (10 March 2021) is outside the limits.R 2 values are 0.94, 0.93, and 0.14 for ash/ice/SO 2 , respectively, but this latter becomes 0.79 just neglecting the nine LFs with ice total mass greater than 100 Gg.The complete list of VPR and LUT p total masses is given in Table A1.
present in the volcanic cloud.In these cases, VPR SO2 retrievals are characterized by bigger uncertainties, as described in Section 2.3, thus probably give out an overestimated value for the SO2 total mass.However, all the VPR and LUTp ash/ice values agree (within the errors), while regarding SO2 only one (10 March 2021) is outside the limits.R 2 values are 0.94, 0.93, and 0.14 for ash/ice/SO2, respectively, but this latter becomes 0.79 just neglecting the nine LFs with ice total mass greater than 100 Gg.The complete list of VPR and LUTp total masses is given in Table A1.In Figure 13, SO2 total masses from SEVIRI and TROPOMI are shown.The temporal correspondence between the various images was obtained considering a maximum time difference of +/−1 h (but usually is less than 10 min), resulting in only 22 LFs.On the secondary axis (right), the sum of ash and ice total mass is reported also to highlight the presence or not of large quantities of particles inside the volcanic cloud.Generally, the agreement is quite good for weak eruptions while it is poor for bigger eruptions where SEVIRI SO2 results are found to be greater than TROPOMI.These paroxysms were characterized by a high particle (ash/ice) content (dashed line) that certainly can affect all the measurements.
Figure 13.SO2 total masses from SEVIRI-VPR, SEVIRI-LUTp, and TROPOMI.The error bars represent the retrieval uncertainties of the three procedures.On the secondary axis (right), the sum of VPR ash and ice total mass is reported also as a reference.The "Date" labels (x-axis) are in the format year-month-day.In Figure 13, SO 2 total masses from SEVIRI and TROPOMI are shown.The temporal correspondence between the various images was obtained considering a maximum time difference of +/−1 h (but usually is less than 10 min), resulting in only 22 LFs.On the secondary axis (right), the sum of ash and ice total mass is reported also to highlight the presence or not of large quantities of particles inside the volcanic cloud.Generally, the agreement is quite good for weak eruptions while it is poor for bigger eruptions where SEVIRI SO 2 results are found to be greater than TROPOMI.These paroxysms were characterized by a high particle (ash/ice) content (dashed line) that certainly can affect all the measurements.present in the volcanic cloud.In these cases, VPR SO2 retrievals are characterized by bigger uncertainties, as described in Section 2.3, thus probably give out an overestimated value for the SO2 total mass.However, all the VPR and LUTp ash/ice values agree (within the errors), while regarding SO2 only one (10 March 2021) is outside the limits.R 2 values are 0.94, 0.93, and 0.14 for ash/ice/SO2, respectively, but this latter becomes 0.79 just neglecting the nine LFs with ice total mass greater than 100 Gg.The complete list of VPR and LUTp total masses is given in Table A1.In Figure 13, SO2 total masses from SEVIRI and TROPOMI are shown.The temporal correspondence between the various images was obtained considering a maximum time difference of +/−1 h (but usually is less than 10 min), resulting in only 22 LFs.On the secondary axis (right), the sum of ash and ice total mass is reported also to highlight the presence or not of large quantities of particles inside the volcanic cloud.Generally, the agreement is quite good for weak eruptions while it is poor for bigger eruptions where SEVIRI SO2 results are found to be greater than TROPOMI.These paroxysms were characterized by a high particle (ash/ice) content (dashed line) that certainly can affect all the measurements.
Figure 13.SO2 total masses from SEVIRI-VPR, SEVIRI-LUTp, and TROPOMI.The error bars represent the retrieval uncertainties of the three procedures.On the secondary axis (right), the sum of VPR ash and ice total mass is reported also as a reference.The "Date" labels (x-axis) are in the format year-month-day.Figure 14 shows the comparison between SO 2 column densities derived from SEVIRI (both with VPR and LUT p procedures) and from TROPOMI measurements.Two LFs were considered: the strong eruption of 23 October 2021 in which SEVIRI SO 2 total mass is much greater that from TROPOMI (first row) and the weak one of 9 August 2021 in which the two sensors almost agree (second row).In these two examples, the spatial extent of the volcanic clouds is almost the same, while in the case of 23 October 2021, the SO 2 column densities are significantly higher for SEVIRI than for TROPOMI, especially in the central part of the volcanic cloud.The reason for that is the presence of a high concentration of ash and ice particles that cause an overestimation and an underestimation of the columnar SO 2 content when retrieved in the TIR and UV spectral ranges, respectively [21,22].
greater that from TROPOMI (first row) and the weak one of 9 August 2021 in which the two sensors almost agree (second row).In these two examples, the spatial extent of the volcanic clouds is almost the same, while in the case of 23 October 2021, the SO2 column densities are significantly higher for SEVIRI than for TROPOMI, especially in the central part of the volcanic cloud.The reason for that is the presence of a high concentration of ash and ice particles that cause an overestimation and an underestimation of the columnar SO2 content when retrieved in the TIR and UV spectral ranges, respectively [21,22].

Conclusions
In this work, the characterization of the 66 paroxysmal explosive episodes that occurred at Etna volcano in the period 2020-2022 has been presented.This analysis was performed by processing the satellite images collected from the SEVIRI sensor on board MSG geostationary satellite every 15 min, to retrieve the volcanic cloud top height (VCTH), the ash, ice, and SO2 total masses' time series as well as the cumulative masses for the whole period.
The VCTHs were computed by exploiting three different algorithms ("Darkest Pixel", "Cloud Tracking", and "HYSPLIT").To compute the most reliable VCTH value for each event, a weighted average procedure has been developed considering the uncertainties and limits of applicability of the different methods.The values obtained vary between 4 and 13 km, and the entire period can be approximately divided into three main time ranges characterized by average VCTH values of 9 (from 13 December 2020 to 19 March 2021), 6 (from 24 March 2021 to 17 June 2021) and 10 km (from 19 June 2021 to 21 February 2022).A very similar trend (Pearson's correlation coefficient of 0.65) is found in Calvari and Nunnari [24] which processed the measurements collected from the ground-based TIR cameras to retrieve the maximum LF heights.The results obtained were compared with the maximum volcanic cloud top altitude reported in the VONA bulletin and with

Conclusions
In this work, the characterization of the 66 paroxysmal explosive episodes that occurred at Etna volcano in the period 2020-2022 has been presented.This analysis was performed by processing the satellite images collected from the SEVIRI sensor on board MSG geostationary satellite every 15 min, to retrieve the volcanic cloud top height (VCTH), the ash, ice, and SO 2 total masses' time series as well as the cumulative masses for the whole period.
The VCTHs were computed by exploiting three different algorithms ("Darkest Pixel", "Cloud Tracking", and "HYSPLIT").To compute the most reliable VCTH value for each event, a weighted average procedure has been developed considering the uncertainties and limits of applicability of the different methods.The values obtained vary between 4 and 13 km, and the entire period can be approximately divided into three main time ranges characterized by average VCTH values of 9 (from 13 December 2020 to 19 March 2021), 6 (from 24 March 2021 to 17 June 2021) and 10 km (from 19 June 2021 to 21 February 2022).A very similar trend (Pearson's correlation coefficient of 0.65) is found in Calvari and Nunnari [24] which processed the measurements collected from the ground-based TIR cameras to retrieve the maximum LF heights.The results obtained were compared with the maximum volcanic cloud top altitude reported in the VONA bulletin and with SO 2 center of mass heights obtained from TROPOMI data.Despite the very different methods, the results show the same trend and are in good agreement.
The volcanic cloud detection and the discrimination between ash and ice in the volcanic cloud have been realized by exploiting the BTD procedure.The thresholds needed for the discrimination are computed, for each LF, by using RTM simulations performed with MODTRAN 5.3 and considering the SEVIRI spectral characteristics, the VCTH, and the specific atmospheric conditions.These values range from 0.48 to 1.83, with a mean value of 1.04, which can be used as a first approximation for Etna to avoid RTM computations.
More than two thousand SEVIRI images were processed and, for each of them, the corresponding ash/ice effective radius, aerosol optical depth, columnar abundance, and SO 2 columnar abundance maps were obtained.Trends of ash/ice/SO 2 total mass are in good correlation with VCTH and in most cases, the ice content is greater than ash.Due to the generally very short, but intense, nature that characterized these LFs, the maximum values of the total masses can be considered a good estimate of the quantities of particles and SO 2 emitted into the atmosphere by the volcano, as the dilution effect is negligible.
It is also interesting to note that the logarithm of the ash/ice/SO 2 total masses' trends are in good agreement (Pearson's correlation coefficient of 0.57) with the lava deposits obtained by processing the same SEVIRI images [61].
A good agreement is also found between SEVIRI VPR and LUT p ash and ice retrievals, while VPR results are significantly bigger than LUT p for SO 2 when high ice content is present in the volcanic cloud.
Table 1 reports the cumulative and the average ash/ice/SO 2 total masses considering all the LF episodes as also the three main periods identified from VCTHs retrievals.It is important to remark that the characteristics of most Etna 2020-2022 LFs (short episodes but sometimes very energetic) make the quantitative estimation of the various parameters critical.The formation of volcanic clouds rich in ash and ice (opaque pixels) makes the retrieval errors greater than in other circumstances.Critical is the cross-comparison with S5P-TROPOMI due to the different spectral characteristics of the two sensors (TIR vs. UV).In this case, the correlation is good only when ash and ice have low value contents.As expected, when high ash and ice values are present, SEVIRI SO 2 retrievals are greater than those found with TROPOMI.The retrieval of ice particles produced in the volcanic cloud is another valuable result of this work.For almost all the LFs, the ice production is relevant and the mass is generally higher than the mass of the ash particles.From the ice mass and knowing the amount of atmospheric water vapor at the height of the cloud, it should be possible to compute the amount of water vapor emitted during the eruption.These results could be used to try to understand, in detail, the mechanisms through which both atmospheric and volcanic water vapor, together with ash particles, contribute to the formation of ice particles in the volcanic cloud.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/rs15082055/s1, Figure S1: The complete ash/ice/SO 2 time series plots for all the 2020-2022 Etna Lava Fountains using the VPR procedure.

Figure 3 .
Figure 3. Example of VCTH computation using the Cloud Tracking (CT) procedure applied on the 16 June 2021 Etna lava fountain: (a) computation of mean wind speed (b) comparison with an atmospheric profile collected in the same region at the closest time.

Figure 3 .
Figure 3. Example of VCTH computation using the Cloud Tracking (CT) procedure applied on the 16 June 2021 Etna lava fountain: (a) computation of mean wind speed (b) comparison with an atmospheric profile collected in the same region at the closest time.

Figure 4 .
Figure 4. Example of VCTH computation using the HY procedure applied on the 9 August 2021 Etna lava fountain.SEVIRI image collected at 14:00 UTC, RGB composition (R: BT 12.0 −BT 10.8 ; G: BT 10.8 −BT 8.7 ; B: BT 10.8 ) with the volcanic cloud outlined in yellow.The output of the HYSPLIT model considering three forward trajectories at different altitudes (8, 9, and 10 km asl) is overlaid on the figure.
(4.5-12.5 km), and the View Zenith Angle (VZA) whose typical values for SEVIRI "Central Mediterranean Area" are 35°-55° plotted as error bars.Specifically, 12 monthly average (years 1981-2010) atmospheric profiles (PTH) and sea surface temperatures (SST) computed from NCEP/NCAR database for the South Italy area (32.5°-42.5°N;10.0°-20.0°E)have been considered and used in MODTRAN 5.3 simulations.The BTD dependence on VZA (not shown in figure) is quite linear and directly proportional: higher VZA produce higher BTD thresholds and vice versa.
shows the trend of the BTD thresholds (obtained as just described) as a function of the months of the year, VCTH (4.5-12.5 km), and the View Zenith Angle (VZA) whose typical values for SEVIRI "Central Mediterranean Area" are 35°-55° plotted as error bars.Specifically, 12 monthly average (years 1981-2010) atmospheric profiles (PTH) and sea surface temperatures (SST) computed from NCEP/NCAR database for the South Italy area (32.5°-42.5°N;10.0°-20.0°E)have been considered and used in MODTRAN 5.3 simulations.The BTD dependence on VZA (not shown in figure) is quite linear and directly proportional: higher VZA produce higher BTD thresholds and vice versa.

Figure 6 .
Figure 6.(a) Inverted arches curves of ice (cyan lines) and ash (black lines) obtained for MSG4-SEVIRI 18 February 2021 Etna LF using MODTRAN 5.3 RTM and varying AOD0.55 and Re.The BTD Figure 6.(a) Inverted arches curves of ice (cyan lines) and ash (black lines) obtained for MSG4-SEVIRI 18 February 2021 Etna LF using MODTRAN 5.3 RTM and varying AOD 0.55 and R e .The BTD threshold (BTD = 0.8) is obtained with a horizontal linear fit (red line) of the points with maximum R e .(b) Theoretical MSG4-SEVIRI BTD thresholds monthly variation as a function of VCTH of Etna and the VZA (identified by error bars).

3 )
Figure7shows the VCTH for all the 2020-2022 Etna LFs using the three different procedures described in Section 2.1 (DP, CT, and HY).Each method has strengths and weaknesses depending on the conditions, so to get the "best estimate", a weighted average is presented (see Equations (1)-(3)).The choice of the weights derives from the intrinsic characteristics of the three methods and from our experience: DP is considered the most reliable for high altitudes (strong eruptions with completely opaque pixels), while HY turns out better for low altitudes (weak eruptions).Because the CT method and the HYSPLIT model use a single atmospheric profile and a complete set of meteorological data, respectively, the latter is considered generally more reliable.The overall VCTH uncertainty was obtained from the combination of the uncertainties (according to Equations (1)-(3)) of the three methods used, resulting in a mean value of +/−0.5 km.VCTH = 0.5•VCTH DP + 0.2•VCTH CT + 0.3•VCTH HY for VCTH DP > 9 km (1)

Figure 7 .
Figure 7. Time series of the VCTH obtained using the DP, CT, and HY procedures for all the 2020-2022 Etna LFs.The error bars represent the uncertainties of the three methods as described in Section 2.1.The red solid line indicates the VCTH weighted mean (Equations (1)-(3)) used as input in the VPR procedure.The "Date" labels (x-axis) are in the format year-month-day, with letters a and b to distinguish the LFs that occurred on the same day.Generally, the three methods agree well with each other: the mean absolute difference is 1.1, 1.2, and 1.4 km for DP vs. CT, DP vs. HY, and CT vs. HY, respectively.The VCTH ranges from 4 km (18 January 2021) to 13 km (23 October 2021).As the figure shows, the results emphasize three main periods (highlighted by the background colors of the plot): after an initial phase with low altitudes volcanic cloud (except for 21 December 2020), the period from 16 February to 19 March 2021 was characterized by high volcanic

Figure 7 .
Figure 7. Time series of the VCTH obtained using the DP, CT, and HY procedures for all the 2020-2022 Etna LFs.The error bars represent the uncertainties of the three methods as described in Section 2.1.The red solid line indicates the VCTH weighted mean (Equations (1)-(3)) used as input in the VPR procedure.The "Date" labels (x-axis) are in the format year-month-day, with letters a and b to distinguish the LFs that occurred on the same day.

Figure 10 .
Figure10.VPR ash/ice/SO 2 total masses and VCTH time series for all the 2020-2022 Etna LFs.The "Date" labels (x-axis) are in the format year-month-day, with letters a and b to distinguish the LFs that occurred on the same day.

Figure 11 .
Figure 11.Cross-comparison between the VCTH obtained from SEVIRI as the weighted mean of the DP, CT, and HY procedures, TROPOMI SO2 altitude, and VONA reports.The error bars represent the height uncertainties of the three sensors.The "Date" labels (x-axis) are in the format year-monthday, with letters a and b to distinguish the LFs that occurred on the same day.

Figure 11 .
Figure 11.Cross-comparison between the VCTH obtained from SEVIRI as the weighted mean of the DP, CT, and HY procedures, TROPOMI SO 2 altitude, and VONA reports.The error bars represent the height uncertainties of the three sensors.The "Date" labels (x-axis) are in the format year-month-day, with letters a and b to distinguish the LFs that occurred on the same day.

Figure 13 .
Figure13.SO 2 total masses from SEVIRI-VPR, SEVIRI-LUT p , and TROPOMI.The error bars represent the retrieval uncertainties of the three procedures.On the secondary axis (right), the sum of VPR ash and ice total mass is reported as a reference.The "Date" labels (x-axis) are in the format year-month-day.

Figure 14 .
Figure 14.SO2 column densities [g m −2 ] from SEVIRI-VPR (a,d), SEVIRI-LUTp, (b,e), and TROPOMI (c,f) for two LF episodes (at the top the 23 October 2021, at the bottom the 9 August 2021).Note that the color ramp for the two LFs is not the same.

Figure 14 .
Figure 14.SO 2 column densities [g m −2 ] from SEVIRI-VPR (a,d), SEVIRI-LUT p , (b,e), and TROPOMI (c,f) for two LF episodes (at the top the 23 October 2021, at the bottom the 9 August 2021).Note that the color ramp for the two LFs is not the same.

Figure A1 .
Figure A1.Comparison between the SEVIRI-VCTH and the maximum LF height obtained from TIR ground-based cameras by Calvari and Nunnari [24].

Figure A2 .
Figure A2.Comparison between VPR Ash/Ice/SO2 total masses with the ground lava deposits obtained by Ganci et al.[61] from SEVIRI images.

Figure A1 .
Figure A1.Comparison between the SEVIRI-VCTH and the maximum LF height obtained from TIR ground-based cameras by Calvari and Nunnari [24].

Figure A2 .
Figure A2.Comparison between VPR Ash/Ice/SO2 total masses with the ground lava deposits obtained by Ganci et al.[61] from SEVIRI images.

Figure A2 .
Figure A2.Comparison between VPR Ash/Ice/SO 2 total masses with the ground lava deposits obtained by Ganci et al.[61] from SEVIRI images.

Table 1 .
Cumulative and average (in brackets) ash/ice/SO 2 total masses in the overall period (OP) December 2020-February 2022 (upper row) and the three specific periods identified (P1, P2, and P3), both from VPR and LUT p procedures.