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

Abstract

A record-breaking dust episode took place in Crete on 22 March 2018. The event was characterized by surface concentrations exceeding 1 mg m⁻³ for a period of 4–7 h, reaching record values higher than 6 mg m⁻³ 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 PM₁₀ 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.

1. 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₁₀ 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₁₀ 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.

2. 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/).

2.1. Modeling

We use the limited area model WRF-Chem in a two-way nesting configuration in order to describe the physical processes at various atmospheric scales. The grid-space for the external domain is set at 12 × 12 km. An intermediate domain at 3 × 3 km is used over Crete and the surrounding region and the inner grid is set at 1 × 1 km covering the entire island of Crete. An additional sensitivity run is also performed with a 24 × 24 km domain covering the same region as the 12 × 12 km domain. We initialize the main run at 21 March 2018, 0000 UTC and run it until 23 March 0000 UTC (T + 60). The internal 3 × 3 km and 1 × 1 km domains are activated for the period 21 March 1800 UTC–22 March 1800 UTC and 22 March 1200–1800 UTC respectively. The model is driven by initial and boundary conditions from the NCEP final analysis (FNL) dataset at 1° × 1° grid space and the sea surface temperature is also the NCEP SST at 1° × 1°. The planetary boundary layer scheme is the Mellor-Yamada-Janjic and the land properties are parameterized by the Noah Land Surface Model [29,30]. Surface heat and moisture fluxes are computed by the Monin-Obukhov-Janjic scheme [31,32]. Cloud convection is parameterized with the Grell and Devenyi [33] scheme for the two outer model grids and convective parameterization is switched off for the finest 1 × 1 km grid. The Lin et al., 1983 [34] scheme is used for the calculation of cloud microphysical properties. Radiative transfer is computed by the Rapid Radiative Transfer Model [35] for both short-wave and long-wave radiation taking also into account the effects of dust particles in radiative transfer.

The emission of dust is computed with the Air Force Weather Agency (AFWA) scheme [36,37] that follows Marticorena and Bergametti 1995 (MB95) [38] for the computation of vertical dust fluxes based on wind friction velocity, soil type and surface moisture. The MB95 dust emission scheme has been widely used and evaluated in global and limited area dust models e.g., DEAD: [39]; GEOS-Chem: [40]; SKIRON: [41,42]; RAMS: [8,18,19]. WRF-Chem includes five dust size bins with effective radii of 0.73, 1.4, 2.4, 4.5 and 8 μm. The modeled dust Aerosol Optical Depth (AOD) at 550 nm (${\tau }_{550}$) is calculated using the corresponding column mass load of dust ($M$) and the extinction efficiencies at 550 nm ($Q_{ext550}$) of the five WRF sub-bins, as in Equation (1):

$${\tau }_{550}={\displaystyle \sum }_1^{50}\frac{3}{4{\rho }_n{r}_{ef{f}_n}}{M}_n{Q}_{ext550}{}_n$$

(1)

where $\rho$ is the particle mass density, $r_{eff}$ is the effective radius and $n$ denotes the corresponding WRF sub-bin. The extinction efficiencies at 550 nm (

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_{eff}$) and the geometric standard deviation (${\sigma }_g$) of the sub-bin lognormal size distributions.

${\mathit{r}}_{\mathit{m}\mathit{i}\mathit{n}}-{\mathit{r}}_{\mathit{m}\mathit{a}\mathit{x}}\left(\mathsf{\mu }\mathbf{m}\right)$	${\mathit{r}}_{\mathit{e}\mathit{f}\mathit{f}}\left(\mathsf{\mu }\mathbf{m}\right)$	${\mathit{\sigma }}_{\mathit{g}}$	$\mathit{\rho }$ (kg m−3)	${\mathit{Q}}_{\mathit{e}\mathit{x}\mathit{t}550}$
0–0.1	0.73	2	2500	2.79851
1.0–1.8	1.4	2	2650	2.26168
1.8–3.0	2.4	2	2650	2.24124
3.6–6.0	4.5	2	2650	2.14511
6.0–10.0	8.0	2	2650	2.09845

) 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.

2.2. Remote Sensing

2.2.1. PollyXT LiDAR

The PollyXT LiDAR [47] is located in Finokalia (35°20′15.4″ N, 25°40′09.8″ E, 238 m ASL). The LiDAR is equipped with three elastic channels at 355 nm, 532 nm and 1064 nm, two rotational-vibrational 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.

2.2.2. 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. Dust RGB channel enhancement. The “night-mode” flag is raised when solar zenith angle exceeds 85°.

