From Tropospheric Folding to Khamsin and Foehn Winds: How Atmospheric Dynamics Advanced a Record-Breaking Dust Episode in Crete

A record-breaking dust episode took place in Crete on 22 March 2018. The event was characterized by surface concentrations exceeding 1 mg m−3 for a period of 4–7 h, reaching record values higher than 6 mg m−3 at the background station of Finokalia. We present here a detailed analysis of the atmospheric dynamical processes during this period, to identify the main reasons for such extreme dust advection over Crete. At the synoptic scale, the weakening of the polar vortex and the meridional transport of polar air masses at upper tropospheric layers resulted in a strong jet streak over north Africa and Central Mediterranean and corresponding tropospheric folding that brought cold stratospheric air in mid and upper troposphere. Cyclogenesis occurred at the Gulf of Sirte in Libya, resulting in strong winds over the north-east parts of Libya, enhancing particle emissions. The dust plume traveled at low altitude (0.5–3 km) along the warm conveyor belt preceding the depression cold front. This type of dusty southerly wind is commonly known as “Khamsin”. As the flow approached Crete, Foehn winds at the lee side of the island favored the downward mixing of dust towards the surface, resulting in local maxima of PM10 in Heraklion and Finokalia. The analysis is based on the combination of high-resolution WRF-Chem simulations reaching up to 1 × 1 km grid space over Crete, ground-based and satellite remote sensing of the dust plumes (PollyXT LiDAR, MSG-SEVIRI, MODIS) and detailed surface aerosol in situ measurements at urban (Heraklion, Chania, Greece) and background (Finokalia) stations in Crete.


Introduction
Mineral dust particles that are emitted from the desert areas in the Sahara and the Middle East travel long distances in the atmosphere and affect air quality, weather, climate and local ecosystems at the Mediterranean basin [1][2][3][4]. Dust affects the radiative transfer in the area [5][6][7], and cloud processes due to the activation of dust particles as cloud condensation nuclei (CCN) and ice nuclei (IN) [8,9]. Moreover, the ocean depositions of dust affect biogeochemical processes and phytoplankton blooming [10,11]. In terms of health impacts, increased concentrations of fine dust aerosols may cause respiratory diseases and other related health issues [12][13][14][15]. The island of Crete is very often affected by windblown dust originating from the Sahara, due to its proximity to the Africa coastline. Regarding health effects, the severity of the dust events is mainly dictated by the amounts of dust aerosol near the surface. Such conditions with increased surface PM 10 of Saharan origin are more often found in Crete when the air masses originate from the nearby dust sources of Libya [16].
In general, the most crucial parameter for the emission of dust is the near-surface wind. Increased wind speeds may occur due to synoptic wind forcing, topographic effects (e.g., valley channeling), low-level jets (LLJ) [17] or squall lines and storm downdrafts [18,19]. However, most of the above processes result in detached elevated dust plumes over the Mediterranean. Most dust layers in the area are observed at heights of 4-5 km in the troposphere and are associated either with Mediterranean low-pressure systems or with the summer anticyclonic circulation over north Africa [20][21][22][23]. For the surface PM 10 concentration in Crete to reach significant levels, the mesoscale circulation needs to favor transport directly from the coastal dust sources of Libya towards Crete. This type of flow is commonly known as "Khamsin" winds in Libya and Egypt and is mostly frequent during spring months [24]. Such conditions result in the generation of dense and low altitude dust clouds that travel inside the lower troposphere and can reach the south parts of Crete in less than 12 h after emission. In such cases Crete's mountain range acts as a perpendicular barrier to the southerly winds and Foehn downslope winds may develop on the lee side of the island causing a significant downward flux. These conditions have significant implications for human health as they produce high near-surface concentrations of dust, raise surface temperatures and reduce relative humidity [25][26][27].
Such meso-γ weather phenomena are not accurately represented in global and limited area atmospheric models mainly due to the poor description of topography at meso-β and meso-α grid increment where most topographic elements are sub-grid features. Even at a grid-space of e.g., 5-10 km which is common for weather forecasting purposes, the actual resolved topography is at a scale of 4∆x [28], implying that the explicit resolving of Foehn winds would require a grid increment of 1 × 1 km or finer. In this work we investigate the weather processes from synoptic down to meso and local atmospheric scales during a record-breaking dust event that took place in March 2018 in Crete based on synergistic analysis of high-resolution modeling, remote sensing and in situ data. The study is divided in 3 sections. Section 2 includes the description of methodological and analysis tools. In Section 3 we present the main findings regarding the multi-scale atmospheric dynamics and the dust measurements during the event. Section 4 is a summary and discussion of the main findings.

