WRF-Chem model simulations of a dust outbreak over the central Mediterranean and comparison with multi-sensor desert dust observations

In this study, the Weather Research and Forecasting model with online coupled chemistry (WRF-Chem) is applied to simulate an intense Saharan dust outbreak event that took place over the Mediterranean in May 2014. Comparison of a simulation using a physics-based desert dust emission scheme with a numerical experiment using a simplified (minimal) emission scheme is included to highlight the advantages of the former. The model was found to reproduce well the synoptic meteorological conditions driving the dust outbreak: an omega-like pressure configuration associated with a cyclogenesis in the Atlantic coasts of Spain. The model performances in reproducing the atmospheric desert dust load were evaluated using a multi-platform observational dataset of aerosol and desert dust properties, including optical properties from satellite and ground-based sun photometers and lidars, plus in situ particulate matter mass concentration (PM) data. This comparison allowed us to investigate the model ability in reproducing both the horizontal and the vertical displacement of the dust plume, as well as its evolution in time. The comparison with satellite (MODIS-Terra) and sun photometers (AERONET) showed that the model is able to reproduce well the horizontal field of the aerosol optical depth (AOD) and its evolution in time (temporal correlation coefficient with AERONET of 0.85). On the vertical scale, the comparison with lidar data at a single site (Rome, Italy) confirms that the desert dust advection occurs in several, superimposed “pulses” as simulated by the model. Crossanalysis of the modeled AOD and desert dust emission fluxes further allowed for the source regions of the observed plumes to be inferred. The vertical displacement of the modeled dust plume was in rather good agreement with the lidar soundings, with correlation coefficients among aerosol extinction profiles up to 1 and mean discrepancy of about 50 %. The model–measurement comparison for PM10 and PM2.5 showed a good temporal matching, although it revealed a marked overestimation of PM10 and PM2.5 (of the order of 70 % during the dust peak). For PM10, it was also possible to investigate the accordance between the modeland the measurement-based dust PM10, this confirming the model PM10 overestimation to be related to over-predicted dust mass up to a factor of 140 %. In all the model–measurement comparisons performed, the enhanced capabilities of the physics-based emission scheme with respect to its simplified, minimal version were evident and are documented.


2
The model-measurements comparison for PM 10 and PM 2.5 showed a good temporal matching, although it revealed a marked overestimation of PM 10 and PM 2.5 (of the order of 70% during the dust peak). For PM 10 , it was also possible to investigate the accordance between the model-and the measurements-based dust-PM 10 , this confirming the model PM 10 overestimation to be related to over-predicted dust mass up to a factor of 140%. In all the model-to-measurements comparisons performed, the enhanced capabilities of the physics-based emission scheme with respect to its simplified, minimal version were evident and 5 are documented.