Component	Scaling Range	Gama Factor
Red	$\{\begin{array}{c}1K\left(day\right)\\ 2K\left(night\right)\end{array}\langle \left(B{T}_{12}-B{T}_{10.8}\right)\rangle -4K$	1.0
Green	$\{\begin{array}{c}15K\left(day\right)\\ 9K\left(night\right)\end{array}\langle \left(B{T}_{10.8}-B{T}_{8.7}\right)\rangle 0K$	$\{\begin{array}{c}3\left(day\right)\\ 2.25\left(night\right)\end{array}$
Blue	$\{\begin{array}{c}289K\left(day\right)\\ 293K\left(night\right)\end{array}\langle B{T}_{10.8}\rangle \{\begin{array}{c}261K\left(day\right)\\ 243K\left(night\right)\end{array}$	1.0

. Each of the components is scaled using the following equation:

$$C_i^{scaled}=\min \left(X\right)+\delta X\times {\left(\frac{C_i-\min \left(C_i\right)}{\delta C_i}\right)}^{\frac{1}{{\gamma }_i}}$$

(2)

where $C_i$ is the value for red, green and blue components; $\delta C_i=\max \left(C_i\right)-\min \left(C_i\right)$ for each component; $X=\left[4,254\right]$ is the color range; $\delta X=\max \left(X\right)-\min \left(X\right)$ is the number of colors; and ${\gamma }_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.

2.2.3. MODIS

Columnar AOD observations are obtained from the MODerate resolution Imaging Spectroradiometer (MODIS) on-board the NASA’s polar satellites Terra and Aqua since 2000 and 2002, respectively [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).

2.3. 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⁻³ (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₁₀ 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].

3. Results

From 1993 to March 2018, the maximum surface concentration of PM₁₀ ever observed at Finokalia was 2793 μg m⁻³ (recorded on 24 February 2006 with a time step of 5 min). Between 2004 and 2017, only six days with 5 min PM₁₀ concentrations larger than 1 mg m⁻³ have been recorded. The mean AERosol Robotic network (AERONET), [61] aerosol optical depth (AOD) in Crete at 500 nm, was 0.19 ± 0.10 (1σ) with maximum values during April and May (0.24 ± 0.15). For the particular months the minimum Angstrom exponents are observed, with averaged value equal to 0.86.

On 22 March, Crete was affected by a particularly severe North African dust event. The advected dust plume is evident at the region between Africa coast and Crete in the daily composites of MODIS visible channel (Figure 1a,b). In the geostationary Dust RGB images of MSG-SEVIRI (Figure 1c,d) the dust plume is indicated by dark magenta colors at 1200 UTC and 1500 UTC 22 March 2018. The presence of dust was initially detected in the western parts of the island (Chania, Greece), where the phenomenon peaked at noon, with maximum recorded PM₁₀ of 0.5 mg m⁻³. In the central and eastern areas (Heraklion and Finokalia), the dust concentrations also begun to increase during the morning hours of 22 March 2018. At these areas, the dust concentration increased rapidly so that at 1500 UTC in the afternoon it reached a record value of 4730 μg m⁻³ in Heraklion. The concentration of dust at the easternmost station (Finokalia), increased sharply in the afternoon, and at 1720 UTC a new record value for Greece was again recorded at 6340 μg m⁻³, exceeding the 5000 μg m⁻³ upper limit of the tested range reported by the instrument’s manufacturer. The average daily values recorded at Chania, Heraklion and Finokalia were 206 μg m⁻³, 1125 μg m⁻³ and 850 μg m⁻³ respectively.

3.1. 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⁻¹ (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⁻⁶ m² s⁻¹ K kg⁻¹). 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⁻¹). 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 at 850 hPa on 22 March 0600 UTC. The minimum temperatures and geopotential heights at this stage of development are −34 °C, 5400 m at 500 hPa (not shown) and −3 °C, 1340 m at 850 hPa respectively (Figure 3a). Moreover, as seen at the 850 hPa map (Figure 3a), a strong baroclinic zone extends from about 20° N, 12° E over the Sahara Desert up to the west coastline of Greece with a horizontal West-East temperature gradient of more than 15 °C between 15°–20° E. This combination of atmospheric dynamics (i.e., tropospheric folding, upper level PV-θ anomaly, baroclinic trough propagation, lower troposphere baroclinicity) is responsible for the formation and intensification of East Mediterranean cyclones and for the generation of dust storms [64,65,66,67].

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² s⁻²) 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⁻³ at certain frontal areas and the low-level convergence inside the warm sector actually traps the dust particles ahead of the cold front.

3.2. 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⁻² inside the core of the system (i.e., Figure 6b). The dusty Khamsin flow affects Crete between 1200–1800 UTC, 22 March 2018.

3.2.1. 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⁻³ 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⁻³. The wind speeds reach up to 30 m s⁻¹ 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).

3.2.2. 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⁻³. 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⁻³ 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.

3.2.3. 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⁻³. 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⁻³. 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⁻²) behind the front (after 1700 UTC). However, the downslope winds are underestimated by almost 10 m s⁻¹ 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⁻¹ 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 coarse-mode dominated, thus the mass-specific extinction coefficient value of $k_{ext}=0.52m^2{\mathrm{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⁻³ while the modeled profile quickly declines to values less than 1000 μg m⁻³ 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 PM₁₀ 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₁₀ is calculated from the five dust bins as: PM₁₀ = 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₁₀ concentrations exceeding 6 mg m⁻³ during 1400–1800 UTC. At the same time, PM₁₀ at Chania station that is located 140 km to the west was below 0.2 mg m⁻³. 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₁₀ 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₁₀ in Finokalia reaches 5.5 mg m⁻³ one hour earlier than the actual recorded peak of 6.34 mg m⁻³. However the model clearly underestimates PM₁₀ in Heraklion since although the simulation produces very high values exceeding 2.5 mg m⁻³ at this station, the actual concentrations were actually two times higher exceeding 4 mg m⁻³. In Chania station, the model fails to predict the very low recorded PM₁₀ 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₁₀ 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.

4. 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–2018). Analysis of the atmospheric processes during this episode suggests that the unique PM₁₀ 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₁₀. 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₁₀ in Finokalia is 5.5 mg m⁻³ and 6.34 mg m⁻³ 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.