Investigation of Volcanic Emissions in the Mediterranean: “The Etna–Antikythera Connection”

: Between 30 May and 6 June 2019 a series of new flanks eruptions interested the south-east flanks of Mt. Etna, Italy, forming lava flows and explosive activity that was most intense during the first day of the eruption; as a result, volcanic particles were dispersed towards Greece. Lidar measurements performed at the PANhellenic GEophysical observatory of Antikythera (PANGEA) of the National Observatory of Athens (NOA), in Greece, reveal the presence of particles of volcanic origin above the area the days following the eruption. FLEXible PARTicle dispersion model (FLEXPART) simulations and satellite-based SO 2 observations from the TROPOspheric Monitoring Instrument onboard the Sentinel-5 Precursor (TROPOMI/S5P), confirm the volcanic plume transport from Etna towards PANGEA and possible mixing with co-existing desert dust particles. Lidar and modeled values are in agreement and the derived sulfate mass concentration is approximately 15 μ g/m 3 . This is the first time that Etna volcanic products are monitored at Antikythera station, in Greece with implications for the investigation of their role in the Mediterranean weather and climate.


Introduction
Volcanic eruptions can inject huge amounts of particles and gases into the atmosphere, which can have an important impact on regional and global scale [1][2][3]. Volcanic emissions are composed by volcanic particles generated by magma fragmentation and gases. The most abundant is the sulfur dioxide (SO2) which is oxidized by the hydroxyl radicals (OH) to form sulfate aerosols in the atmosphere [4,5]. Through kinetic-based simulations, high-temperature chemistry quickly forms oxidants, such as OH, HO2, and H2O2, can lead to the production of sulfate within a few seconds after the emission of SO2 begins.
The oxidation of SO2 and hydrogen sulfide (H2S) gases generates submicron sulfate particles with a lifetime of several years if those are injected in the lower stratosphere or in the upper troposphere, and from days to a few weeks if those are injected in the troposphere [1]. In the absence of clouds, the photochemical SO2 conversion to sulfate particles proceeds more rapidly during the day and in the summer (5-10% h −1 SO2 conversion rate) than during the night and in the winter (0.3-1% h −1 ). The rate of conversion depends on the nature of the surface, the presence of co-pollutants, the temperature, and the relative humidity. In addition, if the plume has a high dust/ash density then heterogeneous surface reactions may play a significant role [6]. Hobbs et al. [7] observed oxidation rates in the Mt. St. Helens plume that were comparable to those observed in power-plant plumes (0.1% h −1 ). A minimum of 13 Tg/y of time-averaged SO2 fluxes has been reported by actual volcanic measurements worldwide during the period 1970-1997 [8]. Τhe global SO2 budget from volcanoes has been recently revised to about 23 ± 2 Tg/yr, for the period 2005-2015, based on spaceborne observations from the Ozone Monitoring Instrument (OMI) onboard NASA's Aura satellite [9]; The contribution of volcanoes to the total sulfur emissions in the atmosphere is up to 10% [10]. Remarkably, volcanic emissions also have a bigger impact on the tropospheric aerosol burden than other sulfur sources [11] because volcanoes tend to emit SO2 at higher altitudes than most other surface sulfur emissions, where the lifetime is longer [12].
The impact of volcanic emissions on the planetary radiative budget is further intensified when volcanic eruptions inject particles in the stratosphere. The long residence time in the stratosphere accompanied by the strong interactions between both ash and sulfates with solar radiation increase the Earth's albedo and can cause ozone depletion on a global scale (e.g., at Mt. Pinatubo massive eruption in 1991, record-low ozone abundances were observed over much of the northern hemisphere) [1,13,14]. Earlier research has broadly considered two major aerosol sources in the stratosphere, namely, the direct injection of ash and sulfates by volcanic eruptions and the isentropic transport of carbonyl sulfide [15][16][17].
Additionally, volcanic material injected and transported at flight altitudes similar to dust and smoke [18] poses a serious hazard to aviation since it can cause damage and loss of power of the aircraft engines [19][20][21]. Ash particles can lead to dangerous aircraft engine damage, and products of SO2 have highly corrosive properties [22]. Despite the importance of this process to emergency responders and aviation, the operational systems capable of monitoring tephra dispersal and fallout in near-real-time, and subsequently provide the expected impact assessment are still limited and not fully adapted to the growing requirements of precision and reliability [23].
Forecasting the transport of volcanic plumes for aviation hazards mitigation requires timely analysis of all available observations to initialize and refine ash dispersion forecasts and issue Volcanic Ash Advisories in the wake of a volcanic eruption [24,25]. A variety of satellite sensors which have been used over the past decades, like the TROPOspheric Monitoring Instrument onboard the Sentinel-5 Precursor satellite (TROPOMI/S5P) [26], can detect volcanic plumes. The TROPOMI imaging spectrometer represents a stepchange in gas monitoring from space as it measures in four different spectral regions (UVvisible, near-infrared, shortwave infrared), providing volcanic SO2 vertical column density (VCD) maps at an unprecedented spatial resolution, used to track the horizontal transport of volcanic SO2 clouds [27]. Observations from both satellites and ground-based lidar instruments have been used to monitor and track the altitude and properties of volcanic aerosols [28][29][30].
In the days following a volcanic eruption, SO2 and ash particles are injected in different heights and may take different transport pathways, as a result of the varying vertical and horizontal wind shear [31,32]. Modeling and forecasting the transport and atmospheric concentrations of volcanic aerosol released during volcanic eruptions (e.g., SO2) depend critically on the knowledge of the eruption source term. Unfortunately, the estimation of the source term is difficult to obtain by direct measurements and often relays on ground-based remote systems as video-surveillance systems (e.g., [33]). The main source term of an eruption includes parameters like the release height, the total mass, the mass eruption rate, and the duration of the eruption [34]. An additional parameter describing ash emissions in particular is the particle size distribution that among the input parameters is one of the most difficult parameters to measure in real time [23].
Long-range transport of volcanic aerosols over the Mediterranean is recorded systematically by the European Aerosol Research Lidar Network (EARLINET) (e.g., [35]). Recent observations from the newly established PANhellenic GEophysical observatory of Antikythera (PANGEA) of the National Observatory of Athens (NOA) during the period 30 May-6 June 2019 depict the transport of elevated layers above the Eastern Mediterranean. The presence of these elevated layers is mostly associated with the continuous Etna volcanic activity.
The aim of this study is to investigate the properties of Etna volcanic emissions and their transport paths in the Mediterranean. Mt. Etna is the largest point-source of particulate matter in the atmosphere of the Mediterranean, affecting the atmospheric levels of airborne particles and their deposition rates at both local and regional scales [36][37][38]. For this purpose, we combine satellite observations and modeling tools as well as groundbased remote sensing near the volcano and at the remote island of Antikythera located 765 km downwind of the volcano.
The manuscript is organized as follows: in Section 2 we present the case study of the 30 May-6 June 2019 eruption of Etna, the atmospheric circulation and the transport pathways between Etna and the Antikythera station. In Section 3 we introduce the methodology used in our analysis to describe the long-range transport of Etna emissions. In Section 4 we present the results obtained by applying the methodology. The transport of the volcanic SO2 plume is simulated with the Lagrangian particle dispersion model FLEXPART and it is validated with data from the Modern-Era Retrospective analysis for Research and Applications version 2 (MERRA-2 reanalysis) and remote sensing observations (Polly XT lidar, TROPOMI/S5P, Astrium, The Hague, The Netherlands). Finally, in Section 5 we present our conclusions and indicate directions for future work.

