Volcanic Emissions, Plume Dispersion, and Downwind Radiative Impacts Following Mount Etna Series of Eruptions of February 21–26, 2021

During the extended activity of Mount Etna volcano in February–April 2021, three distinct paroxysmal events took place from February 21 to 26, which were associated with a very uncommon transport of the injected upper‐tropospheric plumes toward the north. Using a synergy of observations and modeling, we characterized the emissions and three‐dimensional dispersion for these three plumes, monitored their downwind distribution and optical properties, and estimated their radiative impacts at selected locations. With a satellite‐based source inversion, we estimate the emitted sulfur dioxide (SO2) mass at an integrated value of 55 kt and plumes injections at up to 12 km altitudes, which qualifies this series as an extreme event for Mount Etna. Then, we combine Lagrangian dispersion modeling, initialized with measured temporally resolved SO2 emission fluxes and altitudes, with satellite observations to track the dispersion of the three individual plumes. The transport toward the north allowed the height‐resolved downwind monitoring of the plumes at selected observatories in France, Italy, and Israel, using LiDARs and photometric aerosol observations. Volcanic‐specific aerosol optical depths (AODs) in the visible spectral range ranging from about 0.004 to 0.03 and local daily average shortwave radiative forcing (RF) ranging from about −0.2 to −1.2 W m−2 (at the top of atmosphere) and from about −0.2 to −3.0 W m−2 (at the surface) are found. The composition (possible presence of ash), AOD, and RF of the plume have a large inter‐plume and intra‐plume variability and thus depend strongly on the position of the sampled section of the plumes.


Introduction
Volcanic activity, spanning from passive degassing to explosive eruptions, can have relevant downwind impacts on the atmospheric composition (e.g., von Glasow et al., 2009), aerosol properties (e.g., Sellitto et al., 2017), air quality and health (e.g., Michaud et al., 2005), the formation and lifetime of clouds (e.g., Malavelle et al., 2017) and the radiative balance (e.g., Andersson et al., 2015;Kloss et al., 2021;Sellitto et al., 2020). Volcanoes emit both primary and secondary pollutants. The main primary volcanic emissions are gases such as water vapor, carbon dioxide, and sulfur dioxide (SO 2 ), and particles such as ash. Small (submicron) secondary sulfate aerosols (SAs) can form from the gaseous-phase and liquid-phase oxidation/nucleation of primary SO 2 emissions. All these emissions can have environmental and climatic impacts downwind of the dispersion pathways of the volcanic plumes. The SA, in particular, are the main contributors of the volcanogenic modulation of the Earth's radiative balance and can impact the local, regional, or global climate system (Sellitto & Briole, 2015). This is due to their high reflectivity properties in the solar spectral range, mirrored by their very high single scattering albedo (SSA) (Krotkov et al., 1997), and the relatively long atmospheric residence time linked to their tiny typical sizes (Stevenson et al., 2003). While the importance of moderate-to-strong stratospheric volcanic eruptions on the Earth's climate system is now relatively well established (Ridley et al., 2014;Santer et al., 2014), the role of smaller eruptions, with injection in the troposphere (e.g., Mather et al., 2003), on the regional climate system, and the local impacts of persistent passive degassing activity (Sellitto et al., 2020), is not yet well understood and quantified (Oppenheimer et al., 2011).
The downwind impacts of tropospheric volcanic plumes and their spatio-temporal extent depend critically on the interplay of three main aspects: (a) the internal geochemical and geophysical processes linked to magma degassing and plume injection, (b) its atmospheric physico-chemical evolution, and (c) the local and regional atmospheric dynamics driving the plume dispersion. To characterize these three components of plumes release and evolution, and to estimate their downwind impacts, a coordinated synergy of ground-based and satellite observations, and the modeling of transport and physical-chemical evolution, are necessary. Such synergies have been exploited in a number of studies in the past (e.g., Haywood et al., 2010;Kloss et al., 2021;Sellitto et al., 2016;Webley et al., 2012).
In this paper, we couple ground-based and satellite observations with dispersion and impact modeling to study a specific phase of the extended activity of Mount Etna occurring during the period February-April 2021. We focus, in particular, on the series of paroxysmal eruptions during the period from February 21 to 26, 2021. This period was characterized by an uncommon local dynamics that allowed a rare transport of the emitted plumes toward the north, where a number of ground-based observatories could be used to monitor the morphological, optical, and radiative properties of the plumes. Vertically resolved information can be retrieved from stations equipped with LiDAR (Light Detection And Ranging) systems. From different points of view, this is a unique case-study to gather better insights into the evolution of a moderate volcanic eruption. For this series of events, the volcanic activity lasted for a long time, the intensity of the individual eruptions was unusually strong, the atmospheric injection was particularly high and the plumes were transported toward areas with a relatively large density of downwind observatories. We show that a synergistic observational/modeling approach provides a complete end-to-end characterization of the volcanic events, from emissions to dispersion to downwind impacts.
The paper is organized as follows. In Section 2, we introduce and describe the array of observations and modeling tools used in this work. In Section 3, we provide a qualitative description of the eruptive series under investigation, and we frame this volcanic activity into the wider Mount Etna's internal geophysical context. In Section 4, we show results for the characterization of the emission, dispersion, and impacts of this series of events. We draw conclusions in Section 5.