Methodology
The methodology consists of combined observational and modeling approach, using the resources of the Greek National Research Infrastructure (RI) PANACEA (PANACEA-PANhellenic infrastructure for Atmospheric Composition and climate change), operating the national facilities of the ACTRIS RI (Aerosols, Clouds, and Trace gases Research Infrastructure: https://www.actris.eu/). where ρ is the particle mass density, r e f f is the effective radius and n denotes the corresponding WRF sub-bin. The extinction efficiencies at 550 nm (Table 1) are calculated with Mie theory [43,44], using a refractive index of 1.55 + i0.005 for dust [45]. We consider lognormal size distributions for the sub-bins with minimum-maximum radius, effective radius and geometric standard deviations as shown in Table 1. Gravitational settling and surface deposition of dust is based on the parameterization of Wesely, 1989 [46]. Wet deposition of dust is not parameterized in these simulations; however, it is not considered important since there was no significant rainfall along the major transport path. We performed a 3 day spin up run that is considered sufficient for establishing a regional dust background and for providing the initial dust conditions for the subsequent simulations, taking also into account that the main dust event has a relatively limited geographical and temporal extent. Table 1. Extinction efficiencies at 550 nm (Q ext550 ) for the dust sub-bins of WRF-CHEM. The sub-bin characteristics are also provided, in terms of the minimum (r min ) and maximum (r max ) radius, the effective radius (r e f f ) and the geometric standard deviation (σ g ) of the sub-bin lognormal size distributions. . The LiDAR is equipped with three elastic channels at 355 nm, 532 nm and 1064 nm, two rotationalvibrational Raman channels at 387 nm and 607 nm, two linear depolarization channels at 355 nm and 532 nm, and one water vapor channel at 407 nm. The combined use of its near field and far field telescopes provides reliable vertical profiles of aerosol optical properties from 0.25 km to 10 km in height.

MSG-SEVIRI-Dust RGB
The Meteosat Second Generation (MSG) Spinning Enhanced Infrared and Visible Imager (SEVIRI) false color dust red-green-blue (RGB) imagery uses channels 7, 9 and 10 (at 8.7 µm, 10.8 µm and 12.0 µm respectively). The 12 µm and 10.8 µm channels are able to discriminate between land surface temperatures and cloud temperatures. The addition of 8.7 µm channel provides information about trace gases and aerosols. The values assigned to the Red, Green and Blue components for generating Dust RGB composite images are given in Table 2. Each of the components is scaled using the following equation: where C i is the value for red, green and blue components; [4,254] is the color range; δX = max(X) − min(X) is the number of colors; and γ i is the respective gamma factor. Dust RGBs have been successfully used to monitor dust and volcanic ash plumes over land and sea both during day and night. The infrared channels used in the RGB composites are sensitive to various surface and atmospheric properties-the Red and the Blue components in this case determine the presence of dust (i.e., warm with large negative BT difference between 10.8 and 12 micron channels) and therefore depicted as bright magenta (during day) or purple (during night) color over land. A dusty atmosphere can also be tracked over water as a magenta color. The SEVIRI Dust RGB imageries are commonly used as a proxy measure of dust presence in various applications; however, it is worth mentioning that interpretation of the dust signals (dark magenta) can be challenging for different meteorological situations. For example, the RGB signal can be severely suppressed by the presence of high atmospheric column water vapor [48]. Further, the RGB signal can actually be reversed, such that it appears less pink with increasing dust, in the presence of a surface/lower tropospheric temperature inversion. Additionally, dust mineralogy and dust layer height can also modulate the RGB signal to some extent. The strengths and the shortcomings of SEVIRI dust RGB images are highlighted by [48] in detail.  [49,50] (https://ladsweb.modaps.eosdis.nasa.gov/). In the present analysis we use the Level 2 (L2) data of the aerosol optical depth at 550 nm (AOD 550nm ), acquired from the latest version (Collection 6.1, C061) of the MODIS retrieval algorithm, corresponding to satellite overpasses (swaths) of five minutes time interval and their nominal spatial resolution is 10 km × 10 km at nadir. For our purposes, we have used the "combined" AOD produced via the implementation of the Dark Target (DT) algorithm applied based on different assumptions over maritime areas [51,52] and vegetated land [52,53] while above bright surfaces (i.e., deserts) and transition zones situated between arid and vegetated land areas the Deep Blue (DB; enhanced Deep Blue) approach is followed [54,55]. The merging procedure of the DT-Land, DT-Ocean and DB AOD products is described in [56]. Global evaluation studies have shown that the expected error of MODIS DT AOD 550nm above sea lies between +(0.04 + 0.1τ A ) and −(0.02 + 0.1τ A ) [52] while above land the corresponding range is ±(0.05 + 0.15τ A ) [52,57], where τ A stands for AERONET AOD. For the DB AOD, an absolute uncertainty of approximately (0.086 + 0.56τ M )/AMF, where AMF is the geometric air mass factor and τ M is the MODIS AOD 550nm , has been reported by [58]. In contrast to SEVIRI, the frequency sampling of the MODIS spectroradiometer is lower since it flies onboard polar orbit satellites (i.e., Terra and Aqua). More specifically, AOD retrievals are available only during early morning hours (Terra satellite overpass) and late noon (Aqua satellite overpass). In addition, AOD observations are not possible over areas where clouds exist as well as in oceanic regions where the boundary reflectance is maximized (i.e., sun-glint).

In-Situ Measurements
Measurements were performed on the island of Crete at the environmental research station of the University of Crete at Finokalia (35 • 20 N, 25 • 40 E). Finokalia environmental research station was established in 1993 and it has been characterized as representative of the background regional conditions for the eastern Mediterranean basin [59]. The long-range transport of mineral dust from north Africa is the main reason for exceedance of the EU limit value of 50 µg m −3 (Directive 2008/50/EU) at the Finokalia station, when the air mass back trajectories originate from northern Africa [3]. Such events have taken place for around 5% of the days since 2004. Additional measurements were available at the atmospheric quality measurement station of the Region of Crete at the cities of Heraklion (35 • 19 N, 25 • 7 E) and Chania (32 • 5 N, 20 • 16 E). The PM 10 measurements were performed with two identical Thermo ESM Andersen FH 62 I-R beta attenuation particulate monitors with a time resolution of 5 min [14]. Both instruments have been compared to gravimetric techniques in the past with good results [60].

Atmospheric Dynamics
To understand the main mechanisms leading to this unusual event, we examine the major atmospheric processes taking place in the area starting from synoptic scale dynamics and moving to mesoscale and local dynamics as resolved by WRF-Chem model. At first, the weakening of the Arctic polar vortex allows the meridional transport of cold polar air towards Europe on 21 March 2018. This process is evident in the FNL analysis of 300 hPa temperature in Figure 2a. A streak of very low temperatures (lower than −50 °C) is advected across Europe and the upper level polar air mass reaches 30° N in the West Mediterranean by 1200 UTC, 21 March 2018. WRF-Chem fields from the outer 12 × 12 km grid are used in the following sections to analyze the impacts of this intrusion in atmospheric circulation. Associated with the thermal gradient there is a southwesterly 300 hPa jet streak present at 0600 UTC, 22 March 2018 with maximum wind speeds over 70 ms −1 (Figure 2b). The left exit region of the jet streak is found over the Gulf of Sirte and Central Mediterranean, denoting upper level divergence at this area. Tropospheric folding and the intrusion of stratospheric air at Mediterranean and North Africa latitudes is clearly shown by the potential vorticity (PV) and potential temperature (θ) plots in Figure 2c,d. Such PV anomalies have long been recognized as diagnostic quantities for identifying tropospheric folding in the atmosphere e.g., [62,63]. In Figure 2c the stratospheric air mass is spread southwestwards along the isentropic (iso-θ) surface of 315 K

Atmospheric Dynamics
To understand the main mechanisms leading to this unusual event, we examine the major atmospheric processes taking place in the area starting from synoptic scale dynamics and moving to mesoscale and local dynamics as resolved by WRF-Chem model. At first, the weakening of the Arctic polar vortex allows the meridional transport of cold polar air towards Europe on 21 March 2018. This process is evident in the FNL analysis of 300 hPa temperature in Figure 2a. A streak of very low temperatures (lower than −50 • C) is advected across Europe and the upper level polar air mass reaches 30 • N in the West Mediterranean by 1200 UTC, 21 March 2018. WRF-Chem fields from the outer 12 × 12 km grid are used in the following sections to analyze the impacts of this intrusion in atmospheric circulation. Associated with the thermal gradient there is a southwesterly 300 hPa jet streak present at 0600 UTC, 22 March 2018 with maximum wind speeds over 70 ms −1 (Figure 2b). The left exit region of the jet streak is found over the Gulf of Sirte and Central Mediterranean, denoting upper level divergence at this area. Tropospheric folding and the intrusion of stratospheric air at Mediterranean and North Africa latitudes is clearly shown by the potential vorticity (PV) and potential temperature (θ) plots in Figure 2c,d. Such PV anomalies have long been recognized as diagnostic quantities for identifying tropospheric folding in the atmosphere e.g., [62,63]. In Figure 2c the stratospheric air mass is spread southwestwards along the isentropic (iso-θ) surface of 315 K characterized by PV up to 10 PVU (1 PVU= 10 −6 m 2 s −1 K kg −1 ). Stratospheric PV (higher than 2 PVU) is evident at the poleward side of the polar jet and tropospheric PV (less than 2 PVU) is found at the equatorward side. The downward transport of momentum towards the surface is facilitated by the further equatorward penetration of the PV anomaly over north Africa interacting with the developing dry convective boundary layer during the morning hours. As shown by the West-East cross section in Figure 2d, at 0600 UTC the core of the polar jet is at around 11-12 km above 11 • E (black contours higher than 70 ms −1 ). The upward bowing of the 305 K isentropic line west of 11 • E shows the cold core structure and the low static stability at this area that allows mixing of the winds from the top of the boundary layer towards the surface. This downward mixing of high momentum and cold air mass enforces the eastward advection of the surface cold front. The structure of the eastward propagating trough is shown in Figure 3a  Regarding dust emissions, the key factor is the enforcement of the boundary layer turbulent kinetic energy (TKE) [68]. This is evident at the 925 hPa map at 0900 UTC ( Figure 3b) where the highest TKE values (up to 2 m 2 s −2 ) are found at northeast Libya coinciding with the generation area of the dust storm in the model as shown in the following paragraphs.
The location of cyclogenesis is found on 21 March 2018, 1800 UTC over the east part of the Gulf of Sirte, as seen by the sea level pressure (slp) minimum of 1000 hPa at this area ( Figure 4a). Tropospheric folding is also identified in Figure 4 by the eastward propagation of the stratospheric intrusion and the minimum tropopause heights (geopotential at 2 PVU iso-surface) over north Africa and Central Mediterranean. The depression intensifies as it moves northwards towards the Ionian Sea and reaches minimum slp of 990 hPa at the Gulf of Taranto at 1200 UTC, 22 March 2018 ( Figure 4b). The location of the cold front is identified by the abrupt change in wind direction and the temperature gradient along the 15 • E meridian in Figure 5a at 0600 UTC, 22 March 2018. The warm conveyor belt ahead of the cold front is actually the main driving pathway for the mobilization and transport of dust from Libya towards Crete. As seen in Figure 5b, the near-surface dust concentration exceeds 15 mg m −3 at certain frontal areas and the low-level convergence inside the warm sector actually traps the dust particles ahead of the cold front.
Regarding dust emissions, the key factor is the enforcement of the boundary layer turbulent kinetic energy (TKE) [68]. This is evident at the 925 hPa map at 0900 UTC (Figure 3b) where the highest TKE values (up to 2 m 2 s −2 ) are found at northeast Libya coinciding with the generation area of the dust storm in the model as shown in the following paragraphs.

Dust Transport
The low-pressure system is slow moving on 22 March 2018 and the location of the conveyor belt preceding the cold front forms a high-speed pathway connecting the Sahara with Crete. This transport path is shown in Figure 6a

Dust Transport
The low-pressure system is slow moving on 22 March 2018 and the location of the conveyor belt preceding the cold front forms a high-speed pathway connecting the Sahara with Crete. This transport path is shown in Figure 6a-d at 0900, 1200, 1500 and 1800 UTC respectively, on 22 March 2018. The dust plume is constrained inside the warm sector of the depression with dust loads up to 20 g m −2 inside the core of the system (i.e., Figure 6b). The dusty Khamsin flow affects Crete between 1200-1800 UTC, 22 March 2018.

Dust Transport
The low-pressure system is slow moving on 22 March 2018 and the location of the conveyor belt preceding the cold front forms a high-speed pathway connecting the Sahara with Crete. This transport path is shown in Figure 6a

Dust Concentrations
The severity of the event is shown by satellite and model results in Figure 7. MODIS-Terra overpass at 0950 UTC captures an area of AOD up to 5, extending from the coast of Libya towards the cloud-covered region north of 34° N (Figure 7a). Furthermore, AOD of 5 is the maximum value in MODIS classification [69] implying that greater values could be actually present at the area. The corresponding model AOD map in Figure 7b exhibits a good spatial agreement with MODIS in dust dispersion but the maximum AOD over the sea does not exceed 3.5 implying an underestimation of the modeled optical depth with respect to MODIS. This can be attributed to an underestimation of dust emissions or to latency in dust transport since inland model values exceed 5 and during the following hours the modeled AOD over the sea exceeds 6 (not shown).

Dust Concentrations
The severity of the event is shown by satellite and model results in Figure 7. MODIS-Terra overpass at 0950 UTC captures an area of AOD up to 5, extending from the coast of Libya towards the cloud-covered region north of 34 • N (Figure 7a). Furthermore, AOD of 5 is the maximum value in MODIS classification [69] implying that greater values could be actually present at the area. The corresponding model AOD map in Figure 7b exhibits a good spatial agreement with MODIS in dust dispersion but the maximum AOD over the sea does not exceed 3.5 implying an underestimation of the modeled optical depth with respect to MODIS. This can be attributed to an underestimation of dust emissions or to latency in dust transport since inland model values exceed 5 and during the following hours the modeled AOD over the sea exceeds 6 (not shown).
At the time when the dust plume approaches Crete, the cold front has already moved eastwards so that only the central and east parts of the island are affected by the advected dust particles (see also Figure 6). Surface concentrations in the model exceed 6000 µg m −3 especially towards the east parts of Crete and near Finokalia station ( Figure 8). As seen in Figure 8d the maximum concentrations are found at the leeward side of the highest mountain ridges in Crete, suggesting a downward mixing mechanism i.e., Foehn downslope winds at these areas. This hypothesis is confirmed by the vertical model cross sections over Finokalia in Figure 9. The presence of a Foehn mechanism is identified in Figure 9a by the descending of the isentropes at the lee side of the mountain range. The low-level flow blocking produces a barrier jet and the subsequent isentropic draw down of free tropospheric air on the lee side results in downward mixing of the dusty layers increasing the near-surface concentrations up to 6000 µg m −3 . The wind speeds reach up to 30 m s −1 near the surface and the adiabatic warming of the descending air mass leads to air temperatures exceeding 30 • C at the lee side and to horizontal temperature gradients of more than 10 • C between the windward and leeward ridge sides (Figure 9b). in MODIS classification [69] implying that greater values could be actually present at the area. The corresponding model AOD map in Figure 7b exhibits a good spatial agreement with MODIS in dust dispersion but the maximum AOD over the sea does not exceed 3.5 implying an underestimation of the modeled optical depth with respect to MODIS. This can be attributed to an underestimation of dust emissions or to latency in dust transport since inland model values exceed 5 and during the following hours the modeled AOD over the sea exceeds 6 (not shown). At the time when the dust plume approaches Crete, the cold front has already moved eastwards so that only the central and east parts of the island are affected by the advected dust particles (see also Figure 6). Surface concentrations in the model exceed 6000 μg m −3 especially towards the east parts of Crete and near Finokalia station ( Figure 8). As seen in Figure 8d the maximum concentrations are found at the leeward side of the highest mountain ridges in Crete, suggesting a downward mixing mechanism i.e., Foehn downslope winds at these areas. This hypothesis is confirmed by the vertical model cross sections over Finokalia in Figure 9. The presence of a Foehn mechanism is identified in Figure 9a by the descending of the isentropes at the lee side of the mountain range. The low-level flow blocking produces a barrier jet and the subsequent isentropic draw down of free tropospheric air on the lee side results in downward mixing of the dusty layers increasing the near-surface concentrations up to 6000 μg m −3 . The wind speeds reach up to 30 m s −1 near the surface and the adiabatic warming of the descending air mass leads to air temperatures exceeding 30 °C at the lee side and to horizontal temperature gradients of more than 10 °C between the windward and leeward ridge sides (Figure 9b).

Effects of Model Grid Space on Resolved Dynamics
This local-scale phenomenon partially explains the poor representation of this event in operational dust models. A sensitivity experiment, shown in Figure 8, clearly displays the effects of resolved topography in modeled surface dust concentrations. Four different model runs are compared in Figure 8 with identical configuration either than model grid-space. At 24 × 24 km space (Figure 8a), the model hardly captures Crete topography with a maximum height of 400 m (when the actual highest top is at 2400 m). In this run the dust pattern is literally uninterrupted by the island and the maximum concentration hardly exceeds 2000 μg m −3 . At 12 × 12 km grid space ( Figure  8b), the maximum top is at 1000 m, a slight transformation of the flow is evident due to the orographic barrier and an area of maximum dust concentration exceeding 3000 μg m −3 appears at the lee side towards the eastern part of Crete. The 3 × 3 km and 1 × 1 km runs (Figures 8c,d) both resolve the Foehn mechanism and clearly present maximum concentrations at the lee sides of the mountain ridges (notice also the different topographic contours in these plots). However, the dust concentrations in the finest run are even higher than the 3 × 3 km run. These results justify the need for high-resolution simulations in order to resolve similar events over complex topography, but also raise an interesting consideration on how much more could the model improve by further increasing the grid resolution.

Comparison with Observations
To examine the spatiotemporal evolution of the extreme dust episode at the station of Finokalia, we analyze the modeled vertical time-height properties of dust concentration, temperature and relative humidity from the finest 1 × 1 km run over the station (Figure 10a). The peak of the event is between 1615-1715 UTC when the modeled surface concentration exceeds 6000 μg m −3 . At this time the core of the plume is located at 0.9-1.3 km a.s.l. where dust concentration is higher than 8000 μg m −3 . The downward mixing of the dry and warm Saharan air mass is evident in Figure 10 by the abrupt increase in ambient temperature reaching up to 32 °C at 1700 UTC and by the 10% relative humidity line that actually constraints the downslope air mass. Comparison of the modeled against measured meteorological values at Finokalia shows that the simulated surface values of temperature and relative humidity are in good agreement and in general following the time-evolution of the in situ measurements (Figure 10b). An interesting finding arises from the comparison between modeled and measured wind properties. The model reproduces correctly the passage of the cold front. In agreement with the observations the wind is veering from SSE (150°) to WNW (300°) and is also reduced to calm conditions (1.2 m s −2 ) behind the front (after 1700 UTC). However, the downslope winds are underestimated by almost 10 m s −1 at the station location (Figure 10b before

Effects of Model Grid Space on Resolved Dynamics
This local-scale phenomenon partially explains the poor representation of this event in operational dust models. A sensitivity experiment, shown in Figure 8, clearly displays the effects of resolved topography in modeled surface dust concentrations. Four different model runs are compared in Figure 8 with identical configuration either than model grid-space. At 24 × 24 km space (Figure 8a), the model hardly captures Crete topography with a maximum height of 400 m (when the actual highest top is at 2400 m). In this run the dust pattern is literally uninterrupted by the island and the maximum concentration hardly exceeds 2000 µg m −3 . At 12 × 12 km grid space (Figure 8b), the maximum top is at 1000 m, a slight transformation of the flow is evident due to the orographic barrier and an area of maximum dust concentration exceeding 3000 µg m −3 appears at the lee side towards the eastern part of Crete. The 3 × 3 km and 1 × 1 km runs (Figure 8c,d) both resolve the Foehn mechanism and clearly present maximum concentrations at the lee sides of the mountain ridges (notice also the different topographic contours in these plots). However, the dust concentrations in the finest run are even higher than the 3 × 3 km run. These results justify the need for high-resolution simulations in order to resolve similar events over complex topography, but also raise an interesting consideration on how much more could the model improve by further increasing the grid resolution.

Comparison with Observations
To examine the spatiotemporal evolution of the extreme dust episode at the station of Finokalia, we analyze the modeled vertical time-height properties of dust concentration, temperature and relative humidity from the finest 1 × 1 km run over the station (Figure 10a). The peak of the event is between 1615-1715 UTC when the modeled surface concentration exceeds 6000 µg m −3 . At this time the core of the plume is located at 0.9-1.3 km a.s.l. where dust concentration is higher than 8000 µg m −3 . The downward mixing of the dry and warm Saharan air mass is evident in Figure 10 by the abrupt increase in ambient temperature reaching up to 32 • C at 1700 UTC and by the 10% relative humidity line that actually constraints the downslope air mass. Comparison of the modeled against measured meteorological values at Finokalia shows that the simulated surface values of temperature and relative humidity are in good agreement and in general following the time-evolution of the in situ measurements (Figure 10b). An interesting finding arises from the comparison between modeled and measured wind properties. The model reproduces correctly the passage of the cold front. In agreement with the observations the wind is veering from SSE (150 • ) to WNW (300 • ) and is also reduced to calm conditions (1.2 m s −2 ) behind the front (after 1700 UTC). However, the downslope winds are underestimated by almost 10 m s −1 at the station location (Figure 10b before 1500 UTC). This comparison emphasizes on the very local nature of the event and on the large gradients in wind speed induced by the complex topography. For example, values greater than 30 m s −1 are only simulated at about 0.5 km above Finokalia station, implying that in reality the station was inside the Saharan layer and even the 1 × 1 km grid space is not sufficient for accurately reproducing all atmospheric processes near the surface. Additional limitations may be related to the improper reproduction of wind gusts or to assumptions made by the boundary layer parameterization of the model. Unfortunately, the signal of the PollyXT LiDAR is most probably attenuated during the specific event. A time window with reliable LiDAR profiles is only found at 0300 and 0400 UTC on 22 March (Figure 11a). This timeframe corresponds to a different air mass, also carrying Saharan dust that arrives at layers 1.5-6 km over the station. The Klett-Fernald-Sasano method was utilized for the derivation of the extinction coefficient profiles at 532 nm, using a LiDAR ratio of S = 55. The resulting AODs are 2.93 and 1.07 at 0300 and 0400 UTC respectively. The dust particle mass concentration profiles of the model (solid line) and PollyXT LiDAR (dashed line) for 0300 and 0400 UTC (green and blue line respectively) are also calculated from these profiles. The dust plume was assumed to be Unfortunately, the signal of the PollyXT LiDAR is most probably attenuated during the specific event. A time window with reliable LiDAR profiles is only found at 0300 and 0400 UTC on 22 March (Figure 11a). This timeframe corresponds to a different air mass, also carrying Saharan dust that arrives at layers 1.5-6 km over the station. The Klett-Fernald-Sasano method was utilized for the derivation of the extinction coefficient profiles at 532 nm, using a LiDAR ratio of S = 55. The resulting AODs are 2.93 and 1.07 at 0300 and 0400 UTC respectively. The dust particle mass concentration profiles of the model (solid line) and PollyXT LiDAR (dashed line) for 0300 and 0400 UTC (green and blue line respectively) are also calculated from these profiles. The dust plume was assumed to be coarse-mode dominated, thus the mass-specific extinction coefficient value of k ext = 0.52 m 2 g −1 was utilized for the LiDAR mass concentration profiles [70]. Comparison of modeled dust profiles with these LiDAR retrievals ( Figure 11b) indicates a satisfactory agreement between modeled and observed concentrations in terms of the event scale. The model seems to follow also the modification in maximum concentrations between this 2-h period. However the depth of the layer is not accurately reproduced especially at 0300 UTC when the LiDAR backscatter signal shows a uniform layer of 1000-1600 µg m −3 while the modeled profile quickly declines to values less than 1000 µg m −3 above 2 km. Overall, this comparison may not be representative of the main event, but it generally increases our confidence on the model simulation suggesting also an under-prediction of the actual dust concentrations in the model. coarse-mode dominated, thus the mass-specific extinction coefficient value of = 0.52 ⋅ g was utilized for the LiDAR mass concentration profiles [70]. Comparison of modeled dust profiles with these LiDAR retrievals (Figure 11b) indicates a satisfactory agreement between modeled and observed concentrations in terms of the event scale. The model seems to follow also the modification in maximum concentrations between this 2-h period. However the depth of the layer is not accurately reproduced especially at 0300 UTC when the LiDAR backscatter signal shows a uniform layer of 1000-1600 μg m −3 while the modeled profile quickly declines to values less than 1000 μg m −3 above 2 km. Overall, this comparison may not be representative of the main event, but it generally increases our confidence on the model simulation suggesting also an under-prediction of the actual dust concentrations in the model. The surface concentrations of PM10 in the model are compared with the in situ measurements at the stations of Finokalia, Heraklion and Chania ( Figure 12). Dust is the only aerosol considered in these simulations and the modeled PM10 is calculated from the five dust bins as: PM10 = bin1 + bin2 + bin3 + 0.87•bin4. The local nature of the event is clearly illustrated by the in situ measurements (dashed lines) in this plot. Finokalia and Heraklion stations, which are located towards the east part of Crete, recorded remarkably increased PM10 concentrations exceeding 6 mg m −3 during 1400-1800 UTC. At the same time, PM10 at Chania station that is located 140 km to the west was below 0.2 mg m −3 . Taking into account the complexity of the event the model reproduces these record-breaking values in a satisfactory manner, but with several biases. The time evolution of the dust-storm is accurately described and the maximum modeled PM10 are also found during 1400-1800 UTC in Heraklion and Finokalia. The coevolution of measured and modeled values at these two stations is also evident in Figure 12 and the correlation coefficient is 0.43 for Finokalia and 0.56 for Heraklion. The narrow maximum spike of modeled PM10 in Finokalia reaches 5.5 mg m −3 one hour earlier than the actual recorded peak of 6.34 mg m −3 . However the model clearly underestimates PM10 in Heraklion since although the simulation produces very high values exceeding 2.5 mg m −3 at this station, the actual concentrations were actually two times higher exceeding 4 mg m −3 . In Chania station, the model fails to predict the very low recorded PM10 and thus implies a more spatially extended episode than what actually occurred. This is probably explained due to the underprediction of Lefka Ori topography in the model. This mountain range extends south of The surface concentrations of PM 10 in the model are compared with the in situ measurements at the stations of Finokalia, Heraklion and Chania ( Figure 12). Dust is the only aerosol considered in these simulations and the modeled PM 10 is calculated from the five dust bins as: PM 10 = bin1 + bin2 + bin3 + 0.87·bin4. The local nature of the event is clearly illustrated by the in situ measurements (dashed lines) in this plot. Finokalia and Heraklion stations, which are located towards the east part of Crete, recorded remarkably increased PM 10 concentrations exceeding 6 mg m −3 during 1400-1800 UTC. At the same time, PM 10 at Chania station that is located 140 km to the west was below 0.2 mg m −3 . Taking into account the complexity of the event the model reproduces these record-breaking values in a satisfactory manner, but with several biases. The time evolution of the dust-storm is accurately described and the maximum modeled PM 10 are also found during 1400-1800 UTC in Heraklion and Finokalia. The coevolution of measured and modeled values at these two stations is also evident in Figure 12 and the R 2 correlation coefficient is 0.43 for Finokalia and 0.56 for Heraklion. The narrow maximum spike of modeled PM 10 in Finokalia reaches 5.5 mg m −3 one hour earlier than the actual recorded peak of 6.34 mg m −3 . However the model clearly underestimates PM 10 in Heraklion since although the simulation produces very high values exceeding 2.5 mg m −3 at this station, the actual concentrations were actually two times higher exceeding 4 mg m −3 . In Chania station, the model fails to predict the very low recorded PM 10 and thus implies a more spatially extended episode than what actually occurred. This is probably explained due to the underprediction of Lefka Ori topography in the model. This mountain range extends south of Chania and acts as a perpendicular barrier to the flow from Africa. It consists of 58 tops exceeding 2 km and even at this high 1 × 1 km grid space this pattern is not accurately reproduced. As a result, the core of the dust plumes arriving in the model at 1-2 km, instead of being deflected by the mountain range, they actually cross the ridge resulting in downslope winds at the lee side that affect the west parts of Crete. Actually the coarse 24 × 24 km model run (Figure 8a) predicts very well the PM 10 values in Chania but for the wrong reasons. In this run the model hardly sees the Crete Mountains and as a result the dust clouds continue travelling detached towards the north without increasing the surface dust concentration.
Atmosphere 2018, 9, x FOR PEER REVIEW 15 of 20 Chania and acts as a perpendicular barrier to the flow from Africa. It consists of 58 tops exceeding 2 km and even at this high 1 × 1 km grid space this pattern is not accurately reproduced. As a result, the core of the dust plumes arriving in the model at 1-2 km, instead of being deflected by the mountain range, they actually cross the ridge resulting in downslope winds at the lee side that affect the west parts of Crete. Actually the coarse 24 × 24 km model run (Figure 8a) predicts very well the PM10 values in Chania but for the wrong reasons. In this run the model hardly sees the Crete Mountains and as a result the dust clouds continue travelling detached towards the north without increasing the surface dust concentration.

Summary and Conclusions
We analyzed the most intense Saharan dust episode ever recorded at Finokalia station in Crete during its 14 years operational period (2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018). Analysis of the atmospheric processes during this episode suggests that the unique PM10 concentrations recorded in Crete on 22 March 2018 were produced by a combination of meteorological processes taking place at various atmospheric scales. In general, the main driver for strong dust events in Mediterranean is atmospheric dynamics and this particular event is a typical example that combines several of the dynamical processes commonly involved in dust episodes: At synoptic scale, the weakening of the polar vortex allows descending of cold stratospheric air of high potential vorticity towards lower latitudes at Mediterranean and north Africa.
At mesoscale, the tropospheric folding favors cyclogenesis and intensifies the low-pressure systems along the surface baroclinic zones. The formation of mid-latitude cyclones at this area results in mobilization of dust particles from the Sahara and the scale and movement of these Mediterranean depressions determines also the dust pathways towards Europe.
At local scale, the vertical stability of the dusty air mass is disturbed when the plumes reach the orographic barriers of Mediterranean islands (e.g., Crete, Cyprus, Sicily) or the inland continental mountain ranges. Depending on the relevant heights between the mountain ridge and the elevated plumes, topography may act either as a stopping barrier that leads to increased dust concentrations at the windward ridge [22] or as a primary reason for the generation of downslope Foehn winds and downward mixing of dust over the leeward side, which was the case for the current event.
From a modeling point of view, resolving of such complex atmospheric patterns is a very challenging task since a high-resolution nested grid may often be required at both the origin areas of

Summary and Conclusions
We analyzed the most intense Saharan dust episode ever recorded at Finokalia station in Crete during its 14 years operational period (2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018). Analysis of the atmospheric processes during this episode suggests that the unique PM 10 concentrations recorded in Crete on 22 March 2018 were produced by a combination of meteorological processes taking place at various atmospheric scales. In general, the main driver for strong dust events in Mediterranean is atmospheric dynamics and this particular event is a typical example that combines several of the dynamical processes commonly involved in dust episodes: At synoptic scale, the weakening of the polar vortex allows descending of cold stratospheric air of high potential vorticity towards lower latitudes at Mediterranean and north Africa.
At mesoscale, the tropospheric folding favors cyclogenesis and intensifies the low-pressure systems along the surface baroclinic zones. The formation of mid-latitude cyclones at this area results in mobilization of dust particles from the Sahara and the scale and movement of these Mediterranean depressions determines also the dust pathways towards Europe.
At local scale, the vertical stability of the dusty air mass is disturbed when the plumes reach the orographic barriers of Mediterranean islands (e.g., Crete, Cyprus, Sicily) or the inland continental mountain ranges. Depending on the relevant heights between the mountain ridge and the elevated plumes, topography may act either as a stopping barrier that leads to increased dust concentrations at the windward ridge [22] or as a primary reason for the generation of downslope Foehn winds and downward mixing of dust over the leeward side, which was the case for the current event.
From a modeling point of view, resolving of such complex atmospheric patterns is a very challenging task since a high-resolution nested grid may often be required at both the origin areas of dust and at the arrival areas, several hundred kilometers away. Our WRF-Chem simulations clearly identify the importance of local scale dynamics on the surface PM 10 . As expected, the 1 × 1 km simulation outperforms the coarser 3 × 3 km, 12 × 12 km and 24 × 24 km runs in terms of surface dust concentrations. The major spatiotemporal and quantitative properties of the event are reproduced in satisfactory agreement with the observations, considering the complexity of the event. For example the maximum modeled and in situ PM 10 in Finokalia is 5.5 mg m −3 and 6.34 mg m −3 respectively. However, even at the fine grid scale of 1 × 1 km, certain limitations are evident regarding the model performance, mostly related to sub-grid unresolved topographic features and to the associated dynamics induced by terrain variability at the low-level tropospheric flow. In general, the event was even more spatially limited than shown in the model and affected only the eastern parts of Crete, leaving the west parts (Chania) practically unaffected.
In this work we explain the physical processes related to atmospheric dynamics during this particular event. Further analysis of the in situ measurements and of the specific properties of the aerosol particles will follow. Moreover, future work is needed to investigate the impact of climatic forcing to the atmospheric patterns generating such events. Following the consideration that stratospheric intrusions are the primary force for dust-storm generation, any possible future changes in polar vortex properties due to climate change may be relevant to changes in the frequency and strength of dust events in the Mediterranean.