Volcanic Activity/Emissions
The Mt. Etna volcano (37.74° N, 15.00° E, 3300 m above sea level (a.s.l.)) is located in Sicily, Italy, and is one of the most active volcanoes on Earth. This active stratovolcano has historically recorded eruptions over the past 3500 years while it has been erupting periodically since September 2013 to this day. Lava flows, explosive eruptions with ash plumes, and Strombolian lava fountains commonly occur from one or more of its summit craters named Voragine (VOR), Northeast Crater (NEC), Bocca Nuova (BN), and Southeast Crater (SEC). In the latter, recent activity is located on a new cone formed since 2011 on the volcano's eastern flank, namely the New Southeast Crater (NSEC), [39][40][41][42].
On 1 May 2019 the Istituto Nazionale di Geofisica e Vulcanologia, Osservatorio Etneo (INGV-OE) reported Strombolian activity from the BN crater, followed the day after by Strombolian activity at the NSEC. The activity from the summit craters was almost constant, often forming dilute ash emissions that were dispersed quickly by winds.  According to the Volcano Observatory Notice for Aviation (VONA) messages, the volcano observatory issued a warning report, from yellow to orange levels, early in the morning of 30 May 2019 while a red alert was issued from 10:50 UTC to 17:29 UTC when a strong ash emission was observed.
TROPOMI/S5P instrument SO2 retrievals provided the same information ( Figure A1). By this time, lava flows and explosions had produced persistent SO2 plumes that drifted east and north east for over 800 km from the source, as seen by TROPOMI/S5P ( Figure  A1). This volcanic activity did not result in any major air traffic disruption.
The TROPOMI/S5P polar orbiting instrument produces daily global SO2 VCD maps ( Figure A1). It represents the amount of SO2 molecules in a column overhead per unit surface area generally expressed in Dobson Units (1 DU = 2.68 × 10 16 molecules SO2 cm −2 ). From the SO2 VCs, one can easily calculate the total SO2 mass for a given area (e.g., Figure  A1) which is a quantity useful to investigate volcanic activity [43].