Volcanic SO 2 Emission Flux Rates and Total Mass Estimations From Satellite-Based Source Inversion
The SO 2 emission flux rates time series and total mass have been obtained by exploiting the measurements collected from the Spinning Enhanced Visible and InfraRed Imager (SEVIRI). The SEVIRI instrument is a multispectral imager on board the Meteosat Second Generation (MSG) geostationary satellite. It has 12 spectral channels from visible to Thermal InfraRed (TIR) spectral ranges and a spatial resolution of 3 km at subsatellite 10.1029/2021JD035974 3 of 21 point. The temporal resolution varies from 5 (rapid scan mode over Europe and Northern Africa) to 15 min (Earth full disk) (see https://www.eumetsat.int/meteosat-second-generation for more details on MSG platform and SEVIRI instrument). The SEVIRI measurements are collected in real time from the Multimission Acquisition SysTem (MAST) developed at Istituto Nazionale di Geofisica e Vulcanologia (INGV) . All the SEVIRI images used in this work have been resampled in a regular grid of 3 × 3 km 2 and processed every 15 min.
The SO 2 columnar abundance is computed by applying the Volcanic Plume Retrieval (VPR) procedure (Guerrieri et al., 2015;Pugnaghi et al., 2013Pugnaghi et al., , 2016. This approach is based on the computations of the plume transmittances in the SEVIRI TIR bands centered at 8.7, 10.8, and 12 μm, obtained through linear equations based on the original image radiances and a new image with volcanic plume signal removed; the latter is computed by means of a linear regression of the radiance values outside the edges of the plume itself. The SO 2 estimates are obtained from 8.7 μm transmittances, after eliminating the contribution of ash and ice particles computed from 10.8 and 12 μm channels. The main advantages of VPR method are that the only required input is the volcanic cloud top height and that it can be applied to every geostationary or polar multispectral TIR sensor and for every volcano in the world.
From the SO 2 abundance, the SO 2 emission flux rate is obtained by applying the "traverse approach" (Merucci et al., 2011) considering a transect placed at 30 ± 1.5 km from the summit craters and a wind speed derived from the ARPA (Agenzia Regionale per la Protezione Ambientale) database (Scollo et al., 2009). Such a distance for the reference transect has been selected to minimize the retrieval uncertainties induced by both the opacity of the pixels too close to the craters and the dilution of the pixels too far from the emission. Knowing the wind speed, the SO 2 flux at 30 km is then reported to 0 km over the vents (Corradini et al., , 2021. From SO 2 emission rate, the SO 2 total emitted mass is obtained by temporal integration. The volcanic cloud top height is obtained by applying the consolidated "dark pixel" procedure (Corradini et al., 2018), based on the comparison between the minimum SEVIRI 10.8 μm brightness temperature of a pixel contained in a fixed area over the summit craters and an atmospheric temperature profile measured approximately in the same area and at the same time of satellite acquisition (Corradini et al., 2018). In this case, an area of 19 × 19 SEVIRI pixels and ARPA database were used, respectively. This method is based on the assumption that the plume is opaque and in thermal equilibrium with its environment.

Plume Dispersion With the FLEXPART Lagrangian Model
We simulate the dispersion of the SO 2 plumes emanated by the Mount Etna activity during the period February 21-26, 2021, using the Lagrangian dispersion model FLEXPART (Pisso et al., 2019). The simulations are initialized with the emission rates and altitude estimated with SEVIRI, using the methods described in Section 2.1. A vertical emission profile has been modeled based on the parameterization of Mastin et al. (2009), as done by Sellitto et al. (2016). As meteorological inputs, our FLEXPART simulations use European Centre for Medium-Range Weather Forecasts (ECMWF) ERA-5 reanalysis data at 30 km horizontal resolution and 137 height levels (100-200 m vertical resolution in the troposphere, 1-2 km in the lower stratosphere). The FLEX-PART outputs are averaged over 30 min intervals and are given at 0.1° × 0.1° horizontal grid and 9 altitude levels, from the surface to 13 km altitude, with a vertical resolution of 1 km, from 5 to 13 km, and a unique vertically broad layer, from surface to 5 km. Dry and wet deposition, as well as chemical depletion, are considered in the simulations, using the default FLEXPART sink parameterization for SO 2 .

Plume Radiative Forcing Estimations With the UVSPEC Radiative Transfer Model
We estimate the radiative impact of the Mount Etna aerosol plumes for the events under investigation using the radiative transfer model (RTM) UVSPEC and the LibRadtran package (Emde et al., 2016). The equinox-equivalent clear-sky daily average shortwave (integrated between 300 and 3,000 nm) surface and top of the atmosphere (TOA) direct radiative forcings (RFs) are estimated with a similar methodology as done by Sellitto et al. (2016Sellitto et al. ( , 2020; please refer to these previous works for more details on the setup of the radiative model. Only the aerosol component of the volcanic plumes under investigation is considered for these radiative studies. The radiative perturbation associated with the volcanic SO 2 and water vapor, and the gaseous sulfuric acid formed in the plumes, can be significant for major stratospheric eruptions like El Chichon 1982 and Hunga Tonga 2022 (Gerstell et al., 1995;Sellitto et al., 2022). These impacts are considered relatively small for the weaker eruption of the present study and are thus neglected here.
The background atmospheric state is set using the AFGL (Air Force Geophysics Laboratory) winter mid-latitudes climatological standard. The surface albedo is set to 0.15, independent from wavelength, which is a typical average shortwave albedo for the vegetated surfaces underlying the locations associated with our study. Baseline and volcanically perturbed radiative transfer calculations are carried out using different aerosol layers configurations: using aerosol extinction profiles at downwind LiDAR stations, as volcanically perturbed conditions, and with the same aerosol extinction profiles with the identified volcanic plume layer removed, as nonvolcanic baseline (see specific configurations in Section 4.4). Different hypotheses must be considered for the non-measured optical parameters of volcanic aerosols. The spectral variability of the volcanic aerosol extinction, the SSA, and the angular distribution of the scattered radiation are not directly measured and can be quite uncertain due to the possible presence of both SA and a fraction of residual ash. Values of the Ångström exponent (AE, a compact measure of the spectral variability of the particle extinction and linked to the average size of the particles in the plume), of the SSA (linked to the absorption properties of the aerosol layer) and of the asymmetry parameter (g, linked to the angular distribution of the aerosol layer) are selected for each case described in Section 4.4; please refer to this section for the details of these different hypotheses. The AE, SSA, and g values in each experiment are taken as wavelength-independent. Different runs are used, for all experiments, for different solar zenith angles (SZAs) to sample the full day. The daily average shortwave TOA RF for the volcanically perturbed aerosol layer is calculated as the SZA-averaged upward diffuse irradiance for the baseline simulation minus the volcanically perturbed simulation, integrated over the whole shortwave spectral range. The shortwave surface RF is calculated as the SZA-average net (upward minus downward) global (direct plus diffuse) irradiance for the baseline simulation minus the volcanically perturbed simulation, integrated over the whole spectral range.

SO 2 , SA, and Dust Observations With IASI
The Infrared Atmospheric Sounder Interferometer (IASI) is a Fourier transform spectrometer covering the large infrared spectral range between 645 and 2,760 cm −1 , with a relatively high spectral resolution (0.5 cm −1 ) and a relatively small radiometric noise (noise equivalent spectral radiance of about 20 mW/(cm −1 m 2 sr) around 1,000 cm −1 ). The IASI instrument observes the Earth's atmosphere and surface with a downward-looking geometry, with circular footprints of down to 12 km radius at the nadir spaced by 25 km, and with a swath of 2,200 km (Clerbaux et al., 2009). The instrument series is on board MetOp-A, -B, and -C spacecrafts since 2006, 2012, and 2018, respectively. Each instrument of the series provides a near-global coverage every 12 hr. For this work, we have used IASI observations of SO 2 , SA, and dust.
The IASI daily SO 2 data set includes SO 2 columns and volcanic plume altitude. The SO 2 partial columns at six altitude intervals between 5 and 25 km are retrieved from IASI observations by exploiting the correlation between brightness temperature differences and SO 2 total columns with assumptions on the SO 2 plume altitudes (Clarisse et al., 2012). From these partial columns' profiles, the total column is obtained by vertical integration. In addition, the altitude of the SO 2 plume is estimated by using a sensitive trace gas detection method for high spectral infrared measurements (Clarisse et al., 2014).
The detection of SA and its aerosol-type-specific optical depth is provided by the AEROIASI-H 2 SO 4 retrieval algorithm (Guermazi et al., 2021). The AEROIASI-H 2 SO 4 algorithm is based on a self-adapting Tikhonov-Phillips regularization method, built around the RTM Karlsruhe Optimized and Precise Radiative transfer Algorithm (KOPRA) (Stiller et al., 2000). This method iteratively fits each IASI spectrum at 20 spectral micro-windows by adjusting the vertical profile of SA number concentration, jointly with water vapor profiles and surface temperature. The SAs are assumed as spherical droplets of an aqueous solution of sulfuric acid (H 2 SO 4 ), and the corresponding refractive indices are taken from Biermann et al. (2000). The information on the total SA-specific optical depth is mainly provided by the selective absorption of the undissociated H 2 SO 4 in the SA solution droplets. The interference of the co-existent SO 2 is avoided by using an operational spectral micro-window around 900 cm −1 , a spectral region characterized by H 2 SO 4 absorption bands but outside the main SO 2 -sensitive region around 1,100-1,200 cm −1 .
For the characterization of the three-dimensional distribution of desert dust plumes during the event analyzed in the paper, we use the AEROIASI-Dust approach (Cuesta et al., 2015(Cuesta et al., , 2020. The AEROIASI-Dust algorithm derives vertical profiles of desert dust, in terms of extinction coefficient at 10 μm, from individual thermal infrared spectra measured by IASI satellite sensor in cloud-free conditions, both over land and ocean. This method iteratively fits each IASI spectrum at 12 spectral micro-windows by adjusting the vertical profile of dust aerosol abundance jointly with surface temperature. The information on the vertical distribution of dust is mainly provided by their broadband radiative effect, which includes aerosol thermal emission depending at each altitude on the vertical profile of temperature. As for AEROIASI-H 2 SO 4 , the AEROIASI-Dust algorithm uses the RTM Karlsruhe Optimized and Precise Radiative transfer Algorithm (KOPRA) (Stiller et al., 2000). The current paper uses AEROIASI retrievals from the version described by Cuesta et al. (2020) but using a complex refractive index derived by Di Biagio et al. (2017) with a dust sample from a desert in Mali (which shows a smaller difference between IASI and simulated spectra than that used by Cuesta et al. (2020)).

SO 2 and AI Observations With Sentinel-5p TROPOMI
The TROPOspheric Monitoring Instrument (TROPOMI), onboard the Copernicus Sentinel 5 Precursor (S5p) spacecraft, is a joint project of the Netherlands Space Office and the European Space Agency (Veefkind et al., 2012). The TROPOMI is a satellite-based spectrometer operating in the ultraviolet/visible/near-infrared spectral region and at a nadir-viewing geometry. It provides observation of key atmospheric constituents at the unprecedented high spatial resolution of down to 7 × 3.5 km 2 . In the present work, we use the offline Level 2 data sets of Aerosol Absorbing Index and the total SO 2 column.
The UltraViolet Aerosol Absorbing Index (UVAI) is a compact parameter used to visualize the presence of absorbing aerosols in the sampled air masses. It is based on the spectral ratio of the measured TOA reflectance and a pre-calculated theoretical reflectance for a Rayleigh scattering-only atmosphere, at a given pair of UV wavelengths. Positive residuals of the observed and modeled reflectances are linked to the presence of UV-absorbing aerosols, like dust, smoke, and dense ash layers. Negative residuals may indicate the presence of non-absorbing aerosols, while values close to 0 are found in the presence of clouds or no aerosols. As the UVAI is dependent on different aerosol layer optical and morphological properties, and on the underlying surface reflectance, in this paper, we use this parameter only as a general indication of the presence of dust.
The total column SO 2 is retrieved using a DOAS (Differential Optical Absorption Spectroscopy) method. The SO 2 slant column density (SCD), that is, the SO 2 concentration along the mean light path through the atmosphere, is derived by fitting cross-sections in a spectral range characterized by the absorption of SO 2 (in this case, at different spectral fitting micro-window between 310 and 390 nm, depending on the SO 2 burden). After different calibrations and corrections of the derived SO 2 SCD, this latter is converted to vertical columns using the air mass factor from radiative transfer calculations. For more details about the TROPOMI total column SO 2 product, please refer to Theys et al. (2017).

SO 2 Observations With OMPS-NM
The Ozone Mapping and Profiler Suite Nadir Mapper (OMPS-NM) flies on the Suomi National Polar-orbiting Partnership (Suomi-NPP) satellite since 2012 and measures the Earth backscattered UV/visible radiation at a nadir-viewing geometry (Dittman et al., 2002). In this work, the Level 2 height-resolved SO 2 product is used. For this product, a direct vertical column fitting algorithm is used to retrieve the SO 2 column amounts in the lower (centered at 2.5 km), middle (centered at 7.5 km), and upper (centered at 11 km) troposphere, as well as the lower stratosphere (centered at 16 km).

Aerosol Classification With CALIOP
The Cloud-Aerosol LIdar with Orthogonal Polarization (CALIOP) is a LiDAR system onboard the Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observation (CALIPSO) spacecraft; it flies at about 700 km altitude in a sun-synchronous orbit and was part of the so-called A-Train (Winker et al., 2010) before being transferred in September 2018 in the so-called C-Train (see https://svs.gsfc.nasa.gov/31049). The CALIOP space LiDAR is in orbit since April 2006 and is still in operation. It provides vertical profiles of aerosols and clouds at about 01:30 and 13:30 local time. In this work, we use aerosol detections of the Vertical Feature Mask level 2 product version 4.2, at varying horizontal resolution, as well as the related aerosol types classification information.

MPLNET LiDAR Observations at Sede Boker Site
The National Aeronautics and Space Administration (NASA) Micropulse Lidar Network (MPLNET; Welton et al., 2001) project was established in 1999 in support of the NASA Earth Observing System (EOS; Wielicki et al., 1995). During the last two decades, MPLNET observations have significantly contributed to fundamental studies and applications on climate change and air quality studies and in support for NASA satellite and suborbital missions (Bilal et al., 2019;Campbell et al., 2021;Lolli et al., , 2020. This federated global network is constituted by homogeneous commercially available Micropulse LiDARs (MPL) instruments manufactured by Droplet Measurement Technology. The MPLNET instruments are single-wavelength LiDAR systems that use a diode-pumped Nd:Yag laser at 532 nm, with a 1-min temporal resolution and 75-m vertical resolution. The MPLNET network is active during the whole day and night and operates in all weather conditions, within the limit of laser signal extinction. The MPLNET network is deployed at global scale and covers most latitudes bands. For this study, the observations from Sede Boker (30.8°N, 34.8°E, 480 m a.s.l.) permanent MPLNET site are specifically used. To retrieve the vertically resolved optical aerosol and cloud properties from a single wavelength LiDAR, strong assumptions, that limit retrieval accuracy, are needed. To reduce this uncertainty, MPLNET LiDAR network deploys the instruments, if possible, together with an Aerosol Robotic Network (AERONET; Holben et al., 1998) sun-photometer to constrain the LiDAR equation and reduce the retrieval errors (Welton et al., 2000). More details on the full MPLNET standard data product suite, including aerosols, have been provided by Welton et al. (2018). The standard, automated, MPLNET products are not designed to accurately capture elevated volcanic plumes due to their transient nature, high altitude, and unique microphysical properties relative to the local aerosols that limit the usefulness of the constrained aerosol retrieval approach. Instead, MPLNET provides a custom retrieval method that bypasses the automated aerosol height algorithm, and requires manual specification of the plume altitude (for this study based on inspection of the lidar signal and volume depolarization ratio profiles). Upon specification of the plume base and top, the optical depth of the plume is determined using the attenuation of the lidar signal and a molecular scattering (Lewis et al., 2016). The retrieval of plume properties is calculated using the same MPLNET-constrained retrieval algorithm using the plume optical depth. This process often results in few retrievals during the plume's advection over the lidar site as calculation of the plume optical depth is difficult with noisy signals. Instead, the Lidar Ratios from these few retrievals are averaged, and the result is used to determine the plume properties by fixing the Lidar Ratio using the classic Fernald lidar algorithm (Fernald, 1984). The retrieved lidar ratios were in good agreement with several LiDAR studies on volcanic emissions (e.g., Prata et al., 2017) and a fixed value of 50 sr was used.

The ACTRIS-Fr OHP LiDAR Observatory
The Observatoire de Haute-Provence (OHP) located in southern France (43.9°N, 5.7°E, 670 m a.s.l.) is one of the NDACC (Network for the Detection of Atmospheric Composition Change) alpine stations, equipped with a variety of LiDAR instruments for monitoring of the lower and middle atmosphere.
The longest continuous LiDAR data record at OHP is provided by a Rayleigh-Mie-Raman LiDAR for temperature and aerosol measurements (hereafter referred to as LTA) operating at the wavelength of 532 nm. The LTA instrument (Keckhut et al., 1993) has provided routine measurements for over 3 decades with a mean measurement rate of 10-12 nighttime acquisition nights per month and a typical duration of acquisition of 3-5 hr. To retrieve the vertical profiles of aerosol backscatter and extinction, the raw lidar returns from three LTA elastic channels, covering respectively troposphere, lower stratosphere, and upper stratosphere, are merged, and a Fernald-Klett inversion method is applied (Fernald, 1984;Klett, 1985), assuming a constant LiDAR ratio of 50 sr. The retrieved backscatter and extinction coefficients are reported with a time resolution of 1 min and a vertical resolution of 15 m. The scattering ratio is then computed as a ratio of total (molecular plus aerosol) to molecular backscattering, where the latter is derived from ECMWF meteorological data. A more detailed description of the instrument, aerosol retrieval and error budget is provided by Khaykin et al. (2017) and references therein. The ultraviolet tropospheric differential absorption LiDAR (LiO 3 tr) also operated at OHP at the same time as the LTA LiDAR. It records backscatter signal at 316 nm for ozone monitoring (Ancellet & Beekmann, 1997;Gaudel et al., 2015). The 316 nm signal is very weakly absorbed by ozone in the troposphere and can be used for aerosol backscatter monitoring between 3 and 13 km altitude. As for the LTA instrument, the Fernald-Klett inversion method (Fernald, 1984) is applied assuming a constant LiDAR ratio of 50 sr at 316 nm and a 1.01 scattering ratio close to a molecular return at 11 km. Since no overlap correction is applied to match the LTA retrieval at 3 km, a scattering ratio uncertainty of the order of 10% is expected below 4 km.
The OHP observatory is also equipped with an automatic LiDAR (CIMEL model CE376) sounding in the troposphere from ∼200 m to ∼12 km of altitude at the wavelength of 532 nm, named GAIA and operational since 2019. The GAIA system is also equipped with a depolarization channel at the same wavelength to gather more information on the nature and optical properties of aerosols, mostly in the boundary layer and in the free troposphere.
It provides data at a frequency of 1 min and the datasets are then averaged on 10 min slices. Its spatial resolution is 15 m. The scattering ratio is computed as for the LTA LiDAR, but interpolating radiosoundings collected at OHP launched every week in the late morning. The data are corrected for instrumental deadtime, overlap in the first hundreds of meters of altitude (using an algorithm called CHECK delivered by the PHOTONS-AERONET network), sky background, and squared altitude to take into account the quadratic decrease of the signal as a function of altitude. A Haar wavelet method (Brooks, 2003) can then be used to retrieve the boundary layer height. The GAIA data are inverted using the Klett method (Klett, 1981).

The ACTRIS Napoli LiDAR Observatory
The Napoli LiDAR observatory (Naples, Italy, 4.18°E, 40.84°N, 118 m a.s.l.) is part of the ACTRIS (Aerosol, Clouds and Trace Gases Research Infrastructure) research infrastructure, a pan-European distributed research infrastructure for short-lived atmospheric constituents producing high-quality data in the area of atmospheric science (Wandinger et al., 2020). The Napoli observatory includes passive and active remote sensing systems and near-surface sampling systems, for atmospheric studies. Thanks to its quite central location in the Mediterranean basin, it represents a strategic location to study optical and microphysical properties of the aerosols coming from local sources and long-range transport.
The station is equipped with a multiwavelength elastic/Raman LiDAR device able to detect elastic signals at 355, 532, and 1,064 nm, Raman N 2 echoes at 386 and 607 nm, and aerosol depolarization at 532 nm. Details on this LiDAR system are reported by Boselli et al. (2021). Vertical profiles of the aerosol extinction coefficient and linear depolarization ratio are shown and discussed in this paper. The aerosol extinction coefficient profiles were retrieved using the procedure introduced by Ansmann et al. (1990). The calibrated particles' linear depolarization profiles were obtained from the backscattered light polarized along both perpendicular and parallel directions with respect to the laser beam polarization and following the inversion procedure reported by Biele et al. (2000) and Freudenthaler et al. (2009).

The AERONET Sun-Photometric Network
Total column measurements of the aerosol optical depth (AOD), AE, and SSA at several wavelengths, among other aerosol parameters, are carried out globally with Cimel sun photometers, as part of the AErosol RObotic NETwork (AERONET; see the AERONET website at http://aeronet.gsfc.nasa.gov and the description by Holben et al., 1998, for further details). In this manuscript, we use the AERONET observations for these three aerosol parameters at the selected station of Sede Boker, OHP, and Napoli, which are part of the AERONET network.

Summary of Satellite and Ground-Based Observations
The different satellite and ground-based observations used in this paper are listed in the following Table 1.

The Mount Etna Activity in February and March 2021
Between February and April 2021, Mount Etna experienced an intense eruptive activity from the South-East Crater consisting of several astonishing short-lasting lava fountain eruptions coupled with lava effusion and episodic pyroclastic density currents. The activity took place after a period of intense eruptive activity at Bocca Nuova, Voragine, and North-East summit craters which started on September 2019. In contrast to the other three summit craters, over this period, the South-East crater was episodically erupting but featured a slow and gradual energy increase of its eruptive activity, which eventually intensified between December 2020 and January 2021, and abruptly escalated on the afternoon of February 16, 2021, producing the aforementioned sequence of lava fountains. Over these 2 months (February-April 2021), 17 paroxysmal episodes were produced with a recurrent characterization, consisting of the resumption of strombolian activity growing in intensity, lava flow effusion, and transition from discrete explosions to the continuous gas and lava jetting typical of lava fountaining. Overall, each lava fountain lasted a few hours with eruptive columns up to at least 10 km altitude and the atmospheric injection of copious burdens of volcanic ash and SO 2 . At the local scale, this activity produced fall-out deposits in the surrounding towns and disrupted air traffic with Catania airport closed on different occasions. Lava flows propagated mostly in the eastern flank of the volcano in the Valle del Bove and remained confined at a mean altitude of 2,700 m asl. The volcanological framework of this extended activity of Mount Etna volcano has been discussed by Bonaccorso et al. (2021), Calvari and Nunnari (2022), and Corsaro and Miraglia (2022). In the present study, we specifically focus on the three paroxysmal lava fountain events that occurred between February 21 and 26, 2021.
During the same time interval of this study, a major Saharan dust outbreak to central Europe occurred (Text S1 and Figures S1 and S2 in Supporting Information S1). In the region of interest of our study, the dust stays confined between 3 and 6 km altitude. The volcanic and dust plumes remain completely vertically separated in this area, thus facilitating the detection and characterization of the volcanic plumes using vertically resolved observations.

SO 2 Emission Rates for the Three Paroxysms
We estimate the SO 2 emission rate and plume injection altitude for the three individual paroxysms of February 21-26, 2021, with the satellite-based method described in Section 2.1. In Figure 1, the temporal profile of the emission rate and altitude, for Events #1, #2, and #3, are shown. The main characteristics for the three events are summarized in Table 2. All times in this section and the whole paper are reported as UTC (Universal Time Coordinated). We associate with these three events the plumes P1, P2, and P3, respectively (see Section 4.2). Events #1 and #3 have a similar duration (∼4-5 hr) and peak SO 2 emission rate (∼2,000-2,500 kg/s). Event #2, though shorter (∼2 hr), has the largest peak SO 2 among the three events (∼5,000 kg/s). For the three events, the injection reached altitudes as high as ∼12 km.
The total emitted SO 2 mass for the three events is estimated at 19.3, 20.9, and 14.6 kt for Events #1, #2, and #3, respectively. Event #2 injects the largest amount of SO 2 , despite the shortest duration among the three events, due to the significantly larger peak emission rate. The total SO 2 mass emitted from the three added events is 55 kt (0.055 Tg).
Typical values of the emission rates for moderate-to-strong eruption of Mount Etna are generally less than 2000 kg s −1 Sellitto et al., 2016). Then, this specific activity is characterized by relatively brief but very intense paroxysmal events, if compared to more general behavior of Mount Etna. In addition, the altitude of the plume injection is particularly high. Documented extremely high-altitude eruptions of Mount Etna rarely exceed 10-11 km altitude (e.g., Corradini et al., 2020;Sellitto et al., 2016). Thus, this series of three paroxysms qualify as an extreme event in terms of the injection altitude. The total SO 2 mass emitted during these few days of activity is about 50% of the record-breaking sequence of eruptions of February 24-31, 2015, which was estimated, with a similar method as the present paper, at about 100.0 kt (0.1 Tg) overall . Thus, in terms of the SO 2 emission rates, injection altitude, and total SO 2 emitted mass, the sequence of three events during February 21-26, 2021 is in the higher part of the spectrum of eruption strengths of Mount Etna.

The SO 2 Plumes Dispersion Toward the North
The plumes emanating from the three individual paroxysmal events described in the previous section are then characterized in terms of their subsequent dispersion.  Figure S3 in Supporting Information S1. As extensively discussed by Sellitto et al. (2016), the synergy of high-quality regional observations from satellites and realistic transport simulations is crucial toward getting more comprehensive insight into confined plume's dispersion. These two information layers are complementary. The vertical sensitivity of satellite observations of SO 2 is limited, thus satellite observations can only provide a partial information on the height distribution of SO 2 . On the other hand, the accuracy of dispersion simulations is limited by the uncertainties in the input parameters and in the description of physico-chemical in-plume processes. In this work, the dispersion simulations are initialized with the SO 2 emission rates and plume injection altitudes observations described in Section 4.1. In addition, the temporal sampling of satellite observations, one or two overpasses per day for low-Earth orbit spacecraft, is usually not sufficient to smoothly characterize the temporal evolution of the plume's dispersion. On the contrary, dispersion simulations can be realized at subhour temporal resolution, thus filling the temporal gap between satellite overpasses. An animation of the FLEXPART simulations for this case study, with the full temporal resolution of 30 min, can be found in Supplementary video (Data Set S1 in Supporting Information S1). Figures 2-4 show that simulations of the plume dispersion are consistent with observations, thus cementing our confidence in our further interpretation of the plume's dynamics and downwind impacts.
Plume P1 (Event #1) disperses toward the north (February 21) and then northeast (February 22), at altitudes decreasing from 10 to 12 km (February 21) to 8-10 km (February 22). SO 2 columns as large as 15 DU are found from both observations and modeling. Then, the plume flows around an anticyclonic circulation due to the presence of a stable Omega block in Central Europe (Hoshyaripour, 2021), dispersing toward the south-eastern direction (February 23) and then further east, reaching the Eastern Mediterranean (February 24). The SO 2 columns remain relatively large (>10 DU). Plume P2 (Event #2) began on February 23, with a similar initial northern dispersion as P1.
The plume progressively spreads in the east-west direction and forms an elongated structure across southern and central Europe. Plume P2 was much more irregular than plume P1, in terms of both SO 2 column distributions and plume altitude. Values as high as 15 DU are found at the western and eastern ends of the plume, with mean altitudes at about 8-12 km. Smaller column values, around 5-6 DU, and a mean altitude of about 10-12 km are found in the central area of the plume. Two CALIOP intercepts of the early P2 plume show altitudes as large as 12.5 km ( Figures S2 and S5 in Supporting Information S1). Then, the plume structure compresses (February 25) and disperses toward the north-east (February 26). During this phase, the eastern end of plume P2 sweep through southern Italy at 8-9 km altitude. Plume P3 (Event #3) appears on February 25 morning and disperses toward the north-east, quickly joining the dispersion of plume P2. It displays slightly larger values of the SO 2 column (>15-18 DU) and lower altitude (7-10 km) than plumes P1 and P2. This overall description of the dispersion of the plumes P1-3 in the upper-troposphere (UT) is confirmed by the independent observations of the upper-tropospheric SO 2 column from OMPS-NM ( Figure S4 in Supporting Information S1).
The dispersion of plumes P1-3 was consistently directed toward the northern quadrant from the perspective of the Mount Etna volcanic source. Using a more-than-one-decade-long dispersion simulations data set, Sellitto et al. (2017) have estimated the prevalent direction of dispersion of the plumes generated from Mount Etna emissions. They have found that for more than 80% of the time, Mount Etna's plumes disperse toward the eastern quadrant. The dispersion toward the northern and north-western quadrant, as is the case for these events, is a rare occurrence, accounting only for less than 5% of the cases. On the other hand, active continental long-term observatories, that can allow the downwind observation and characterization of Mount Etna's plumes, are only present in the northern and north-western sectors. Thus, this series of paroxysms is a rare occasion to study more in details the dispersed plume of Mount Etna's volcanic emissions.

Volcanic Plumes Observations at Downwind Observatories
Plumes P1 and P2 have been observed downwind at different observatories. In the following, we analyze the plume's observations at three stations: Sede Boker, Israel, OHP, France, and Napoli, Italy. Figure 5a shows the volcanic plume vertical distribution over Sede Boker station, in terms of the SO 2 volcanic tracer, from FLEXPART simulations. Figure 5b shows the series of LiDAR aerosol observations at Sede Boker, for the same period. Plume P1, after a relatively long dispersion phase over central-easter Europe, driven by anticyclonic flow of the Omega block discussed in Section 4.2, overpasses the Sede Boker station starting from February 25 early in the morning and, with discontinuous broken dispersed plumes sections, until February 26 in the early afternoon (Figure 5c). More dense sections of the plume P1 overpass the station during the period spanning from February 25 night to February 26 early morning. During this plume's overpass, P1 has aged 4-5 days since their volcanic emissions. The LiDAR observations show a vertically confined aerosol layer in the same period indicated by the FLEXPART simulations. Even if a part of the plume observation is missed by the LiDAR, during the night between February 24 and 25 (due to low LiDAR signal quality), a consistent vertical behavior is shown by simulations and observations. Both information layers show an initially higher-altitude plume, located at 9-10 km, then descending to 5-6 km altitude. The plume looks filamentary, with a geometrical vertical thickness of less than 500 m (Figures 5a and 5b) and a horizontally confined shape (Figure 5c). The mean aerosol depolarization for the plume is 32.3 ± 1.4% (Table 3). This value is very large for volcanic aerosols and might indicate the presence of a large fraction of fine ash. The three-dimensional morphology of this aerosol layer and the consistency of observations with volcanic tracer dispersion simulations, suggest that a significant presence of dust can be excluded at these altitudes.
The western and eastern peripheries of the plume P2 (hereafter referred to as P2/W and P2/E, respectively) are almost simultaneously observed at OHP (P2/W) and Napoli (P2/E) ground stations. Starting from February 24, the plume P2 takes an elongated and then folded-shaped geometry, with the western periphery P2/W slowly moving south-westwards over Southern France and Eastern Spain, and the eastern periphery P2/E sweeping Southern Italy in a north-western direction (Figures 3 and 4 and Figure S2 in Supporting Information S1).
According to FLEXPART simulations, the P2/W section arrives at OHP station in two separate phases: during the evening/night of February 24 and during the night/early morning between February 25 and 26 (Figure 6c). The  2. IASI SO 2 column observations for the morning (approximate overpass at ∼9:30 LT, panels (a, e)) and afternoon (approximate overpass at ∼21:30 LT, panels (b, f)) MetOp-A spacecraft overpass; corresponding FLEXPART SO 2 column simulations (at 10:00 LT, panels (c, g), and at 22:00 LT, panels (d, h)), for February 21 and 22, 2021. Plume P1 is individuated and its approximate direction of dispersion is indicated as a blue arrow in panels (a, b, e, and f). In panels (c,d, g, and h), black crosses individuate the position of OHP and Napoli observatories.   FLEXPART simulations over OHP show a denser volcanic plume, in terms of the SO 2 tracer, at altitudes of about 9-10 km altitude, with less dense volcanic plume at lower altitudes (2-4 km). This is generally consistent with the observations of the LTA LiDAR system at OHP (Figures 6d and 6e). The LTA LiDAR was operational from about 18:30 to 22:30 on February 24 and February 25. It observes an aerosol layer at about 10-11 km altitude on both days and during the full period, with a maximum value of the aerosol extinction between 19:00 and 21:00 on February 24, consistently with FLEXPART simulations. The very weak volcanic aerosol signature at 10-11 km for LTA LiDAR observations on February 25 is not present in FLEXPART simulations. A previous observation session was carried out during the late afternoon/evening of February 23 (not shown here); a volcanic-related aerosol plume is not observed during this period, consistently with FLEXPART time-series analyses. While the LTA LiDAR system was operational only during nighttime, the co-located GAIA system observed the atmosphere over OHP during daytime as well (data not shown here), even if with a significantly lesser signal/noise ratio at those altitudes. GAIA observations confirm the nighttime volcanic plume overpass. In addition, GAIA observations did not detect any volcanic plume overpass during most of daytime on February 24 and 25, thus corroborating the description given by the FLEXPART simulations. Besides the relatively high-altitude aerosol layer observations, the LTA LiDAR observed a lower-altitude aerosol plume at 3-5 km (Figure 6a). This aerosol   signal can be easily associated with the dust event described in Text S1 in Supporting Information S1, even if a component of volcanic aerosol cannot be excluded based on the FLEXPART vertical profiles of Figure 6c. For a limited time, the LTA and the LiO 3 t LiDARs operated simultaneously at OHP. Using a combination of the aerosol backscattering profiles retrieved by LTA (at 532 nm) and LiO 3 t (at 316 nm) LiDARs, a color ratio has been calculated in the free troposphere on 24/02 between about 19:00 and 20:30. The aerosol backscattering profile observations at both wavelengths show a vertically thin plume at about 10-11 km and a wider plume between 3 and 5 km. Besides this quite distinct vertical structure, the color ratio shows very different values for the two layers: 0.1-0.2, for the highest first layer, and up to 1.0-1.2 for the lowest layer. In addition, depolarization ratio observations with the GAIA system for the highest layer display near-zero values. These observations indicate the presence of distinctly smaller mean particle size, with a possible spherical shape, for the highest layer with respect to the lowest, thus corroborating the possibility that volcanic aerosols, with an SA prevalence, and dust, respectively, dominate the two layers. Vertically separated layers of dust and volcanic aerosols are observed in the area with CALIOP, at altitudes consistent with the previous discussion ( Figure S2 in Supporting Information S1). A confined dust layer is also observed with IASI, through the AEROIASI-Dust retrieval algorithm, between 3 and 5 km ( Figure S1 in Supporting Information S1). Spectroscopic evidence of the presence of SA at OHP is found using the AEROIASI-H 2 SO 4 retrieval algorithm ( Figure S6 in Supporting Information S1). This suggests that the volcanic plume P2/W is dominated by SA, even if the presence of an ash fraction cannot be excluded.
The P2/E section is observed at Napoli station. Based on the FLEXPART simulations, the P2/E volcanic plume overpasses Napoli during the period spanning from late morning to early evening of February 25, at an altitude of 7-10 km, with a maximum density between 15:00 to 18:00 ( Figure 7b). Consistently, the LiDAR system in Napoli observed a vertically confined aerosol plume in the afternoon of February 25, with a maximum at around 17:30 and 7-8 km altitude (Figure 7a). A second aerosol plume, with a wider vertical structure, is observed between 2 and 3 km. The LiDAR system in Napoli operates at two wavelengths and has depolarization information (depolarization ratio δ). Using the two wavelengths, the bulk Ångström exponent (AE) has been calculated. Both AE and δ show distinct values for the two layers: AE of 0.85 and 0.24, and δ of 10.1 ± 5.8% and 26.6 ± 1.7%, for the higher and lower layers, respectively. This indicates significantly smaller and more spherical-shaped particles for the higher layer, thus corroborating the identification of the highest layer as volcanic and the lowest as dust. The AE and δ for the highest layer are nevertheless quite different from what expected for a pure SA-dominated layer, indicating the possibility of the presence of a fraction of larger and aspherical ash particles. Using LiDAR-based aerosol size distributions inversion methods, two different plumes are found at Napoli station, the higher dominated by volcanic aerosols and the lower by dust (Sannino et al., 2022). It is also important to notice that: (a) the SO 2 plume tracer of P2/E in Napoli ( Figure 7b) has peak values nearly 2 orders of magnitude larger than for P2/W at OHP (Figure 6c), and (b) the aerosol extinction is almost 10 times larger for P2/E (Figure 7a) than P2/W (Figure 6d). This clearly indicates the large variability of the morphology and properties of one volcanic plume, for example, when looking at different plume's sections. Based on the results shown in this section, we can conclude that: (a) precise (e.g., in terms of the input parameters) and detailed dispersion simulations, as well as complementary satellite observations at the regional scale, and (b) synoptic downwind observations at different locations, can be used to characterize in details confined plumes of known source and their variability in terms of three-dimensional morphology, composition, and properties. In the present case, for example, two separate and non-mixing vertical layers of volcanic aerosols and dust could be identified. For the volcanic aerosol layer, the variability of the vertical structure and composition at different sections of the plumes (see the case of plumes P2/W and P2/E), as well as the different properties of individual plumes for a given geophysical event (see the comparison of P1 and P2), could be characterized.

Optical Properties and the Impact on the Radiative Balance
We estimate the optical and radiative properties of plumes P1 and P2/E and W using the observations and atmospheric simulations of the previous sections.
The volcanic-plume-specific aerosol optical depth (AOD V ), Ångström exponent (AE V ), and depolarization ratio (δ V ) are listed in Table 3. These are obtained by averaging the LiDAR observations at Sede Boker, IHP, and Napoli observatories during the core plume overpass periods (see details in the caption of Table 3). These values are also compared with the total column (i.e., volcanic plus dust and other possible aerosols in the line of sight of the instrument) AOD T , AE T , and SSA T estimated with co-located AERONET Cimel sun photometers.
For the plume P1, moderate AOD T is found, with relatively large AE T , which points at a limited impact of dust in the column. The AOD V is relatively large, 0.0076 ± 0.0015. A large δ V (32.3 ± 1.4%), with a relatively large AOD V , may point at a consistent presence of ash within the plume P1. This is somewhat surprising for a plume that, at the time of downwind observation, has already aged about 4-5 days since emissions. It is worth noticing that this plume shows a consistent descent during the dispersion (emissions at about 12 km and detection at Sede Boker at 5-10 km altitude) which can be linked to the larger gravitational settling velocity of the massive ash particles. Nevertheless, it is important to note that the plume descent is also modeled for the SO 2 tracer in FLEX-PART simulations. Thus, the relatively fast plume descent seems dominated by regional dynamics, with possible additional gravitational settling for larger ash particles.
The plume P2 is very inhomogeneous, when comparing the western (P2/W) and easter (P2/E) ends. This is mirrored by inhomogeneities in the measured AOD T . The plume P2/W has almost 10 times smaller AOD T with respect to P2/E (P2/W: 0.0036 ± 0.0014, P2/E: 0.031 ± 0.006). As discussed in Section 4.3, this large difference in the AOD T of P2/W and E can be attributed to the fact that P2/E is largely dominated by SA while P2/W likely presents a fraction of ash. With AOD T between 0.25 and 0.30, the total column is largely dominated by dust (Text S1 in Supporting Information S1), with the volcanic component taking only a small part of the column. Relatively small AE T seems to confirm this hypothesis through the indication of the dominance of very large particles (i.e., larger than a few μm) in the overall vertical column.
Using the LibRadtran-UVSPEC RTM, driven by aerosol extinction profiles measured by LiDARs at the three stations, the daily average RF of the isolated volcanic plumes has been estimated with the methodology described in Sect. 2.2.2, and is reported in Table 4. The aerosol extinction profiles used as input to the RTM are the average of the observations at the core overpass time, that is, the same period used to calculate the AOD V in Table 3 and detailed in its caption. Based on the considerations above, we have made the following hypotheses for nonmeasured optical parameters of the plumes. For plume P1 at Sede Boker, a fraction of fine ash was probably present in the volcanic plume. Thus, we have run the RTM using an interval of AE typical of small-sized particles (1.8-2.0, three runs with 0.1 AE increase) and an SSA typical of partially absorbing particles (two runs with 0.90 and 0.95 SSA). A moderately dusty background lower-tropospheric atmosphere is considered for altitudes between surface to 3 km. For plume P2/W at OHP, an SA-dominated plume is considered. Thus, we have run the RTM using an interval of AE typical of small-sized particles (1.8-2.0, three runs with 0.1 AE increase) and an SSA typical of very reflective aerosols (two runs with 0.95 and 1.00 SSA) (Table 4). A very dusty background lower-tropospheric atmosphere is considered for altitudes between surface to 5 km, based on LiDAR observations of the lowest layer of, for example, Figures 6a and 6b. For plume P2/E at Napoli, a fraction of ash in the volcanic plume is considered. Thus, we have run the RTM using the observed AE (0.85), which is typical of larger particles than what found for P1 and P2/W, and SSA typical of partly absorbing aerosols (two runs with 0.90 and 0.95 SSA). In a similar manner as for plume P2/W, a very dusty background lower-tropospheric atmosphere is considered for altitudes between surface to 5 km, based on LiDAR observations of the lowest layer of for example, Figure 7a. With these assumptions, a TOA RF of −0.22 ± 0.06, −0.17 ± 0.03, and −1.19 ± 0.29 W m − 2 , and a surface RF of −0.56 ± 0.08, −0.18 ± 0.04, and −2.99 ± 0.40 W m −2 are found for P1, P2/W, and P2/E, respectively. Associated with a larger AOD V , the RF of P2/E is significantly larger than for P1 and P2/W. The surface RF is much larger than the TOA RF, for both P1 and P2/E, due to the presence of absorbing aerosols in these plumes. This means that part of the incoming shortwave radiation does not reach the surface due to the in-plume absorption, which adds to the radiation scattered back to space. For plume P2/W, the TOA and surface RF have a very similar value, which points at a very limited absorption within the plume due to the likely dominant presence of purely reflective SA.
The RF values obtained for these plumes can be compared with RFs estimated in the past for proximal and distal observations of Mount Etna's emissions. Using ground LiDAR observations as input to the RTM, Sellitto et al. (2020) have estimated the RF of a typical passive degassing plume at a proximal location, about 7 km downwind Mount Etna craters and have found values of about −4.5 and −7.0 W m −2 , at TOA and surface. Mild explosions and a small fraction of ash were present during these observations. For an almost purely SA, after a dispersion of about 350 km downwind Mount Etna, Sellitto et al. (2016) have found an RF efficiency, that is, the TOA RF per AOD unit, of about −40 to −50 W/m 2 /AOD. The RF efficiency of P2/E (about 38 W/m 2 /AOD) and P1 (about 29 W/m 2 /AOD) are consistent with estimations of Sellitto et al. (2020) for an ash-containing passive degassing plume (about 35 W/m 2 /AOD). The RF efficiency of P2/W (about 47 W/m 2 /AOD) is consistent with estimations of Sellitto et al. (2016) for purely SA plume. Compared with larger volcanic eruptions, our estimates are of the same order of magnitude or even larger than the impact of the Raikoke eruption in 2019 (global average RF: −0.3 to −0.4 W m −2 ) (Kloss et al., 2021). The Raikoke had the largest volcanic global-average TOA RF documented in the last 30 years. Of course, Raikoke eruption is much more relevant than the activity of Mount Etna discussed in the present study, due to its larger scale (hemispheric) impact. This demonstrates that, even if Mount Etna eruptions are of a limited global climate interest, they can have a quite relevant impact at the regional scale.

Conclusions
The synergy of volcanic plume's observation and modeling at different spatio-temporal scales is crucial to characterize their emission, dispersion, and downwind impacts. Mount Etna experienced an intense eruptive activity between February and April 2021. During this phase, three peculiar extreme events, in terms of the large SO 2 emission rates (55 kt of SO 2 emitted overall) and injection altitude (up to 12 km), took place between February 21 and 26, that we have studied with a combination of regional-scale satellite observations and Lagrangian dispersion modeling, and local downwind measurement at selected ground stations coupled with offline RTM modeling. These three events and the subsequently formed volcanic plumes (plumes P1-3) displayed a very uncommon dispersion toward the north, which allowed the downwind observation and characterization of the plumes thanks to the presence of different observatory sites in this quadrant. Ground-based LiDAR observations and complementary information from observations and modeling at Sede Boker (plume P1), OHP (western end of plume P2 -P2/W), and Napoli (eastern end of plume P2 -P2/E) revealed a complex inter-plume (P1 vs. P2) and intra-plume (P2/W vs. P2/E) variability of the morphology, composition, optical properties and radiative impacts of Mount Etna plumes. Plume P1 was detected at about 10 then descending to 5-6 km after only 4-5 days atmospheric aging, with depolarization ratio consistent with a possible presence of ash. Plume P2 displayed a dramatic variability from east to west ends, with a thin and likely SA-dominated plume to the west and a dense ash-bearing plume to the east. Both the plumes' AOD (0.004-0.03) and local clear-sky daily average shortwave RF (−0.2 to −1.2 W m −2 , at TOA and −0.2 to −3.5 W m −2 at the surface) point at a relevant impact on the upper-tropospheric aerosol layer and the regional climate at the Mediterranean scale. The synergy of observations and modeling presented in this work allowed to empirically disentangle the information about volcanic aerosols and other aerosol sources, like mineral dust. FLEXPART (FLEXPART, 2022) and LibRadtran (LibRadtran, 2022) are freely downloadable softwares. IASI (AERIS, 2022) and TROPOMI (TROPOMI AAI, 2022; TROPOMI SO2, 2022) SO 2 data are freely available. LiDAR data for Napoli and OHP stations (ACTRIS, 2022) and Sede Boker station (MPLNET, 2022) are freely available. CALIOP data are freely available (LARC).