Introduction
One of the main sources of uncertainty in our understanding of long-term climate variability is the role played by aerosols, since the related uncertainty greatly exceeds that of the other mechanisms combined all together (Houghton et al., 2001;IPCC, 10 2007). Among aerosols of natural origin, mineral dust is the foremost specie, comprising as much as 75% of the global aerosol mass burden, as estimated by satellite products (Ginoux et al., 2012). The role of mineral dust in the Earth system includes the interactions with other physical, chemical and biogeochemical processes at all scales (Shao et al., 2011b). It affects the earth's climate in many different ways, which are not completely understood and predictable, and influences the atmosphere-Earth balance, directly by scattering and absorbing short-and long-wave radiation with consequences on the net heating rates (e.g., 15 Alpert et al. 1989, Balkanski et al., 2007. The uncertainties in the direct radiative forcing are primarily attributed to the mineral aerosol shape (Kalashnikova andSokolik 2002, Haapanala et al., 2012), but also to their optical properties (Sokolik and Toon 1999;Bi et al., 2011) and their chemical composition (e.g., Claquin et al. 1998). In addition to these direct effects, aerosol indirectly affects the radiative balance by modifying cloud properties (e.g., Rosenfeld et al., 2001;Ghan and Schwartz, 2007;Wang et al., 2010;Karydis et al., 2011;Huang et al., 2014). 20 It has been estimated that about half of the global total natural dust emissions are generated in the Sahara Desert and its surroundings (Goudie, 2009;Huneeus et al., 2011;Ginoux et al., 2012;Shao et al., 2011b). Deep convection produced by the strong surface heating can uplift mineral dust particles for several kilometers into the free troposphere, where they are finally advected over large distances at the continental and intercontinental scales (Goudie and Middleton, 2001;Engelstaedter et al., 2006). Saharan dust is mainly transported along four trajectories patterns (D'almeida 1986;Shao 2011b). The largest fraction 25 (60%) of the dust loaded from Saharan sources remains in Africa, being transported and deposited in the Sahelian countries along the well-known "meningitis belt" (Molesworth et al., 2003). Another significant fraction (25%) is transported eastward across the Atlantic Ocean (e.g. Prospero and Mayol-Bracero, 2013;Yu et al, 2015), but a relevant (10%) Saharan dust amount is also carried across the Mediterranean Sea to Europe (Moulin et al., 1998, Israelevich et al., 2012 in episodic storms and/or following seasonal patterns Gobbi, 2004, Pey et al., 2013). The remaining 5% is 30 transported eastward to the middle east. During such outbreak events, mineral dust may be considered as the largest PM 10 3 source at urban and rural sites in the Mediterranean basin (Kaskaoutis et al., 2012;Pey et al., 2013;Salvador et al., 2014;Kabatas et al., 2014, Barnaba et al., 2016, contributing to a relevant percentage of the episodes of PM 10 daily limit exceedance (50 µg m -3 ) registered at these sites (Salvador et al., 2014;Barnaba et al., 2016), with peaks of contribution up to 80% of the total mass (Kaskaoutis et al., 2012).
During the year, the transport pathway of Saharan dust towards the Mediterranean is mainly determined by low-pressure 5 systems over the Atlantic or North Africa, high pressure over the Mediterranean region and/or high pressure at upper levels over Africa (Moulin et al., 1998;Barkan and Alpert, 2008;Querol et al., 2009;Pey et al., 2013, Salvador et al., 2014. Using Meteosat retrievals of dust optical depths, Moulin et al. (1998) showed that the northwards transport of dust follows a seasonal pattern, being eastward when associated with the Sharav cyclones (Alpert and Ziv, 1989), and directed toward the western Mediterranean basin from March to August, caused by the coupling between a Saharan low and a Libyan high or by a 10 cyclogenesis in the Atlantic coasts of Spain.
Modeling the transport of desert dust is receiving increasing attention from the scientific community, allowing to better ascertain its impact on radiation budget (Hsu et al., 1999), clouds (Bangert et al., 2011), as well as on air quality (Goudie and Middleton, 2001;Pey et al., 2013, Barnaba et al., 2016 and human health (e.g., Mallone et al. 2011, Stafoggia et al., 2016. Despite many improvements in characterizing dust source regions thanks to satellite products (Ginoux et al., 2012, Schepanski 15 et al., 2012, modeling dust emission and transport is still challenging due to the high uncertainties associated to the diffuse character of the emissions, re-suspension processes, the inherent complexity of aerosol chemistry and meteorological conditions, which strongly influence dust outbreaks and their spatio-temporal fields (e.g., Knippertz and Todd, 2012). This was evident in the intercomparison performed among 15 different global models in the framework of the Global Aerosol Model (AeroCom) initiative (Schulz et al., 2009) as well as in a recent intercomparison study among 9 European regional dust model 20 simulations (Basart et al., 2016).
Aim of this study is to evaluate the capability of the WRF-Chem model using a physical-based desert dust emission scheme to properly simulate an episode of mineral dust long-range transport occurred over the Central Mediterranean in May 2014.
This aim is pursued taking advantage of the operational aerosol and dust observations available from satellite and AERONET (Aerosol Robotic Network) sunphotometers (Holben et al., 1998) plus additional ground based observations carried on in 25 Central Italy within the EC-LIFE+ DIAPASON Project , Struckmeier et al., 2016Barnaba et al., 2016).
The desert dust event actually consisted of a series of dust plumes generated in the Northwest Sahara by strong winds associated to an omega-like circulation, characterized by a low-pressure system localized in the Atlantic coasts of Spain. Dust plumes were transported northward, resulting in an intense dust event over the Mediterranean (Aerosol Optical Depth, AOD, at 550 nm > 1) impacting mostly Italian and French sites, with maximum desert dust loads between May 21 and 23, 2014. 30 4 A preliminary study by Rizza et al. (2016) used WRF-Chem with a dust emission scheme based on a semiempirical dependence between the horizontal and the vertical dust fluxes (Marticorena and Bergametti, 1995), highlighting a large over-prediction of the AOD dust outbreak over Europe. In this work we test the use of a more advanced physics-based dust emission scheme proposed by Shao (2001, hereinafter denoted as S01), which explicitly considers the two major emission mechanisms for mineral dust, namely the saltation bombardment (Marticorena and Bergametti, 1995) and the aggregate disintegration, and a 5 refined four-classes texture soil type. To provide further insight into the advantages of the S01 scheme, in this study we also make a comparison with the model outcomes using its 'minimal' version (described by Shao et al., 2011a, and referred to as S11 in the following), in which the dust emission is independent from the sand particle size (results of this additional 'sensitivity' test are included in a specific Appendix).
The outline of the work is as follows. The setup of the WRF-Chem model used here is described in Section 2. Data and methods 10 used for the comparison with meteorological and aerosol fields are described in Sections 3. Results are shown and discussed in Section 4, in which we evaluate first the ability of the WRF-Chem model in reproducing the synoptic situation (4.1), and then the dust field in the horizontal (4.2) and vertical (4.3) scale exploiting the multi-sensor aerosol and desert dust observational dataset. Concluding remarks appear in Section 5.

15
2 The Wrf-Chem model WRF/Chem is a fully coupled online community model for the prediction and simulation of weather, dispersion, air quality, and regional climate (Grell et al., 2005). The chemistry model has been built to be consistent with the Weather Research and Forecasting (WRF, http://www2.mmm.ucar.edu/wrf/users/) modeling package. Possible applications of the current modeling system concern: (i) the prediction and simulation of weather, or regional and local climate; (ii) the release and transport of 20 constituents through coupled weather prediction/dispersion model simulations; (iii) the analysis of the full interaction of chemical species as well as particulate matter through coupled weather/dispersion/air quality model; (iv) the study of processes that are important for global climate change issues, including the aerosol direct and indirect forcing.

25
In this study the WRF-Chem version 3.6.1 was used. Figure 1 shows the model domain, which covers North Africa, Southern Europe and the western part of Asia, with 160 x 90 grid points centered at lat = 30.6° N, lon = 18.7° E. In the same figure the location of the six AERONET stations used in the analysis is reported. The horizontal grid spacing is 50 km for both directions with 40 vertical levels up to 50 hPa. The simulation lasted 10 days, starting on May 16, 0000 UTC. Boundary and initial conditions were extracted from NCAR/NCEP Final Analysis (FNL from GFS) (ds083.2), with 1-degree resolution, available 30 5 every 6 hours. An idealized vertical profile for each chemical species is provided to start the model simulation. This vertical profile is based upon northern hemispheric, mid-latitude, clean environment conditions.

Physical parameterizations
As summarized in Table 1, the following physical schemes are used. The Mellor-Yamada-Nakanishi and Niino (MYNN) 2.5 5 level turbulent kinetic energy (TKE) parameterization is used to describe the planetary boundary layer (Nakanishi et al., 2009).
The MM5 similarity scheme (Paulson, 1970) and the RUC Land Surface Model (Benjamin et al., 2004) are chosen to represent the surface layer physics and the land surface interaction. The Rapid Radiative Transfer Model (RRTMG) for both short-wave (ra_sw_physics = 4) and long-wave (ra_lw_physics = 4) radiation is used for the aerosol direct radiative effect (Mlawer et al., 1997). The Purdue Lin scheme (mp_physics = 2) is used for the treatment of the microphysics processes, with all 10 parameterization production terms based on Lin et al. (1983) and Rutledge and Hobbs (1984), with some modifications, including saturation adjustment following Tao et al. (1989) and ice sedimentation. This setup is compatible with the shortwave radiative feedbacks (or what is known as the 'direct effect'), which are included with the chemistry.
Aerosol optical properties are derived using the Maxwell-Garnett mixing rule (aer_op_opt = 2 within the model configuration, see Table 1) in its approximate parameterization, that is considering small spherical randomly distributed black carbon cores 30 in particle (Bohren and Huffman, 1983). 6 Aerosol removal processes include both dry and wet deposition. In particular, a dry deposition scheme accounting for gravitational settling and surface deposition is used (Wesely, 1989) to simulate the dry removal of desert dust, while a simple wet deposition scheme that considers rainout/washout in large-scale precipitation (Balkanski et al., 1993) is used for both seaspray and desert dust aerosols (Su and Fung, 2015). Note that in this scheme only non-parameterized (resolved) precipitation is active in the aerosol removal. 5

Dust emission parameterization
In any dust emission model, the basic parameters to be considered are: 1) the threshold friction velocity at which dust particles begin to move, 2) the horizontal and vertical sand dust fluxes. The emission of dust particles can be classified considering the 10 wind conditions at the surface. In particular, under strong wind conditions the surface wind-shear is the principal dynamic parameter and the dust emission is generally a function of the threshold friction velocity (see further details in Section 4.2).
Under these conditions two main dust emission mechanisms have been recognized: saltation bombardment (Marticorena and Bergametti, 1995) and aggregate disintegration (Shao, 2001). Another important mechanism is the direct aerodynamic lifting (Klose and Shao, 2012), which is effective when the lower troposphere is in the free-convective regime.

15
The WRF-Chem model (version 3.6.1) includes three alternative packages for mineral dust emission, two from the GOCART model ('DUST-GOCART' and 'DUST-GOCART/AFWA'), and a third ('DUSTUOC') from the University of Cologne. This latter is further divided into three emission parameterizations with a progressive level of simplification (Shao 2001(Shao , 2004(Shao , 2011a. A preliminary comparison between the DUST-GOCART/AFWA and the DUSTUOC-Shao 2001 (S01) schemes is discussed in Rizza et al. (2016), showing that the GOCART/AFWA emission scheme produces an important over-prediction 20 of the dust concentration. This is also in agreement with recent findings by Fountoukis et al., (2016). This may in part be explained by the fact that the AFWA scheme considers vertical dust flux only related to the clay content, while the S01 scheme considers a more realistic soil texture type.
In this work we have therefore preferred the S01 scheme, and tested its performances in comparison to observations. As mentioned above, results from the 'minimal' version of DUSTUOC-Shao et al., 2011a (here S11), are also included in the 25 Appendix A of the present study to highlight the advantages of the S01 scheme. In S11, the size distribution of the airborne sand and dust particles is only constrained by the minimally disturbed particles size distribution. The expression 'minimally disturbed' refers to the case in which the disturbance is so weak that the disintegration of aggregates almost does not occur (Shao, 2004).
Both S01 and S11 are based on a dust emission parameterization that considers explicitly the two major dust emission 30 mechanisms described above. In particular, the aggregate disintegration is modeled following the hypothesis that dust 7 aggregates fragment as they hit the surface. Both schemes may be considered as spectral emission schemes, because they are based on a size-resolved dust emission equation by supposing that particles are divided into n=4 particle-size intervals. In both cases the total dust flux (F) is expressed as an integral of the dust emission rate for particles of size d i by saltation of particles of size d s : where f(d i ,d s ) is the dust emission rate for particles of size d i generated by the saltation of particles of size d s . The difference between the S01 and S11 schemes is in the way the dust emission rate is calculated (see eq. 52 of S01 and eq. 34 of S11). The quantity p(d i,s ) can be regarded as a combination of two idealized particle size distributions, known as minimally disturbed particle size distribution, p m (d i,s ), and fully disturbed particle size distribution, p f (d i,s ), whose values are provided in a look-up table for each soil category. As p m and p f are functions of land surface properties, the soil data used in WRF-Chem play an 10 important role in dust emission simulation. In this study, the default soil categorization data set from the United States Geological Survey (USGS) with 5′ geographic resolution was selected. The effective soil texture is obtained from the USGS 12 classes considering only 4 types classes, and namely: sand, sandy clay loam, loam and clay. More details of this formulation can be found in Shao (2001), Shao et al. (2011a), Kang et al. (2011), and Su and Fung (2015). The chemistry and dust emission parameterizations adopted here are reported in Table 2. 15 3 Observational Dataset

Meteorological fields from NCEP/NCAR
Geopotential height maps at 500, 700 and 850 hPa are obtained using the daily mean composites of the NCEP/NCAR (National 20 Centre for Environmental Prediction/National Centre for Atmospheric Research) reanalysis (Kalnay et al., 1996). Composites (averages) of the daily-mean variables over several days are created from the NCEP/NCAR Reanalysis (http://www.esrl.noaa.gov/psd/data/).

25
The observation-based characterization of the aerosol field over the horizontal scale is made here by using both a network of sunphotometers located at multiple sites and measuring synchronously, and satellite retrievals capturing wider areas within a single passage. In particular, AERONET sunphotometers (Holben et al., 1998)

AERONET AOD dataset
AERONET is a federation of ground-based sunphotometers established by the USA NASA and currently led by NASA and the French CNRS. It includes nearly 1000 sunphotometers (CIMEL ®) spread worldwide, whose data are processed following the same aerosol retrieval procedures (Dubovik et al, 2000a, Dubovik et al, 2000b, Dubovik et al, 2006 and made available 5 in quasi real-time through the dedicated NASA portal (aeronet.gsfc.nasa.gov).
The main quantity measured by sunphotometers is the aerosol optical depth (AOD), an optical parameter quantifying the aerosol load in the whole atmospheric column. The AOD is unit-less, and represents the integral over altitude of the aerosol extinction coefficient (units of length -1 ).
In this study we use Level 2 (L2, i.e. cloud screened and quality assured) AOD measurements in the visible spectrum performed

MODIS AOD dataset
The MODerate-resolution Imaging Spectroradiometer (MODIS, Salomonson 1989) instrument flies on board the NASA TERRA and AQUA spacecrafts since December 1999 and May 2002, respectively. It has 36 wavelength bands spanning from the visible to the infrared, high spatial resolution, and near daily global coverage. TERRA, whose data are used here, overpasses the equator at 1030 LT. Aerosol characterization was and currently is a core MODIS mission (Kaufman et al., 1997) and the 25 AOD is still the most robust aerosol physical parameter derived from space. Two different approaches are used to retrieve the AOD from MODIS data. These are commonly referred to as 'Dark Target' (DT, Kaufman et al. 1997) and 'Deep Blue' (DB, Hsu et al., 2004). The algorithm at the basis of the 'Dark Target' approach is further differentiated when applied over ocean (Remer et al., 2005) or land (Levy et al. 2007a, b), and it is not suitable to be applied over bright surfaces (deserts, snow, sun glint). The 'Deep Blue' approach was developed to fill this gap (Hsu et al., 2004) and well complements the 'Dark Target' 30 retrievals. The most recent collection 6 (C006) of MODIS AOD data provides a single AOD product combining both the DT 9 and the DB AOD retrievals and, as considered the "best-of" product for most quantitative purposes (Levy et al., 2013), it was used in our study (in particular here we use the MODIS daily product MOD08_D3 v6).

Aerosol altitude-resolved view over Rome (Italy)
The characterization of the aerosol field over the vertical scale is made here employing continuous (h24) lidar/ceilometer 5 measurements performed in the CNR-ISAC Rome Atmospheric Supersite (CIRAS) in Rome -Tor Vergata, which hosts one of the six AERONET sites considered in this study (point A in Figure 1). The site is frequently affected by Saharan dust (e.g., Gobbi et al., 2004;, and lies just in the middle of the area impacted by the desert dust event under examination Barnaba et al., 2016). The Rome-Tor Vergata lidar and ceilometer measurements are therefore used here to evaluate the model capability to reproduce the dust plume vertical structure and its transport timing, as well as to provide 10 further insight into the model-vs-measurements AOD comparison.
In May 2014, the lidar and ceilometer measurements at CIRAS were part of a larger set of aerosol observations performed in the framework of the EC-LIFE+ project DIAPASON ('Desert-dust impact on air quality through model-predictions and advanced sensors observations', www.diapason-life.eu; more details on the project and relevant results can be found in Barnaba et al., 2016;Strukmeier et al., 2016). Lidar/ceilometer instruments characteristics and relevant dataset 15 used in this study are described hereafter.

Lidar datasets
The aerosol vertical profiles were collected by two different co-located instruments: a commercial CHM15K ceilometer (Lufft Mess und Regeltechnik GmbH, www.lufft.com), and a research-type lidar (ATLAS) developed at the ISAC-CNR laboratories.

20
The former is now part of an in-progress Italian network of such systems (the Automated Lidar-Ceilometer Network, ALICENET, www.alice-net.eu), which are conversely already widely employed in Germany, where the national meteorological service (DWD) operates over 50 of these instruments (e.g., Flentje et al. 2010a, Flentje et al., 2010bWiegner and Geiß, 2012). The CHM15K instrument uses a pulsed Nd:Yag laser source at 1064 nm with an output laser energy of about 8 µJ, a pulse repetition rate of 5-7 kHz and a vertical resolution of 15 m. Its configuration allows to sound the aerosol load in 25 the atmosphere in the range 150 m -15 km. As for all lidar systems, the signal in the lowermost atmospheric levels has to be corrected due to the incomplete superposition of the laser and the receiver field-of view (FOV). For this system the overlapping correction function is provided by the manufacturer, which determines the correction in the factory using a reference instrument.
The ATLAS system is a further miniaturization of a previous mobile, polarization-sensitive lidar system (VELIS) developed 30 by ISAC-CNR (Gobbi et al., 2000). ATLAS maintains most of the VELIS characteristics although it uses a different 1 kHz, 30 µJ/pulse laser source and reaches full overlap at approximately 500 m. As VELIS, ATLAS has two receiving channels, collecting respectively the light backscattered by particles in the parallel and perpendicular polarization planes with respect to the laser emitted one. Since spherical particles do not change the polarization plane of the incident light while non-spherical particles do, the comparison of the two lidar channels allows to detect the presence of non-spherical aerosols (as mineral particles) in the atmosphere (e.g., Gobbi 1998). An example of this capability is provided in Section 4.3.

5
Both the CHM15K and ATLAS are able to work unattended in continuous mode (h24), and their measurements are therefore used here to investigate the capability of the model to reproduce the dust plumes over Rome in terms of both temporal matching and vertical extent. In particular, we use both the (qualitative) range-corrected lidar signal (RCS), to a first approximation related to the aerosol amounts, and the (quantitative) aerosol extinction profiles from ATLAS (see Section 4.3). In fact, the inversion of the lidar RCS into aerosol optical properties (aerosol backscatter and extinction coefficients, β a and α a , 10 respectively) requires the employment of the backward solution of the Klett inversion algorithm (Klett, 1981) to the data. In addition to the estimation of the molecular backscatter and extinction coefficients (β m and α m , respectively, calculated from climatological monthly air density profiles), the solution requires two assumptions: a boundary value at a reference height z 0 where β a (z 0 ) = 0 (Rayleigh calibration) and a so called 'lidar ratio' (S a = α a /β a ). In our case, a calibration constant was derived applying the Rayleigh calibration to nighttime and cloud-free signals averaged over 1 h at 75 m height resolutions. For the 15 second assumption, we used an approach based on numerical simulations of aerosol scattering (e.g., Gobbi, 2001, Barnaba et al., 2004), and widely validated elsewhere (Gobbi et al. 2003;. In particular, in this study α a is computed using a functional relationship α a = α a (β a ) derived by Barnaba and Gobbi (2001) assuming non-spherical desert dust particles. The expected error on α a is of the order of 30%. This approach requires an iterative inversion technique to correct the backscatter signal for extinction losses until convergence in the integrated aerosol backscatter (IAB=∑ 0 zcal β a (z)) is reached.

20
The estimation of the aerosol extinction coefficient at altitudes below complete superposition of the laser and telescope FOV is obtained from a linear fit of the first two valid lidar points.

In situ PM 10 data
To complement the column-integrated and the vertically-resolved aerosol optical properties described above, for the Rome 25 site the observational dataset used to test the model also includes the standard particulate matter (PM) metrics regulated by the EU Air Quality legislation (i.e., the daily average PM 10 and PM 2.5 data). In this case hourly resolved measurements at the Rome-Castel di Guido site (about 15 km W of the city center) were used. These were collected using a SWAM dual channel instrument (FAI, Italy, http://www.fai-instruments.com/images/img/pdf%20eng/DOC401.PDF), providing mass concentration measurement on hourly basis thanks to a specific application of ß technology including information about 30 atmospheric mixing ratio. Relevant results are provided in Section 4.3.1.

Model capability to correctly reproduce the meteorology driving and associated to the dust event
Several authors evidenced that the northwards dust transport pathway from Sahara follows a seasonal pattern, changing from 5 the eastern to the western Mediterranean basin during spring and summer (Moulin et al., 1998;Engelstaedter et al, 2006).
As 'case study' representative of the springtime conditions, we selected an intense dust episode affecting the Central Mediterranean between 19 and 24 May 2014. This case corresponds to one of three different major cyclogenesis situations that are thought to be responsible for the northwards transport of Saharan dust toward the Mediterranean (e.g., Engelstaedter 10 et al., 2006), that is the cyclogenesis in the Atlantic coasts of Spain.
The synoptic analysis of the dust event is described using the NCEP/NCAR reanalysis (Kalnay et al., 1996). In Figure 2 we show the geopotential height in the middle (500 hPa) and lower (850 hPa

25
The ability of the model to reproduce the meteorology driving/associated with the dust event is evaluated in terms of geopotential height (Figure 3) during selected dates (21, 22 and 23 May). The geopotential field is obtained from the reanalysis at 700 hPa (panels a, c, e, in Figure 3) and compared to the corresponding WRF-Chem simulations (panels b, d, f) for the three selected dates. This quantity is important because the geopotential at 700 hPa gives an indication of the circulation pattern associated with the dust transport in the low-middle troposphere.

-Identification of desert dust source areas
As described above, an important condition for the existence of a dust source is the availability of fine grained material, which can be lifted from the ground when the surface wind speed exceeds a definite threshold. The threshold wind velocity depends on the surface roughness and grain soil size and in literature it is found to vary from about 6 to 9 ms -1 (e.g., Chomette 15 et al., 1999).
In order to evaluate the location of the dust sources that are directly connected with the investigated dust intrusion, in Figure   4 we superimpose the modeled AOD at 550 nm (shaded contours) and the total dust flux calculated with the S01 scheme (black The first dust plume (Figure 4a) was generated by easterly surface winds of approximately 20 m s -1 speed (black arrows) in the region marked by the ellipsoid S1. The peak value of about 100 µg m -2 s -1 was located at about (lat, lon) = (34° N, 8° W) that roughly correspond to the source area of Chott el Jerïd in Tunisia (Ginoux et al., 2012). This is a large endorheic Salt Lake (Chott), which becomes salt flats as it dries. The emission took place between 0600 UTC and 1600 UTC on May 18. To give 13 a reference for the following analysis, the model-based temporal evolution predicts this dust plume originated in S1 to reach the Rome area approximately at 0800 UTC on May 19 and above 2 km altitude.
The second dust plume (Figure 4b) was generated by westerly surface winds of approximately 20 m s -1 speed (black arrows) in the region marked by the ellipsoid S2. The peak value of about 100 µg m -2 s -1 was located at about (lat, lon) = (30° N, 5° W) that roughly corresponds to the source region of the Grand Erg Occidental Desert (Ginoux et al., 2012). The emission took 5 place between 1200 UTC and 1800 UTC on May 20. According to the model simulations, the core of this second dust plume reached Rome at about 2200 UTC on May 21, above 3 km altitude.
The third and most intense dust plume was generated on May 21. The source region is shown in Figure 4c and it is delimited by the ellipsoid S3 whose peak value of roughly 120 µg m -2 s -1 is located at ( A fourth, weaker dust impulse is produced on May 24. It was generated on May 24 at 1400 UTC on the region delimited by the ellipsoid S4 in Figure 4d, and reached Rome in the evening of the same day. Overall this result highlights that the investigated dust event was actually characterized by multiple, superimposed dust impulses, a pattern that will be confirmed 15 by the lidar record of desert dust profiles over Rome (Section 4).
All the identified source areas (S1, S2, S3 and S4) are located within a persistently active source region situated south of the Atlas Mountains and characterized by a system of ephemeral salt lakes that stays dry in summer, but receives some water in winter. This system may play an important role in modulating dust emissions (Engelstaedter et al., 2006;Salvador et al., 2014).
In Appendix A ( Figure A1), a view similar to that of Figure 4 is reported except that the S11 emission scheme is used. It is 20 evident that the two source regions evidenced in Figure 4 by the ellipsoids S2 and S3 are not present in the latter case, thus producing a dust emission pattern not compatible with the multiple dust plumes observed by lidar (Section 4).

Comparison of WRF-Chem and MODIS AOD over the Mediterranean
In Figure 5 we report the comparison between the AOD at 550 nm retrieved by MODIS-TERRA (left column) and that 25 predicted by WRF-Chem (S01 scheme, right column). In particular, each row of Figure

20
To complement the comparison with the satellite observations and better follow the temporal evolution of the AOD field over the Central Mediterranean, we used the AERONET measurements at the six sites shown in Figure 1. the measured AOD is well reproduced for the whole period. The corresponding results for the S11 scheme simulation are reported in Figure A3 in Appendix A. As expected from the previous comparison to MODIS AOD data, the results confirm the relevant S11 AOD under-prediction at the six AERONET sites. To quantify the model vs AERONET agreement, the mean bias and the temporal correlation coefficient for the six AERONET stations are reported in Table 3 (for both S01 and S11 5 cases). Overall, the average S01 AOD for the whole period for the six stations is 0.15 corresponding to a bias of -0.06 and an 85% of temporal correlation coefficient with respect to AOD measurements. For comparison, the S11 prediction produces a larger AOD negative bias (-0.12) and a lower correlation (71%).
This perspective shows that the maximum of a (first) dust plume reached southern Sardinia on 21 May, while in the easternmost sites (Ersa, Ispra, Rome) the maximum AOD is detected on 22 May. This result confirms and integrates the analysis of satellite 10 data, while the lidar view in Rome (next section) will further show that the latter maximum comes from the superposition of a second plume travelling above the first one.

Model capability to reproduce the desert dust plume vertical patterns;
The vertical evolution of the desert dust plume over Rome is revealed by the lidar and ceilometer measurements (see Section   15 2). Figure 7 shows the altitude (0-6 km) versus time (h24) cross section of the aerosol field as detected by both the CHM15k ceilometer ( Figure 7a) and the ATLAS lidar systems (Figure 7b), and forecasted by the model (Figure 7c) for the period 19-25 May 2014. In particular, Figure 7a shows the logarithm of the range-corrected signal (RCS) of the CHM15k system, which, to a first approximation, is proportional to the amount of particles loaded into the atmosphere. Figure 7b shows the particle depolarization ratio, derived from the two ATLAS receiving channels. This parameter is a 'marker' for the presence of non-20 spherical particles, and is therefore particularly suitable to follow the desert dust plume evolution in space (over the vertical scale) and time.
The ceilometer measurements (Figure 7a) also allows to follow the temporal evolution of the (aerosol-tracked) planetary boundary layer (PBL) in each day of the period considered. This can be identified by the green, bell-shaped areas reaching a maximum altitude of about 2 km (particularly evident on May 20,21,22). On May 19 and 23, the boundary layer signal is 25 somehow 'perturbed' by the presence of rain and clouds (red color in Figure 7a).
Above the PBL, elevated aerosol layers are clearly visible in the ceilometer trace. These have been highlighted by dotted, white oval shapes in Figure 7a. Although the (elastic) ceilometer signal of Figure 7a does not allow to discriminate the aerosol type, these elevated layers over Rome are typically associated to Saharan dust (e.g. Gobbi et al., 2004;Gobbi et al., 2013). To prove these layers are actually composed of mineral (non-spherical) particles, Figure 7b shows them to produce a 30 depolarization signal typical of long-range transported desert dust (volume depolarization > 8%). To facilitate the spatiotemporal comparison between the two lidars measurements and the model outcome, the same dust-identification oval shapes have been reported in each panel of Figure 7.
Overall, the lidar measurements in the period 19-25 May reveal the desert dust advection to occur in several, superimposed desert-dust plumes, thus confirming the 'pulsed' nature of this event simulated by the S01 configuration (Figure 7c, see also Section 4.2.1). In most of the cases both measurements and simulation show the desert dust to travel above the PBL and then 5 to descend and mix with the local aerosols within it. In particular, four main different desert-dust plumes, identified to originate over different source regions in Section 4.2.1, can be detected in the lidar/ceilometers records (oval shapes in Figure 7a, b). A first plume arrives over Rome on May 19 (although the presence of clouds prevents from establishing the exact time of its arrival over the city); then it progressively descends and is firstly detected at the ground on May 20. This is compatible with the modeled plume originated in the source area S1 in Figure 4. A second plume, compatible with the one originated by the 10 model in S2 (Figure 4), reaches the atmosphere near Rome above 2 km height in the evening of May 21 and then descends towards the ground. This plume superimposes to the previous one, and to a third major plume arriving in the afternoon of May 22 and extending from the ground up to 6 km (this is the major plume the model identifies to originate in the source region S3, Figure 4). The mixing of the three plumes is observed down to the ground until at least May 24. As predicted by the model, a fourth and weaker pulse arrives aloft in the night between 24 and 25 May and superimposes to the previous ones until the end 15 of the analyzed period.
Although qualitatively, lidar measurements at high vertical and temporal resolution allow to evaluate the model capability to reproduce the desert dust plume vertical patterns and timing. In particular, the comparison of Figures 7a, b with Figure 7c shows the model is able to reproduce well the 'pulsed' pattern of this desert dust advection, reproducing pretty well both its timing and vertical extent. Some differences are found with the timing and vertical location of the fourth plume, which the 20 measurements indicate to arrive at about 2 km around noon of May 24 and the model predicts at lower altitudes and with some hours of delay. Nevertheless, the model is still capable to reproduce the second peak of this fourth plume observed by the lidar systems aloft (above 3 km).
A more quantitative validation of the vertically-resolved model output is provided in Figure 8, in which the ATLAS lidar range-corrected signal (RCS) is inverted to derive the aerosol extinction coefficient (see Section 3.3.1). For this purpose, night-25 time/early-morning profiles have been selected to improve the signal-to-noise ratio in the measurements and thus facilitate the lidar signal inversion. Figure 8 shows that the model mostly reproduces the general vertical pattern of the desert dust plume, with a double layer structure clearly evident on May 22 (Figure 8a, b). The elevated layer is likely uniquely composed of desert dust particles (as revealed by the lidar depolarization trace in Figure 7), while in the PBL aerosol layer, desert dust is mixed with (spherical) particles of local origin. Overall, for this date, which corresponds to the maximum desert dust load over Rome 30 (see also Figure 6b), the model reproduces pretty well the associated aerosol extinction along the vertical profile (Figure 8a,b), with an estimated Normalized Mean Bias (NMB) with respect to the lidar ranging from -50% of the midnight profile ( Figure   8a) to +50% of the early morning one (Figure 8b).
As expected, the model is however unable to reproduce the 'fine' structure of the desert dust plume (minor thin layers within the main ones) revealed by the high-vertical resolution of the lidar trace. This translates into a moderate correlation (R ~ 0.6, 0.7) between the modeled and the observed aerosol extinction profiles. This vertically-resolved perspective also provides 5 further insight in the good matching between the model and the measured AOD over Rome on May 22 (Figure 6b), showing that it partly derives from 'compensation effects' between a modeled underestimated aerosol extinction in the lowermost levels, and a model overestimation aloft. Similarly, the vertically-resolved comparison for May 25 (Figure 8c, d) allows us to better understand why the model underestimates the AOD (Figure 6b) in such lower desert dust loads. In fact, although the model is still able to reproduce the shape of the aerosol profile (moderate-to-excellent correlation in Figure 8c, d), in this case it clearly 10 under-predicts the aerosol extinction (NMB of about -60%). This view, and particularly the steeper decrease of the aerosol extinction with height in the lowermost levels, also points to some underestimation of the PBL height by the model. The same model vs. lidar quantitative comparison employing the S11 simulation is reported in Figure A4, once again highlighting the worse performances of this configuration with respect to the S01 one. In fact, in this case a very low aerosol extinction coefficients is associated to the desert dust plume, and the vertical structure observed by lidar is completely lost in this 15 simulation.

Comparison to ground-level PM values
A further quantitative evaluation of the model ability to reproduce the observed aerosol/dust load is given in Figure 9, where the aerosol mass concentration predictions at the particular vertical level coincident with the ground are compared to in-situ 20 measured PM 10 and PM 2.5 data, these being the metrics regulated by the European Air Quality Directive 2008/50/EC (EC, 2008). In particular, in Figure 9 model-simulated (red curve) and in situ-measured (blue curves) PM 10  In fact, the model predicts desert dust to contribute up to 90% of the total PM 10 In quantitative terms, the comparison to the observed desert dust mass concentration reveals that the model overestimates the daily-average desert-dust-PM 10 by a factor of 2-3, as the observed desert-dust contributes to less than 50% of the (total aerosol) daily PM 10 .

5
As the comparison of the aerosol extinction profiles does not show a similar model overestimation in the lowermost levels (e.g., Figure 8), this result points to a probable misrepresentation of the desert dust size distribution within the model, with over-predicted role of large particles (impacting total aerosol mass more than aerosol extinction coefficient at 550 nm).
Furthermore, we speculate that an important factor in determining the poor model performance in predicting the PM fields could be related to an inadequate representation of the desert dust wet removal. In fact, the wet removal by convection-driven 10 (parameterized) precipitation is neglected in the simulation (Section 2.1.2) while it contributes to about 40-50% of the total rainfall simulated by the model (see Section 2.1.2). Actually, the event under investigation was associated to precipitation both during the northward advection of desert dust (particularly over the Mediterranean Sea between Sardinia and Central Italy, not shown) and during its transit over the Rome observational site (rain detected on May 19 and May 23, as visible from Figure   6a). This aspect would certainly merit further investigation by including in the model set up a new dust wet deposition scheme 15 (for example the one recently developed by Tsarpalis et al. (2017), which explicitly considers the rate of dust scavenged by precipitation inside and below the cloud, based on the parameterization proposed by Seinfeld and Pandis, 1998). This further analysis is however beyond the purpose of the current study.

20
This study evaluates the performances of a physics-based desert-dust emission scheme within the WRF-Chem model. In particular, we used the physics-based dust emissions scheme proposed by Shao et al., 2001 (here S01) to simulate an important dust outbreak occurred over the Central Mediterranean regions in the period 19-24 May 2014, and we compared the simulations to a large set of aerosol/desert-dust observations. To highlight the advantages of using a detailed physically-based scheme, the performances of the S01 configuration were also compared with the outcome of its simplified (minimal) version (Shao et al., 25 2011, here S11), in which the emission is independent from the sand particle size. In all the comparisons S01 outperformed the S11 scheme as described in the text and documented in Appendix A.
In the case study considered here, the intrusion of mineral dust in the Mediterranean was associated to a synoptic-scale omega-like pressure configuration with a cyclogenesis in the Atlantic coasts of Spain. In fact, the cyclone was responsible for strong westerly Atlantic winds affecting the northern Sahara, while the northward transport was made possible by 30 southwesterly currents on the west side of the ridge associated with the omega-like pattern. In general, the synoptic conditions for the geopotential height at 700 hPa were well reproduced by WRF-Chem. This allowed us to simulate with good confidence the path of the dust during the northward intrusion.
A first comparison between the modeled aerosol optical depth (AOD) field and source emission functions allowed us to identify four different source regions of desert dust for the investigated northward intrusion. In particular, we recognized a persistently active source region located south of the Atlas Mountains and between Algeria and Tunisia, confirming some 5 recent findings from Ginoux (2012). This region is in fact characterized by a system of ephemeral salt lakes that stays dry in summer, but receives some water in winter playing an important role in modulating the dust emissions.
A multi-platform observational dataset (including satellite, AERONET, lidar and in situ PM data) was used to test the ability of the model to reproduce the spatiotemporal pattern of the dust intrusion, which was found to be composed by several, superimposed, time-shifted dust-pulses.

10
On the horizontal scale, the comparison between the modeled AOD field with the corresponding satellite retrievals (MODIS-TERRA) showed the WRF-Chem simulation to satisfactorily resolve the arrival, the temporal evolution and the extent of the plumes over the Central Mediterranean. Results also showed a good agreement between the modeled AOD and the one measured from the ground at six selected AERONET sites in the Mediterranean region (correlation coefficients, R = 0.85).
The combined analysis of AERONET data and simulations showed that the first dust intrusion occurred on 21 May reaching 15 southern Sardinia, the second and most intense dust plume occurred on 22 May, penetrating up to northern Corsica and Central Italy. This result confirms and complements the analysis from the polar satellite data (TERRA), which is necessarily limited in time.
The characterization of the aerosol field over the vertical scale was made here by employing continuous (h24) lidar/ceilometer measurements performed in a single observational point that lies just in the middle of the area investigated 20 (Rome, Italy), and was therefore expected to be particularly suitable to evaluate the model capability to reproduce the dust plume vertical extent and its transport timing. Overall, the lidar measurements in the period 19-25 May clearly highlighted the 'pulsed' nature of the event examined. In most of the cases the desert dust is shown to arrive above the PBL and then to descend and mix with the local aerosols within it. The comparison with lidar measurements also highlighted that the good matching between the model and measured AOD comes from a rather good reproduction of the aerosol extinction coefficient along the 25 profile (Normalized Mean Bias, NMB, of about 50%), with the best performances in terms of aerosol optical properties obtained during the maximum of the dust event (minimum bias between the modeled and the measured aerosol extinction coefficients). During the weaker desert dust conditions registered at the end of the event, the model is shown to reproduce well the shape of the observed vertical extinction profile (correlation coefficient R up to 1), although with a marked underestimation (NMB of about -60%). When the model-measurement comparison was done at ground-level in terms of aerosol mass (PM 2.5 30 and PM 10 data are used for this purpose), a tendency to overestimate the desert dust aerosol mass was conversely revealed.
Such an overestimation reaches 70% for PM 10 and PM 2.5 , respectively during the dust peak, reduced to 10-60% in weak-dust or no-dust conditions (before May 20 and after May 23). For PM 10 it was possible to show that the total mass overestimation is driven by an overestimated dust contribution of the order of 140%. This result points to a possible over-prediction of the number of large dust particles by the model (affecting dust mass more than optical properties). Additionally, we speculate that at least part of the model PM 2.5 and PM 10 overestimation might be related to the simplified wet removal scheme adopted, which only considers non-convective (resolved) precipitation as active in the desert dust removal processes. In fact, particularly in 5 the central phase of the dust event recorded in Rome, convective precipitation was registered over the Central Mediterranean, between Sardinian and the Italian Tyrrhenian coast. This aspect deserves further investigation and will be addressed in a future work together with the tuning of the several model internal parameters that characterize this kind of size-resolved dust fluxes.

10
Analyses and visualizations used in this paper were produced with the Giovanni online data system, developed and maintained by the NASA GES DISC (http://giovanni.gsfc.nasa.gov/). We also acknowledge the MODIS and AIRS mission scientists and associated NASA personnel for the production of the data used in this research effort. We thank the principal investigators and their staff for establishing and maintaining the AERONET sites, the data from which have been used in this study. The aerosol and desert dust observations in Rome (Italy) employed in this work were supported by the EC-LIFE+ DIAPASON project 15 (LIFE+10 ENV/IT/391).   values obtained from measurements have been reported (dashed blue line) see text for details); (d) as in (b) but for but for daily-averaged data.
10 Figure A1: As in Figure 4 except that the S11 rather than the S01 configuration is used Figure A2: As in Figure 5 (right panel) except that the S11 rather than the S01 configuration is used 15 Figure A3: As in Figure 6 (right panel) except that the S11 rather than the S01 configuration is used Figure A4: As in Figure 8 except that the S11 rather than the S01 configuration is used 20 33 Tables with Table captions