Atmospheric Circulation and Transport Pathways
The atmospheric circulation over the eastern Mediterranean is dominated by persistent northerly and westerly winds, favoring the advection of volcanic products from Etna to Greece [44]. In order to demonstrate the strong connection between Mt. Etna volcano air masses and the PANGEA-NOA station HYSPLIT cluster analysis was performed for a five-year period. The HYSPLIT (Hybrid Single-Particle Lagrangian Integrated Trajectory) model was used to compute air parcel trajectories of long-range transport [45] driven by the 6-hourly meteorological dataset Global Data Assimilation System (GDAS) at a resolution of 1° × 1° for a five-year period (2015-2019). During this five-year period, the Mt. Etna volcano has erupted about 50 times according to INGV reports, but this analysis is assessed by the overall contribution of air masses and not only when Etna erupted within this time.
A two-step cluster analysis was applied to the 48-h forward trajectories starting from Mt. Etna volcano to investigate the emissions transport pathways. The first step involves the clustering of 48 h of forward trajectories at six height levels, 3300 m, 5000 m, 7500 m, 10,000 m, 15,000 m, and 18,000 m above ground. In the second step, clustering analysis was applied for the selected period at each of the height levels, to 1816 in total forward trajectories.
The cluster analysis indicates that due to the synoptic circulation characteristics of the Mediterranean, a significant proportion of tropospheric air masses from Etna volcano are transported eastward towards the island of Antikythera, where the PANGEA station is located. For the five years period examined here the percentage of eastward trajectories starting at each HYSPLIT height are 42% for 3.3 km, 43% for 5 km, 65% for 7.5 km, 63% for 10 km and up to 70% for 15 and 18 km. Thus, the predominant transport pattern highlights the strategic location of PANGEA observatory in monitoring volcanic emissions and establishing the so-called here "Etna-Antikythera connection" (Figure 2).

Experimental Design
Our methodology for the analysis of Etna volcanic plume transport for the specific case study of May-June 2019 flank eruption, incorporates the use of a number of numerical modeling systems, namely, the Lagrangian particle dispersion model FLEXPART coupled with the state-of the-art Weather Research and Forecasting model (WRF-ARW), along with MERRA-2 Reanalysis simulations. In addition, Polly XT NOA lidar data from the PANGEA station and satellite observations from the TROPOMI/S5P are utilized.

Modeling
The transport of volcanic ash and SO2 plumes was simulated with the Lagrangian particle dispersion model FLEXPART [3,46,47] in a forward mode. Additionally, for the characterization of air masses and for the observed aerosol layer identification, the FLEXPART-WRF model was used in backward mode for the computation of source-receptor relationships, as well as back-trajectories, with a total of 10,000 particles released at heights 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, and 4.5 km over the Antikythera station. FLEXPART has been used in a large number of similar studies on long-range atmospheric transport including volcanic plumes [3,48,49].
For the description of the event, the dispersion simulations were driven by hourly meteorological fields from the Advanced Research WRF (ARW) model version 4 [50]. The WRF-ARW spatial set up was at 9 × 9 km resolution domain (Europe shown in Figure 4) with 600 × 370 grid points and 33 vertical levels. Simulations were initiated at 00:00 UTC on 30 May 2019 and were completed at 00:00 UTC on 4 June 2019. Table 1 summarizes the Physics Parameterizations (PP) schemes for the WRF-ARW simulations. The initial and boundary conditions for the offline coupled FLEXPART-WRF runs are taken from the European Centre for Medium-Range Weather Forecasts (ECMWF) Integrated Forecast System (IFS) operational analyses with 0.18° × 0.18° resolution and 91 model levels. Sea Surface Temperature (SST) analysis data were provided by the Copernicus Marine Environ-ment Monitoring Service (CMEMS) at a resolution of 1/12°. The use of 1-hourly WRF meteorological fields at a 9 × 9 km spatial resolution allows a more detailed representation of the volcanic plume dispersion.
For SO2, a range of different values for the mass released were used in the model simulations following previous studies in order to select the more realistic one. Coppola et al. [51] incorporating satellite-based observations showed that the average SO2 flux of the Mt. Etna volcanic activity between 2004 and 2010 was equal to 3259 tonnes per day [51]. Granier et al. [52] estimated the average SO2 emissions for 2005-2010 of Mt. Etna volcano to 3456 tonnes per day using data from the NOVAC (Network for Observation of Volcanic and Atmospheric Change) network. In addition, Queißer et al. [27] made a comparison of SO2 fluxes, of one month from TROPOMI/S5P from the automated scanning FLux Automatic MEasurments (FLAME) ground-network, and fluxes from driving traverse measurements underneath the volcanic plume. An excellent agreement was shown for most of the days demonstrating that reliable, nearly real-time, high temporal resolution SO2 flux time series from TROPOMI/S5P measurements are possible for Mt. Etna. The average monthly SO2 flux from TROPOMI/S5P was 2.83 ± 1.66 kt/day while for FLAME 2.39 ± 1.09 kt/day. On 3 June 2019, at ~11:43 UTC, TROPOMI/S5P recorded total SO2 VCD values of approximately 1.17 ± 1.02 DU (about 33.48 ± 29.19 mg/m 2 ), over Antikythera (see Section 3.3 for details). In this study, in order to calibrate the volcanic SO2 mass flux according to the TROPOMI/S5P values measured over Antikythera, we did recursive runs with the FLEXPART-WRF model with the emission rates over Etna ranging between 2 and 5 kt/day. A value of 4 kt/day which is an average of SO2 emissions of Mt. Etna volcano, was finally selected as in this case the model simulations matched the TROPOMI/S5P SO2 plume patterns over the area and the measured SO2 columnar concentrations over the Antikythera station. We present here the results obtained by simulations with SO2 emissions rates equal to 4 kt/day, emitted from the reported start time of the eruption until 4 June 2019, 00:00 UTC, when the volcanic activity on the SE base of Etna's NSEC started to decline in frequency. Land Surface (LSM) NOAH [58] Ash particles simulations were initiated at the reported start time of the eruption on 30 May 2019 and were completed at 17:30 UTC on the same day, at the eruptive stage. Neither TROPOMI/S5P nor other satellite instruments such as the Ozone Monitoring Instrument [59] on board NASA's Aura spacecraft [60] and the Ozone Mapping and Profiling Suite Nadir Mapper (OMPS-NM) on board the NASA-NOAA Suomi National Polarorbiting Partnership (Suomi-NPP) [61] could detect any clear ash cloud. The detection of ash by satellite instruments is more complicated than that of SO2 because of the varying physicochemical properties of the ash particles (different size, geometrical shape, and composition) [62]. The mass eruption rate is usually estimated from the observed injection height with empirical relationships [23,34,63].
In our study we estimate the mass eruption rate (MER) for ash particles following Scollo et al. [23] study, by inverting the observed plume height using the 1-D plume model of Degruyter and Bonadonna [63]. Similar models are commonly used, as they can capture the first-order physics of a volcanic plume rising in the atmosphere, while remaining computationally efficient [64]. The Degruyter and Bonadonna [63] plume model assumes that (i) the plume is in a steady-state, (ii) the solid particle and gas phase in the plume are well mixed, such that they have the same bulk velocity and temperature, (iii) differences in pressure between the plume and the atmosphere are negligible and the velocity and temperature distributions through a cross-section of the plume follow a top-hat profile, which remains self-similar along the plume trajectory [23,63]. The mass eruption rate from plume height was calculated according to Equation (1): where for a weak plume: The initial injection height in the model is set to the surface of the Etna crater (i.e., 3.3 km a.s.l. up to 4 km a.s.l., based on VONA reports). A total of 10,000 particles were released in each model run, for SO2 tracer and another 10,000 particles for ash. For both ash and SO2 simulations dry and wet deposition processes are also enabled. For ash, the gravitational particle settling [65] was determined assuming spherical particles with a density of 3000 kg/m 3 . For SO2, the oxidation by the OH radical was considered as a sink process, similar to earlier studies (e.g., [66]). The size distribution of volcanic ash particles was described using four size bins (5,9,13, and 21 μm diameter) as the particles size distribution relevant for long-range transport refers to the smaller particles (≤25 μm diameter). Considering only the long-range transport of the smallest particles, the total MER calculated by the 1-D plume model of Degruyter and Bonadonna [63] (which consider all the emitted material) is scaled to account for the near-source fallout. Only 5% of the total MER was used for volcanic ash simulations [67][68][69]. Mastin et al. [34] provided estimates of the fraction of emitted mass carried by small ash grains (<63 mm diameter) present in the proximal ash deposits that range from 2 to 60% depending on the type of volcano. Hourly averaged total column values (ash and sulfates mass concentrations, μg/m 3 ) were produced at the output.
Furthermore, the emission of dust is calculated by the Air Force Weather Agency (AFWA) scheme [70]. WRF-Chem includes five dust size bins to represent the evolution of dust, with effective radii of 0.73, 1.4, 2.4, 4.5, and 8 μm.
Modeling synergy such as MERRA-2 reanalysis was also employed in our study in order to further describe the volcanic plume transport from Etna towards Antikythera. The MERRA-2 (Modern-Era Retrospective analysis for Research and Applications, version 2) meteorological and atmospheric composition reanalysis data set used here [71][72][73] is provided by the NASA Global Modeling and Assimilation Office (GMAO) (publicly available at https://disc.gsfc.nasa.gov/datasets/M2T1NXAER_5.12.4/summary). MERRA-2 is produced using the GEOS-5 atmospheric model [74] which integrates a radiatively coupled version of the GOCART model (Goddard Global Ozone Chemistry Aerosol Radiation and Transport) [75] to simulate concentration profiles of five aerosol components, namely, dust, sea salt, black carbon, organic carbon, and sulfate. The horizontal resolution is 0.5° × 0.625° (latitude × longitude) with 72 hybrid-eta model layers in the vertical up to 0.01 hPa. Along with meteorological observations, aerosol AOD at 550 nm is also assimilated in MERRA-2 from AVHRR, MISR, MODIS, and AERONET. Anthropogenic emissions of SO2 are obtained from the EDGAR-4.2 emission inventory, while volcanic emissions of SO2 are derived from the AeroCom Phase II project.

The PANGEA EARLINET Station of Antikythera
The PANGEA observatory of NOA in the remote island of Antikythera, is located across the travel path of different air masses, providing continuous monitoring of essential climate variables in the Eastern Mediterranean. The Aerosol Remote Sensing facility of ACTRIS (Aerosol, Clouds, and Trace Gases Research Infrastructure; actris.eu) is collecting continuous observations of aerosols and clouds since 2018. The observations presented herein, were acquired from 30 May to 4 June 2019, and are used in order to validate the FLEXPART simulation results.
In particular, the parameters used are collected from the CIMEL sun photometer (part of Aerosol Robotic Network; AERONET) and the Polly XT -NOA lidar (part of EAR-LINET). From CIMEL, the variables used are columnar measurements of the aerosol optical depth (τ) and the Angstrom exponent (a), and the particle size distribution. The Polly XT -NOA lidar [76,77] is a multi-wavelength Raman-polarization system with 24/7 operational capabilities, providing vertical distributions of the particle backscatter coefficient (β) at 355, 532, and 1064 nm, the extinction coefficient (α) at 355 and 532 nm and the particle depolarization ratio (δp) at 355 and 532 nm.
With the aforementioned observations, and using well known methodologies, we can separate spherical and non-spherical particles in mixed aerosol layers (e.g., [78]), towards aerosol characterization (e.g., [77]). In addition, using parameterizations based on AERONET retrievals we can estimate vertical distributions of aerosol concentrations for different aerosol species (e.g., dust, smoke, and marine) from 200 m above the ground up to ~16 km (e.g., [79][80][81]). Recently, this approach was also applied in spaceborne lidar observations [82].
The method applied for the separation between volcanic ash and sulfates and the estimation of their concentrations is the so-called "POlarization-LIdar PHOtometer Networking" (POLIPHON) method [83][84][85]. POLIPHON is based on polarization lidar observations and sun-photometer climatology and is a two-step approach to derive the mass concentration profiles of two different aerosol components: a strongly depolarizing component (i.e., volcanic ash, ) and a non-depolarizing component (i.e., sulfates, ). First, the particle linear depolarization ratio (δp); derived by the volume linear depolarization ratio δv and the molecular and particle backscattering profiles [86], is used to separate the non-depolarizing and depolarizing particles contribution to the particle backscatter coefficient (β) for a given wavelength according to Equation (2). Second, the extinction to volume conversion factor derived from AERONET climatology on different aerosol species, are used to convert the backscatter profiles to mass concentration profiles for each component, according to Equation (3): where and stand for the particle linear depolarization ratio of the depolarizing and non-depolarizing component respectively, stands for = depolarizing component or = non-depolarizing component, is the particle mass density in g cm −3 , / is the extinction to volume conversion factor in 10 −12 mm (i.e., the ratio of the aerosol optical depth τ to the column-integrated volume concentration u, both derived from AERONET measurements), and S is the lidar ratio in sr (the ratio of α to β). Values used in the present study are summarized in Table 2. The overall uncertainty of this methodology for the mass concentration of the depolarizing component can reach up to 30-60%, while for the non-depolarizing component slightly larger uncertainties are expected for aerosol layers with pronounced fine-mode particle concentrations (i.e., sulfate layers of volcanic origin) [84,85].

Satellite Observations: TROPOMI/S5P
In this work, satellite-based observations of the SO2 VCD from the state-of-the-art TROPOMI/S5P sensor [26] are used as an additional source of validation for our results (publicly available via the Copernicus S5P hub: https://scihub.copernicus.eu/). TROPOMI with a spatial resolution of 3.5 × 7 km 2 (3.5 × 5.5 km 2 after August 2019) allows for retrievals with 12 (16) times higher resolution than its predecessor Ozone Monitoring Instrument (OMI).
Details about the TROPOMI algorithm can be found in [43] and in [89]. As the retrievals are sensitive to the height of the SO2 plume, an information that is missing, the SO2 VCDs are reported for three different hypothetical SO2 profiles, for a 1 km thick box located at ground level (centered at 0.5 km), centered at 7 km and centered at 15 km a.s.l. Here, to obtain the shape of the SO2 plume emitted from Etna on 3 June 2019 we follow a filtering method based on the study of Theys et al. [90]. First, pixels with cloud fraction smaller than 0.5, with solar zenith angle smaller than 70° and located at central acrosstrack positions (TROPOMI/S5P rows 50-400) are selected. Then, we identify only pixels with slant column density (SCD) larger than three times the uncertainty on the fitted SO2 slant column (SCDE). For each pixel, a box area of 75 × 75 km 2 is defined. The pixel is assumed to be part of the plume only if the total number of pixels with SCD > 3× SCDE within the box area is larger than two and larger than 0.04 times the number of all the pixels (regardless the SCD value). In this work, for demonstration of the volcanic SO2 plume we plot the filtered 7 km product (see Figure 13). However, the TROPOMI/S5P SO2 values reported for the Antikythera Island (Section 3.1) are calculated by averaging all the retrievals within a 15 km radius circle centered in the island for each one of the three products (0-1, 7, and 15 km) and by interpolating linearly the three products at the plume's central height which is simulated by FLEXPART.

Transport of SΟ2 and Volcanic Ash
To understand the main mechanisms driving the long-range transport of volcanic aerosol over the Mediterranean during the event, we examined the major atmospheric processes taking place in the area on a synoptic scale, as resolved by the WRF-ARW model [50]. This transport episode is associated with a closed long wave trough [91][92][93] affecting the North Italy on the isobaric level of 500 hPa. The modeled geopotential height at 500 hPa at 02:00 UTC, 30 May 2019, (which is around the starting time of the eruption) shows a low-height center (at 5600 gpm) located Northeast of Italy (Figure 3a). At 700 hPa the trough, extended from the Tyrrhenian to Ionian Sea (Figure 3b). The corresponding westerlies suggest an eastward transport of material emitted from Etna towards the Eastern Mediterranean and the PANGEA observatory, at altitudes approximately 3 km height a.s.l. (Figure 3b); at this height (~2 to 3 km) the volcanic particles arrived above Antikythera island. During the following four days this upper-level trough with the closed low circulation system at 500 hPa moved eastwards from the Adriatic Sea to the northwest Greece.
FLEXPART was driven with two different source terms for SO2 and ash in order to investigate their atmospheric transport pathways and dispersion. Simulated column concentrations of sulfates and ash are shown in Figure 4  Vertical cross sections over Antikythera on 3 June 2019, at 00:00-02:00 UTC ( Figure 5) reveal that the sulfates plume is located between 1.5 and 3 km height, in the lower troposphere, whereas the ash plume does not reach Antikythera island with the simulated mass concentrations being small (<5 μg/m 3 ). These separations are mainly due to the different emission rates and duration time for SO2 and ash introduced by FLEXPART model in the control run, while the gravitational settling of the large ash particles also contributes to the different transport of sulfates and ash. For volcanic ash mass flux, the knowledge from satellites and ground-based measurements is limited mainly because this eruption was characterized by low-intensity ash emission, whereas for sulfates apart from satellites, in situ measurements and FLEXPART simulations, there is supporting evidence for the existence of the sulphuric plume by stateof-the-art global atmospheric composition products such as MERRA-2 reanalysis. Figure 6 shows simulated mass density (mg/m 2 ) of sulfates originating from Etna as depicted by MERRA-2 reanalysis on 3 June 2019, 00:00 UTC. The overall location of the SO2 cloud stretching across the Ionian Sea seems to be well captured by MERRA-2 confirming the volcanic plume transport from Etna towards Antikythera (see also Figure 4b). Both MERRA-2 reanalysis ( Figure 6) and FLEXPART-WRF sulfate column mass density values (Figure 4b) are approximately 16 mg/m 2 , on 3 June 2019 at 00:00 UTC, above Antikythera island.

Comparison with Ground-Based and Satellite Remote Sensing Observations
The FLEXPART model simulations (Figure 4)   The modeled arrival height of the first plume above Antikythera is centered approximately at 1.5-2.5 km height during 17:00-18:00 UTC on 2 June 2019 (Figures 7a and 10d (magenta line)), whereas the second sulfates plume is centered at 2-3.5 km height from 00:00 to 02:00 UTC on 3 June 2019 (Figures 7a and 11d (magenta line)). In addition, modeled volcanic ash fluxes are simulated near the surface (below 1 km) at very small values on 3 June 2019 after 18:00 UTC (Figure 7b).
Two-time windows of 1 h, and 2 hours and 35 minutes, were selected to calculate the aerosol optical properties using the Raman method [83], from 17:00 to 18:00 UTC and from 00:00 to 02:35 UTC on 2 and 3 June 2019, respectively. These windows were selected based on the FLEXPART simulations showing that both components, ash and sulfates are present over the station during this time. However, these time frames correspond also to a different type of air mass, advected from the Sahara Desert, carrying dust particles. According to WRF-Chem model simulations, these dusty layers extend from close to the surface up to almost 10 km in height, perhaps co-existing with the particles of volcanic nature (Figures 9, 10c and 11c (purple line)). Dust and volcanic ash, being both large non-spherical particles are difficult to distinguish based on lidar observations alone. Both aerosol species produce comparable values of their intensive optical parameters used for aerosol characterization (i.e., lidar ratio (S) and δp, (see [94,95])).
The time-height evolution of Polly XT -NOA lidar observations at PANGEA is seen in Figure 8 in terms of the attenuated backscatter coefficient at 532 nm (top panel), indicative of particles' concentration above the station, and the volume linear depolarization ratio (δv) at 532 nm (bottom panel), indicative of the particles' shape. Non-spherical particles such as dust or volcanic ash produce δv values larger than 10%, while more spherical particles such as marine aerosols or anthropogenic pollution are expected to produce negligible δv values. Taking advantage of this we apply a simple methodology (presented in Section 3.2) to disentangle the contribution of depolarizing and non-depolarizing components of the aerosol layers observed above Antikythera. For 2 June 2019, the lidar derived profiles (Figure 10a,b) suggest that the depolarizing aerosol component being confined between 1.3 and 2 km, with δp values of (~10 to 20%). Whereas, for 3 June 2019, the lidar derived profiles (Figure 11a,b) suggest that the depolarizing aerosol component being confined between 2 and 3 km, with δp values of (~25 to 30%). By comparing the modeled ash ( Figure 7b) and dust (Figures 9, 10c and 11c (purple line)) mass concentrations to these profiles we can see that the agreement with the lidar is more satisfactory for dust particles, being most probably the non-spherical particle component. Modeled ash particles lie below 1 km, at which height lidar δp values are lower than 5% which corresponds to almost spherical particles, most probably sulfates from the volcanic eruption.    Figures 7a and 11d (magenta line)). Moreover, the depth of the sulfate layer seems to be well reproduced by the model, since the modeled concentrations extend at 2-3.5 km for the selected time interval. For the sulfate plume found between 1.5 and 2.5 km on 2 June 2019, at 17:00-18:00 UTC, the sulfate mass concentration values are about 18 μg/m 3 (Figures 7a and 10d (magenta line)), while the corresponding mass from lidar observations is calculated to be 10-25 μg/m 3 (Figure 10d (blue line)) again in good spatio-temporal agreement with the model.
FLEXPART provides similar outputs and conclusions as the lidar measurements as it considers the classification of aerosols with respect to their source regions and age. Specifically, for the observed aerosol layer, the Lagrangian dispersion model FLEXPART was used for a three-day backward simulation. Figure 12a,b show the source-receptor relationships for air masses inside the 0-4 km and 0-500 m a.s.l. layer accordingly, before reaching the study area. The model output is given in terms of the decimal logarithm of the integrated residence time in seconds in a grid box. The most probable aerosol source region and the aerosol type were assigned accordingly. Volcanic particles that were emitted from Mt. Etna along with a second dusty air mass from northern Africa (Tunisia and northern Algeria) merged over the Mediterranean Sea as a result the final air mass that arrives at Antikythera island include a mixture of volcanic, marine and dust particles (Figure 12a) (see also the possible dust-prone areas in in Figure 12b and in Figure A2). In conclusion, the combined information of the backward trajectory analysis at Antikythera station ending at 00:00 UTC, on 3 June 2019 (Figure 12c) and the source-receptor relationships indicate the presence of volcanic and dust particles. In addition, the FLEXPART forward run simulation ( Figure A3) shows the transport path of the sulfates plume which are compared mostly qualitatively to the satellite-based observations of SO2 VCD from the TROPOMI/S5P. The TROPOMI/S5P patterns for the 7 km product on 3 June 2019, at 12:00 UTC, shown in Figure 13 are similar to the simulated ones, the SO2 cloud stretching across the Ionian Sea being well captured by the model. The SO2 vertical column mass density values derived from TROPOMI/S5P instrument on 3 June 2019, at 12:00 UTC above Antikythera is approximately 1.17 ± 1.02 DU which are about 33.48 ± 29.19 mg/m 2 while, the modeled values are also about 25-30 mg/m 2 ( Figure  A3). Figure 13. SO2 vertical column (in DU) from the TROPOMI instrument aboard Sentinel-5P across Etna's plume on 3 June 2019 at 12:00 UTC (see Section 3.3 for details). The RGB background image is from VIIRS aboard Suomi-NPP (NASA worldview; https://worldview.earthdata.nasa.gov/) which has an overpass time very close (less than 5 min) to Sentinel-5P.

Conclusions
Volcanic emissions may remarkably impact the atmospheric composition in regional and global scales. This is the first time that volcanic aerosol layers originating from Etna on 30 May 2019 were observed at the PANGEA Observatory of Antikythera. A synergy of satellite, ground-based and model data for monitoring the transport pathways and quantifying the amount of volcanic emissions reaching Antikythera is used. FLEXPART-WRF sulfate and volcanic ash simulations described the volcanic plume transport from Etna towards Antikythera. Satellite-based SO2 observations from the TROPOMI/S5P instrument and ground-based Polly XT lidar observations in Antikythera, confirmed the modeled simulations and patterns for the volcanic ash. Two sulfate plumes originating from Mt. Etna arrived over Antikythera: The first plume is centered approximately at 1.5-2.5 km height during 17:00-18:00 UTC on 2 June 2019, whereas the second sulfates plume is centered between 2 and 3.5 km on 3 June 2019, 00:00-02:00 UTC, mixed with dust particles from North Africa and the already existing dust particles that were located above the Mediterranean Sea since 31 May and 1 June 2019. On 3 June 2019 after 18:00 UTC, volcanicash particles emitted from Mt. Etna volcano were also transported above the PANGEA station at altitudes 0-1 km, but at small concentrations. Cluster analysis with the HYSPLIT model was performed in the study to showcase the strong connection between Mt Etna and the PANGEA observatory. Conclusions from this study are summarized as follows: 1: The HYSPLIT cluster analysis indicates that the upper tropospheric-lower stratospheric air masses from Etna are mainly transported eastwards over the Mediterranean and are detected at the new PANGEA observatory of NOA at the island of Antikythera establishing the Etna-Antikythera connection. The PANGEA observatory is located 765 km downwind the volcano and presents an important infrastructure for the monitoring of volcanic emissions from Etna.

2:
The long-range dispersion of sulfates from Etna was simulated with FLEXPART-WRF using an emission rate of 4 kt/day. The sulfates plume spreads mainly northeastward from the volcano and takes a circular shape due to a passing cyclone, while on 3 June 2019, the plume has covered the northeastern parts of Greece, with the main part shifting towards the southern parts of the country and reaching the Antikythera station. This is consistent with the observed movement of the simulated sulfates plume as depicted from TROPOMI, as well as with the MERRA-2 reanalysis products. 3: FLEXPART-WRF simulations for volcanic ash were performed with mass values of the order of 10 5 kg/s following the studies of weak plumes with wind conditions [23,64]. The transport and the different followed paths of sulfates and volcanic ash driven by FLEXPART-WRF show that both plumes move eastward but only sulfate plumes reach the southern parts of Greece. The modeled sulfate mass concentration is approximately 18 μg/m 3 for the first plume and 20 μg/m 3 for the second plume on 2 and 3 June 2019, respectively, whereas the ash mass concentration is below 5 μg/m 3 on 3 June 2019, at 18:00 UTC. 4: The height and mass concentration of the simulated two sulfate plumes were evaluated in a qualitative manner using PANGEA measurements. According to the Polly XT -NOA observations, the two elevated plumes are located between 1 and 2 km and 2 and 3 km above the local PBL on 2 June 2019, 17:00-18:00 UTC and 3 June 2019, 00:00-02:00 UTC respectively, three and four days after the eruption. For this time windows, the possible contribution of ash and sulfate particles to the lidar backscatter coefficient profile β at 532 nm was separated based on the POLIPHON technique. For the sulfates mass concentration, the agreement between the model and the lidar is satisfactory, with the depth of the two sulfate layers to be well captured by the model. On the contrary, the volcanic ash plume is not accurately reproduced by the model for the selected time interval. By comparing the modeled ash and dust mass concentrations to lidar profiles the agreement with the lidar is more satisfactory for dust particles. 5: Finally, the combined information of the backward trajectory analysis, the sourcereceptor relationships and the results of the WRF-Chem model at Antikythera station, on 3 June 2019, at 00:00 UTC indicate the presence of a mixture of volcanic sulfates and dust particles.
The above complex interactions between different aerosol and gas species and the meteorological patterns of the Mediterranean indicate the need to intensify the relevant research activities in this area. As an important step in this direction, the development of PANGEA observatory in Antikythera offers the chance for continuous monitoring of atmospheric tracers. As shown in this study, volcanic ash and sulfate particles may coexist with dust, marine, biogenic, and anthropogenic aerosols over the Mediterranean. This gives us the opportunity to extend our future research by applying also quantitative inversion algorithms (e.g., [96]) to largely improve the predictions of volcanic ash fluxes from Etna eruptions and also the separation between different species constrained by ground-based and satellite measurement data. Funding: This research was funded by EU H2020 E-shape project (Grant Agreement n. 820852) and by the project "PANhellenic infrastructure for Atmospheric Composition and climatE change" (MIS 5021516) which is implemented under the Action "Reinforcement of the Research and Innovation Infrastructure", funded by the Operational Programme "Competitiveness, Entrepreneurship and Innovation" (NSRF 2014-2020) and co-financed by Greece and the European Union (European Regional Development Fund)". This research was also funded by "Stavros Niarchos Foundation (SNF)".

Special Issue Statement:
This article is part of the special issue: "Evaluation and Optimization of Atmospheric Numerical Models".