Interactive comment on “ Atmospheric dust modeling from meso to global scales with the online NMMB / BSC-Dust model – Part 1 : Model description , annual simulations and evaluation ” by C

Response: From this general comment, we understand that the referee defines “offline” as a simulation not allowing for dust radiative feedbacks. The NMMB/BSC-Dust is a coupled “online” model. An “online” model does not necessarily imply that feedbacks between radiation and aerosol are allowed. For the definition of “online” and “offline” models we refer to Grell et al. (2005): Until recently, the chemical processes in air quality modeling systems were usually treated independently of the meteorological model, i.e., “offline”, except that the transport was driven by output from a meteorological model, typically available once or twice per hour. In “online” systems the air quality component of the model is fully consistent with the meteorological component; both components use the same transport scheme (mass and scalar preserving), the same grid (horizontal and vertical components), and the same physics schemes for subgrid-scale transport. The components also use the same timestep, hence no temporal interpolation is needed. This is exactly what the NMMB/BSC-Dust is about. It is an online model and the simulations shown in the paper are done in an online mode.


Introduction
Dust aerosol particles are produced by wind erosion of arid and semi-arid surfaces.The major sources of contemporary mineral dust production are found on the desert regions of the Northern Hemisphere, in the broad dust belt that extends from the Eastern Subtropical Atlantic eastwards through the Sahara Desert to Arabia and Southwest Asia (Goudie and Middleton, 2001;Prospero et al., 2002;Luo et al., 2004;Badarinath et al., 2010).Other significant large-scale sources can be found within the Australian, North and South American, and South African deserts (Formenti et al., 2001;Washington et al., 2003).Model estimates of the amount of dust exported annually are uncertain and range from 1000 to 2150 Tg yr −1 with the largest contributions emitted from the North African (50-70 %) and Asian deserts (10-25 %).
C. Pérez et al.: Dust modeling from meso to global scales -Part 1 Atmospheric burden estimates range from 8 to 36 Tg, an uncertainty factor that exceeds 4 (Zender et al., 2004) .
Strong dust events can reduce visibility to near zero at source regions.While visibility improves downstream of sources, dust particles are mixed vertically, reaching up to several kilometres, from where they are carried over distances of thousands of kilometers by strong winds aloft.With the possible exception of sea-salt aerosol, dust is globally the most abundant of all aerosol species (IPCC, 2001) and is the dominating component of atmospheric aerosol over large areas of the Earth.Dust impacts, interactions and feedbacks within the Earth System span a wide range of spatial and temporal scales.On the short term and close to sources, dust storms represent a serious hazard to health, property, environment and economy.At longer time scales, atmospheric dust affects the climate through direct and indirect radiative effects, and ocean biogeochemistry which in turn may influence the carbon cycle and the climate itself (e.g., Miller and Tegen, 1998;Jickells et al., 2005;Mahowald et al., 2009).
Dust models mainly predict dust emission, transport within the atmosphere and deposition.In the last two decades, several dust cycle models have been developed and coupled online or offline with short and medium-range weather forecast models or with climate models.Dust forecast models such as the regional BSC-DREAM (Nickovic et al., 2001;Pérez et al., 2006b), SKIRON (Kallos et al., 2006) and CHIMERE-Dust (Bessagnet et al., 2004) and the global NAAPS model (Westphal et al., 2009) have focused on the accurate representation of the short-term variability of the dust to provide air quality/visibility forecasts up to 3 to 5 days.Studies using dust climate models (e.g., GISS ModelE; Miller et al., 2006) have mainly focused on the accurate climatological representation of the seasonally dependent dust cycle as well as the study of the dust radiative forcing and its effects on climate at the global scale.
The inclusion of the radiative effects of dust and other aerosols is more unusual in regional and global weather forecast models.Studies have suggested that the inclusion of mineral dust radiative effects can improve the radiative balance of short and medium-range forecast models and thus can increase the overall accuracy of the weather prediction itself (e.g., Pérez et al., 2006b;Helmert et al., 2007).Recently, the ECMWF IFS medium-range forecast model has included a prognostic representation of aerosols including data assimilation of aerosol-related observations in a fully interactive way (Morcrette et al., 2009).
One main focus of dust modeling research is upon the dust emission process which is highly sensitive to meteorological, surface and soil features.While dust distribution and dust effects are important at global scales, dust emission is a threshold, sporadic and spatially heterogeneous phenomenon (Laurent et al., 2009) that is locally controlled on small spatial and temporal scales (Tegen and Schepanski, 2009).Therefore, predicting the magnitude and the spatio-temporal patterns of dust emission with sufficient accuracy remains among the crucial challenges for global and regional dust models.Current theoretical knowledge predicts the vertical dust flux in models if the required input parameters -surface, soil and meteorological features -are accurately determined (Laurent et al., 2009).However, the application of complex emission schemes in global and to a lesser extent regional models is hampered by a lack or strong uncertainty of the required input data at the scales of application as well as the inaccuracies of the driving meteorological/climate model.
Global models have assumed varying degrees of simplification in the dust emission schemes as a function of the availability and accuracy of the input data (e.g., Ginoux et al., 2001;Tegen et al., 2002;Zender et al., 2003a;Miller et al., 2006).Although simplifications are also applied to regional models, some of them include specific surface and soil datasets developed for the most well-known sources regions (e.g., Chatenet et al., 1996;Callot et al., 2000) and complex dust emission schemes.In most cases the emission is tuned to match quantitative dust observations that are mainly available far away from sources.
In this paper, we present and evaluate the NMMB/BSC-Dust, a new online multi-scale atmospheric dust model designed and developed at the Barcelona Supercomputing Center (BSC) in collaboration with NOAA/National Centers for Environmental Prediction (NCEP), NASA Goddard Institute for Space Studies and the International Research Institute for Climate and Society (IRI).The dust model is embedded into the Non-hydrostatic Multiscale Model NMMB developed at NCEP (Janjic, 2005;Janjic and Black, 2007;Janjic et al., 2011) and is intended to provide short to medium-range dust forecasts for both regional and global domains.We also expect the model to become a useful research tool that will improve our understanding of the dust cycle by bridging the gap among the multiple scales involved.In this sense, for example, the model will provide a common modeling framework to simulate dust emission at the local scale with very high resolution and explicit convection, and at the global scale with lower resolution and parameterized convection.Also, these developments represent the first step towards a unified multiscale chemical-weather prediction system (Jorba et al., 2011).
Section 2 summarizes the main characteristics of the NCEP/NMMB model.Section 3 describes the BSC-Dust model, including dust source estimation, emission scheme, dry deposition, wet scavenging, convective mixing and interaction with radiation.In Sect.4, we assess the regional and global configurations of the model.For the regional, we evaluate the daily variability of the dust aerosol optical depth (AOD) in 2006 for a domain covering Northern Africa, Middle East and Europe.The global configuration of the model is compared against the AEROCOM (Textor et al., 2006) dust benchmark dataset (Huneeus et al., 2011), making use of the tools developed at the Laboratoire du Climat et Sciences de l'Environnement in the framework of the AE-ROCOM project.In the global case, we compare monthly and annual means of variables that are related to the direct radiative effect, the dust impact on the ocean biogeochemical cycle and air quality, i.e. aerosol optical properties, dust deposition and surface concentration.A companion paper (Haustein et al., 2011) evaluates and analyses the behavior of the regional configuration of the model during two experimental campaigns in Northern Africa: SAMUM I (Heintzenberg, 2009) and BODEX (Washington et al., 2006).

The NCEP non-hydrostatic multiscale model
The Non-hydrostatic Multiscale Model NMMB (Janjic, 2005;Janjic and Black, 2007;Janjic et al., 2011) is a new unified atmospheric model for a broad range of spatial and temporal scales.Its unified non-hydrostatic dynamical core allows for regional and global simulations.The NMMB has been developed within the Earth System Modeling Framework (ESMF) at the National Centers for Environmental Prediction (NCEP) following the general modeling philosophy of the NCEP regional WRF Non-hydrostatic Mesoscale Model (WRF-NMM) (Janjic et al., 2001;Janjic, 2003) which was operationally used at NCEP as the regional North American Mesoscale (NAM) model.The regional NMMB became the NCEP NAM in october 2011.
The numerical schemes used in the model were designed following the principles set up in Janjic (1977Janjic ( , 1979Janjic ( , 1984Janjic ( , 2003)); Janjic et al. (2001) and Janjic et al. (2011).Isotropic horizontal finite volume differencing is employed so a variety of basic and derived dynamical and quadratic quantities are conserved.Among these, the conservation of energy and enstrophy (Arakawa, 1966) improves the accuracy of the nonlinear dynamics.In the vertical, the general hybrid pressure-sigma coordinate (Simmons and Burridge, 1981) is used.The forward-backward scheme is used for horizontally propagating fast waves, and an implicit scheme is used for vertically propagating sound waves.The Adams-Bashforth scheme is applied for horizontal advection of the basic dynamical variables and for the Coriolis force.In order to eliminate stability problems due to thin vertical layers, the Crank-Nicholson scheme is used to compute the vertical advection tendencies.
In high resolution Numerical Weather Prediction (NWP) applications, the computational efficiency of the model has been higher than the efficiency of most established nonhydrostatic models.The high computational efficiency of the NMM is primarily due to the design of the time-stepping procedure (Janjic, 2003).
In very high resolution two-dimensional runs, the model formulation successfully reproduces a number of classical nonhydrostatic tests.In three dimensional atmospheric runs the model dynamics demonstrates the ability to develop the observed −3 and −5/3 spectral slopes that are induced by the model physics, and not by computational noise (Janjic, 2004).In a decaying turbulence case on the cloud scales with 1 km resolution, the model dynamics develops the −5/3 spectrum consistent with the 3-D turbulence theory.
The unified version of the model was developed for a broad range of spatial and temporal scales extending from large eddy simulations (LES) to global (Janjic, 2005).The global model was formulated on the latitude-longitude grid with polar filtering of the tendencies of the basic model variables.In contrast to the WRF-NMM, which is defined on the Arakawa E grid, the dynamics of the NMMB were reformulated for the Arakawa B grid.The non-hydrostatic component of the model dynamics is introduced through an addon module that can be turned on depending on resolution.Conservative, across the pole polar boundary conditions are specified in the global limit (Janjic, 2009).In regional applications, a rotated longitude-latitude system is used.With the Equator of the rotated system running through the middle of the integration domain, more uniform grid distances are obtained in this way.
The version of the NMMB/BSC-Dust model presented in this paper can simulate weather and dust at global and regional scales.The global model can supply lateral boundary conditions for the regional version of the model run on any regional domain using the rotated latitude-longitude coordinate.The latest version of the NMMB that is being developed at NCEP can run simultaneously with multiple stationary and moving nests, including several levels of nest telescoping on the latitude-longitude grid.The 2-way interaction between the nests and their driving regional and/or global model is under development.
A variety of WRF physical parameterizations are available within the model.This variety is expected to be further extended in the future within the NOAA Environmental Modeling System (NEMS) under development at NCEP.The standard operational (and thoroughly tested in NWP and regional climate applications) physical package used includes the nonsingular Mellor-Yamada-Janjic (MYJ) level 2.5 turbulence closure for the treatment of turbulence in the planetary boundary layer (PBL) and in the free atmosphere (Janjic, 2001), the surface layer scheme based on the Monin-Obukhov similarity theory (Monin and Obukhov, 1954) with introduced viscous sublayer over land and water (Zilitinkevich, 1965;Janjic, 1994), the NCEP NOAH land surface model (Ek et al., 2003) or the LISS model (Vukovic et al., 2010), the GFDL longwave and shortwave radiation (Lacis and Hansen, 1974;Fels and Schwarzkopf, 1975), the Ferrier gridscale clouds and microphysics (Ferrier et al., 2002), and the Betts-Miller-Janjic convective adjustment scheme (Betts, 1986;Betts and Miller, 1986;Janjic, 1994Janjic, , 2000)).The vertical diffusion is handled by the surface layer scheme and by the turbulence scheme.The lateral diffusion is formulated following the Smagorinsky non-linear approach (Janjic, 1990).
As detailed in Sect.3.4 we additionally coupled the rapid radiative transfer model (RRTM) (Mlawer et al., 1997)  treat dust as a radiatively active substance interacting with shortwave and longwave radiation.We note that for this contribution, simulations were run with the GFDL radiation scheme which is used in the current operational configuration of the atmospheric model.

The dust model: BSC-Dust
The BSC-Dust is embedded into the NMMB model and solves the mass balance equation for dust taking into account the following processes: (1) dust generation and uplift by surface wind and turbulence, (2) horizontal and vertical advection, (3) horizontal diffusion and vertical transport by turbulence and convection (4) dry deposition and gravitational settling and (5) wet removal which includes in-cloud and below-cloud scavenging from convective and stratiform clouds.
Transport of dust by advection and turbulent diffusion is analogous to those of moisture transport in the NMMB (Janjic et al., 2009).The model includes 8 dust size bins with intervals taken from Tegen and Lacis (1996) and Pérez et al. (2006a).Within each transport bin, dust is assumed to have a time-invariant lognormal distribution (Zender et al., 2003a).The mass of each bin depends on model processes while the shape of the distribution is fixed to a mass median diameter of 2.524 µm (Shettle, 1986) and a geometric standard deviation of 2.0 (Schulz et al., 1998).The sub-micron particles correspond to the clay-originated aerosol (bins 1-4) and the remaining particles to the silt (bins 5-8).

Dust emission
Required input data for dust emission schemes are surface wind speed and turbulence, land use type, vegetation cover, erodibility, surface roughness, soil texture and soil moisture.This section describes the treatment of the sources and the emission scheme in the model.

Dust sources and soil size distribution
Traditionally, models have used bare ground categories of land cover maps to locate dust sources.Recently, by means of the Total Ozone Mapping Spectrometer (TOMS) aerosol index (AI) satellite retrieval, Prospero et al. (2002) showed that enclosed basins containing former lake beds or riverine sediment deposits are preferential sources that dominate the global dust emission.There are several different model representations of preferential sources (Zender et al., 2003b;Cakmur et al., 2006), based on topography (Ginoux et al., 2001), hydrology (Tegen et al., 2002), geomorphology (Zender et al., 2003b), surface reflectance retrieved from Moderate Resolution Imaging Spectroradiometer (MODIS) (Grini et al., 2005), frequency of high TOMS AI values (Westphal et al., 2009), and UV-visible albedo (Morcrette et al., 2009).In this model, we use the Ginoux et al. (2001) approach by which all topographic lows with bare ground surface are assumed to have accumulated sediments that are potential dust sources (Fig. 1).As in Ginoux et al. (2001), we are assuming a source erodibility factor S defined as representing the probability to have accumulated sediments in the grid cell i of altitude h i .Here, h max and h min are the maximum and minimum elevations in the surrounding 10 • × 10 • topography.
We also account for the seasonal changes in vegetation over semi-arid areas.Indeed, Tegen et al. (2002) have shown for example that Asian dust source strengths are particularly sensitive to the seasonality of vegetation cover.In this regard, dust emission will take place at the fraction of bare soil exposed in a grid cell A which is expressed as A = 1 − V , where V is the vegetation fraction.We use a global 0.144 • monthly climatological vegetation fraction (1985)(1986)(1987)(1988)(1989)(1990) estimated from AVHRR (Ignatov and Gutman, 1998).
To specify the soil particle size distribution required by the emission scheme we use the soil textures of the hybrid STATSGO-FAO soil map.In this database, the FAO twolayer 5-min global soil texture is remapped into a global 30-s regular lat-lon grid.Within continental United States, the soil texture is then replaced by the 30-s STATSGO data.Since soil depth significantly depends upon soil type, the dominant soil texture from 0-30 cm from multi-layer STATSGO soil was selected to match the FAO soil depths and to produce the soil texture.Table 1 displays the mass fractions of clay, silt and sand for each soil texture class estimated from the textural triangle.We use 4 soil populations in our model distinguishing among fine-medium sand and coarse sand according to criteria detailed in Tegen et al. (2002).In this sense for example, clay loams are highly unlikely to contain coarse sand while sandy clay loams could contain both coarse and medium/fine sand.It should be noted that the textural triangle is based on measurements performed by wet sedimentation techniques which break the soil aggregates leading to high amounts of loose clay particles that generally form aggregates of larger size and that may not be encountered in natural soils (Bergametti et al., 2007;Laurent et al., 2009).Following Tegen et al. (2002) we assume that clay is in the form of aggregate and it is reassigned as silt in loamy sands.

Horizontal flux and threshold friction velocity
Sandblasting and disaggregation of clay and silt particles by large particles in saltation dominate the vertical flux of dust which is strongly sensitive to the size distribution of saltating particles (Shao et al., 1993).At the sources the direct emission of small dust particles by wind is negligible.The threshold wind friction velocity of soil particles (i.e. the wind friction velocity above which soil particles begin to move in horizontal saltation flux) shows an optimum particle size Table 1.Top soil texture classes from STASGO-FAO database, mass fractions (%) of coarse sand (CS), fine-medium sand (FMS), silt (S) and clay (C) for each soil texture class, and soil gravimetric water content threshold (w ).Sand is separated into coarse sand and fine-medium sand according to Tegen et al. (2002) range for uplifting between 60 and 80 µm (Iversen and White, 1982).For smaller and larger particles the threshold friction velocity increases due to inter-particle cohesion forces and gravity, respectively.As shown by many studies (e.g., Bagnold, 1941;Gillette and Stockton, 1989;Shao et al., 1993) saltation can be considered as approximately proportional to the third power of the wind friction velocity.NMMB/BSC-Dust simulates the horizontal saltation flux H according to White (1979): where ρ a is the air density, u * is the wind friction velocity, g is the gravitational constant, u * t i is the threshold wind friction velocity and s i is the relative surface area of each soil population i.
The wind friction velocity u * is predicted by the NMMB's surface layer scheme which follows Monin-Obukhov similarity theory and includes a viscous sublayer over land.The aerodynamic roughness used in the calculation of u * is predefined for different land use categories.As discussed in Darmenova et al. (2009) the problem in using a mesoscale model predefined aerodynamic roughness instead of a aeolian roughness in the dust emission scheme is that the two roughnesses reflect processes on different scales.Laurent et al. (2006) show that for bare ground, the typical aerodynamic roughness lengths used in mesoscale models are 2 to 4 orders of magnitude higher than satellite estimates of aeolian roughness.At present, we keep the scales consistent with the mesoscale model and we do not recalculate u * based on satellite estimates of aeolian roughness length.We plan a specific study to analyze the sensitivity of the emission es- timates when using satellite-derived aeolian versus aerodynamic roughness lengths at global and regional scales.The threshold wind friction velocity is expressed as where u * tsd (D s ) is the threshold friction velocity of a smooth and dry surface with soil particles of diameter D s which is parameterized with the semi-empirical relationship of Iversen and White (1982)  where ρ p is the density of the soil particles and Re * = u * tsd D s /ν is the Reynolds number where ν is the kinematic viscosity of air.Because Eq. ( 4) is implicit, we use an alternative expression for the Reynolds number that only depends D s based on Marticorena and Bergametti (1995), i.e.Re * = 1331D 1.56 s + 0.38 where D s is in cm.We consider four mean soil particle diameters representing the soil populations, i.e., clay with 2 µm, silt with 15 µm, fine-medium sand with 160 µm and coarse sand with 710 µm.
Soil water can inhibit dust emission by increasing the threshold friction velocity of soil particles.In NMMB/BSC-Dust soil moisture effects are included following Fécan et al. (1999) through the soil moisture correction parameter f h in Eq. (3).
where w is the topsoil layer gravimetric water content.
The NMMB provides volumetric water content which is converted to gravimetric water content assuming that soil is saturated when all interparticle pores are filled with water as described in Zender et al. (2003a).The maximum amount of adsorbed water or soil gravimetric water content threshold w is an increasing function of the clay fraction in the soil.
On the basis of empirical data, Fécan et al. (1999) derived where %Clay is the percentage of clay in the soil.Table 1 displays the values of w calculated for each soil texture class in the model.Non erodible elements (pebbles, stones or vegetation) dissipate a part of the wind momentum that will not be available to promote saltation.As in the case of soil moisture, this effect can be parameterized by increasing the threshold friction velocity.Here we use the so-called drag partition correction parameter f e in Eq. ( 3) following Marticorena and Bergametti (1995) where f e expresses the efficiency with which drag is partitioned between the roughness elements characterized by an aerodynamic (or aeolian) roughness length (z 0 ) and the erodible surface characterized by a smooth roughness length (z 0s ).The latter is calculated as z 0s = D max /30 where D max is the median diameter of the coarser population of the soil size distribution (Marticorena et al., 1997).

Vertical flux
In complex sandblasting schemes, either the kinetic energy of the saltating grains exceeding a threshold of disruption of soil aggregates (e.g., Alfaro and Gomes, 2001), or the volume of soil removed by saltating grains when they impact the surface (e.g., Lu and Shao, 1999;Shao, 2001) are used to simulate the vertical dust flux and its size distribution explicitly.However, the application of such complex schemes is mainly hampered by the lack of required input data at global and regional scales (Laurent et al., 2009).In view of these large uncertainties, current models assume varying degrees of simplification in the dust emission scheme as a function of available data and in most cases several parameters are tuned to match dust observations that are mainly available far away from sources.In our model we follow the empirical relationship of Marticorena and Bergametti (1995) and Marticorena et al. (1997)  where α = 4 i=1 m i α i is the so called horizontal to vertical flux ratio reflecting the availability of dust in the soil, which is calculated as the sum of the vertical to horizontal flux ratio α i of each soil population class i weighed by their mass fraction m i in the soil.In particular we take the values given in Tegen et al. (2002) for clay (α 1 = 10 −6 cm −1 if m clay < 0.45 or α 1 = 10 −7 cm −1 if m clay ≥ 0.45), silt (α 2 = 10 −5 cm −1 ), fine-medium sand (α 3 = 10 −6 cm −1 ) and coarse sand (α 4 = 10 −7 cm −1 ).
Equation ( 8) provides a size-integrated vertical flux.We assume that the vertical dust flux is size distributed according to the 3 lognormal background source modes of D'Almeida (1987) and then distributed over each size transport bin.The vertical flux of dust F k for each size bin k is where f i = sm i F is the vertical flux for source mode i with sm 1 = 0.036, sm 2 = 0.957, sm 3 = 0.007, and M i,k is the mass fraction of source mode i carried in each transport bin k which is calculated following Zender et al. ( 2003a) where erf is the standard error function, d max,k and d min,k are the minimum and maximum diameters of transport bin k, and Dv,i = (0.832 µm, 4.82 µm, 19.38 µm) and σ g,i = (2.1,1.9,1.6) are the mass median diameter and geometric standard deviation of the source modes (D'Almeida, 1987;Zender et al., 2003a).Summarizing Sect.3.1, the total vertical mass flux of dust F k into transport bin k at each gridcell is where C is a global tuning factor.Finally, the NMMB/BSC-Dust assumes a viscous sublayer between the smooth desert surface and the lowest model layer (Nickovic et al., 2001;Janjic, 1994).In this regard, the model diagnoses the dust concentration on the top of the viscous sublayer, based on F k and the turbulent regime, which represents the lower boundary condition of the NMMB vertical diffusion scheme.More details on this scheme can be found in Nickovic et al. (2001).

Sedimentation and dry deposition
Sedimentation or gravitational settling is the most efficient removal process for large aerosols.We solve implicitly the sedimentation for the dust mixing ratio (χ) from column top (L = 1) to bottom (L = LM) as  where n is the time index, Z is layer depth, L the layer number, and v g is the gravitational settling velocity at each layer, which is calculated following the Stokes-Cunningham approximation: where v gk is the gravitational settling velocity for dust size bin k, d k is the dust diameter, ρ k is the dust density, ρ a is the air density and g is the gravitational constant.The Cunningham correction factor C c accounts for the reduced resistance of viscosity as particle size approaches the mean free path of air molecules λ.
The parameterization of the dust dry deposition at the bottom layer of the model is based on Zhang et al. (2001) which includes simplified empirical parameterizations for the deposition processes of Brownian diffusion, impaction, interception and gravitational settling detailed in Slinn (1982).Dust rebound at the surface is not taken into account due to limited knowledge of this process.
The dry deposition velocity v dk at the bottom layer is expressed as where the gravitational settling velocity v gk is calculated according to Eq. ( 13).The aerodynamic resistance R a is calculated as where z 1 is the height of the lowest model level at which the deposition is evaluated, z visc is the viscous sublayer depth, and ϕ are empirically derived stability functions, κ is the Von Karman constant and u * is the friction velocity.This term is calculated using the NMMB surface layer scheme which is based on the well established Monin-Obukhov similarity theory (Monin and Obukhov, 1954).The similarity theory requires prescription of boundary conditions at two levels in the air, i.e., z 1 and z visc .The relevant variables at z 1 are used as the upper boundary condition.In order to specify the lower boundary z visc , the model includes parameterizations of a viscous sub-layer for land (Zilitinkevich, 1965) and water (Janjic, 1994).
The surface resistance R s is calculated as where E B = Sc −γ is the collection efficiency from Brownian diffusion in the viscous sub-layer that depends on the Schmidt number Sc and an empirical constant γ varying with land use categories;  2.

Wet scavenging and convective mixing
Wet scavenging of dust by precipitation is computed separately for convective and grid-scale (stratiform) precipitation.It represents the most efficient process for the deposition of the smallest dust particles.The model includes parameterizations for in-cloud scavenging, i.e., the process by which particles rainout after nucleation scavenging; and for subcloud or below cloud scavenging, i.e., the process by which particles are washed out through collection by precipitation.
The standard cloud and precipitation schemes of the NMMB model are the grid-scale cloud microphysical scheme of Ferrier et al. (2002), and the convective adjustment scheme of Betts-Miller-Janjic (BMJ) (Betts, 1986;Janjic, 1994), both operational at NCEP.As we will detail in this section, con-vective mixing of dust closely follows the principles of the BMJ scheme.

Grid-scale deposition
The grid-scale cloud microphysical scheme in the NMMB (Ferrier et al., 2002) contains some of the functionalities of more sophisticated microphysics packages used in cloud-resolving models and high-resolution mesoscale models (e.g., Rutledge andHobbs, 1983, 1984;Reisner et al., 1998) while remaining computationally efficient.The prognostic variables in the microphysics are mixing ratios of water vapor, (nonprecipitating) cloud water, rain, and ice.The ice is a composite category composed of small, nonprecipitating ice crystals and precipitating ice particles.Throughout the rest of the NMMB, the prognostic variables are specific humidity and total condensate.This approach assumes that changes due to advection in the relative composition of cloud water, rain, and ice from the previous time step are small within each grid column (Ferrier et al., 2002).
In the NMMB/BSC-Dust, the dust wet deposition is calculated sequentially from model column top (L = 1) down to the surface (L = LM) where F w k (L) is the liquid phase wet deposition flux of dust bin k at level L, F w k (L−1) is the liquid phase wet deposition flux of dust arriving at level L from above (i.e. from level L − 1), F w k (L) is the input flux of scavenged dust at level L due to in-cloud scavenging or below or sub-cloud scavenging, f evp is the fraction of precipitation flux lost to the air column due to precipitation evaporation and α w (fixed to 0.5) is a tuning parameter to account for the fact that not all of the rain droplets will evaporate.
The input flux of scavenged dust bin k by in-cloud scavenging F w k | in is calculated as where k is a solubility parameter which expresses the fraction of dust contained in cloud and ice water that can eventually be precipitated, M k is the mass loading of dust for bin k in the gridcell, P CR is the conversion rate of cloud water to rain by autoconversion, accretion, and shedding of accreted cloud water, P IR is the conversion rate of cloud ice to precipitation through melting, f liq is the fraction of cloud water, f ice is the fraction of cloud ice, QW is the cloud water mixing ratio and QI is the cloud ice mixing ratio.
In our model scavenging is defined as the removal of dust from the grid cell, i.e. in our case activation and collection in cloud droplets or activation in ice particles do not contribute to the model in-cloud scavenging unless they rain out as precipitation.Since there are strong uncertainties related to the activation properties (solubility) of mineral dust, we use intermediate values between purely hydrophobic and purely hydrophilic assumptions found in the literature (Zakey et al., 2006).The values decrease with the increasing particle size as the small particles are more likely to form cloud condensation nuclei.Dust coating with soluble material such as sulphates would increase in-cloud scavenging.We do not consider in-cloud scavenging in the form of snow.
Below cloud scavenging for rain F w k | sub is calculated for each layer according to Slinn (1984) where c is a numerical factor = 3/2 (Loosmore and Cederwall, 2004), P l is the liquid precipitation rate, E l k is the capture efficiency of water droplets, D is the raindrop diameter and d k is the aerosol diameter.E l k incorporates the effects of directional interception, inertial impaction and Brownian diffusion.These terms are detailed in Appendix A1.
For snow, the wet deposition scheme calculates sequentially from column model top (L = 1) down to the surface (L = LM) where S k (L) is the ice phase (snow) deposition flux of dust at level L, S k (L − 1) is the ice phase (snow) deposition flux of dust arriving at level L from above, S k (L) is the input flux at level L due to snow sub-cloud scavenging (Slinn, 1984) where r is 0.6, P s is the snow rate, E s k is the capture efficiency for different types of snow, λ is the characteristic capture length, D m is the characteristic length.Details of E s k , D m and λ are given in Appendix A2.

Convective mixing and scavenging
The NMMB employs the Betts-Miller-Janjic (BMJ) convective parameterization scheme developed by Betts (1986), Betts and Miller (1986) and Janjic (1994).Deep convection is viewed here as a thermodynamically driven process that transports the heat and moisture in order to reduce and eventually remove conditional instability.In the BMJ scheme, subject to several constraints, the equilibrium deep convective clouds are represented by reference temperature (Betts, 1986) and humidity profiles (Janjic, 1994) that depend on convective regime.The shallow convection uses temperature profiles defined following Betts (1986), while the moisture profiles are specified from the requirement that the second principle of thermodynamics be satisfied (Janjic, 1994).The actual model temperature and moisture profiles are relaxed toward these reference profiles.If the deep convection cannot be sustained, the deep convection algorithm is replaced by the shallow one.Let T and q denote the changes of temperature and specific humidity within a convection time step where the subscript ref indicates the equilibrium reference profiles (Betts, 1986), the superscripts n denote the values of temperature and specific humidity at the model levels at the beginning of the time step, τ is the minimum allowed relaxation time, and F (E) (Janjic, 1994) is the cloud efficiency that depends on the convective regime which is proportional to a nondimensional combination of the entropy change over the time step, precipitation over the time step, and the mean temperature of the cloud (Janjic, 2000).
During moist convection dust particles are vertically mixed while scavenged by convective precipitation.BMJ is an adjustment scheme that describes the change in the total moisture at each layer in the column and does not describe the vertical moisture flux or entrainment within the convective plume.The BMJ scheme has been optimized over years of operational application in the NCEP Meso Eta model (Janjic, 2000) and in the WRF-NMM model for precipitation forecasts over North America.Rather than undertaking implementation of a cumulus parameterization that estimates the convective mass flux explicitly, we have implemented a new convective mixing and scavenging scheme for dust following the principles of the BMJ scheme.
In the BMJ, the precipitation produced in a convective cloud is proportional to the total change of humidity during the time step.In order to account for in-cloud scavenging we remove the dust in the convective cloud proportionally to the release of the total moisture as precipitation within the deep convective cloud during the convective time step.The input flux of dust for bin k by incloud scavenging is calculated as where Q tot is the total change of moisture in the convective cloud between calls to convection, and Q tot and M k are the total moisture and the total dust mass loading within the convective cloud at the beginning of the time step.The negative sign in the right side of Eq. ( 25) is due to the fact that Q tot is always negative (reduction of total moisture and production of precipitation in the convective column).Note that the deep convective cloud is treated here as a single layer extending over the depth of the cloud.We assume that the remaining dust is mixed vertically analogously to moisture, so that the reference vertical profile for dust preserves similarity to that of moisture.The change of dust concentration C within a convection time step t is thus where C n is the dust concentration at the beginning of the time step.C ref fullfills two conditions: (1) it follows the shape of the reference moisture profile and (2) the total mass loading of the dust reference profile equals the eventual remaining mass of dust after the in-cloud scavenging during the convective timescale (τ ).Below cloud scavenging is performed following Slinn (1984) (Eq.20) assuming a typical raindrop diameter D of 1 mm for convective precipitation.The terminal velocity of raindrops V t (D) is calculated after Willis (1984) For shallow convective clouds without precipitation we mix the dust homogeneously within the cloud.

Radiation
In order to couple aerosol and radiation processes, the RRTM radiative transfer model (Mlawer et al., 1997) including aerosol effects has been implemented into the model as an alternative option to the operational Geophysical Fluid Dynamics Laboratory radiation package (Lacis and Hansen, 1974;Fels and Schwarzkopf, 1975).RRTM is used in the GFS global operational model at NCEP.With the new radiation module dust can be treated as a radiatively active substance interacting with both short and longwave radiation.For each size bin and wavelength we specify the extinction efficiency, single-scattering albedo and asymmetry factor with a Mie-algorithm based on the work of Mishchenko et al. (2000).Each particle is assumed to be homogeneous and spherical.Although there is sufficient experimental evidence that nonsphericity of desert dust can result in significantly different scattering properties than those predicted by the Mie theory (Mishchenko et al., 2000), the effect of nonsphericity upon radiative fluxes and albedos is small (Lacis and Mishchenko, 1995).The refractive indices are taken from the Global Aerosol Data Set (GADS) (Koepke et al., 1997) modified using Sinyuk et al. (2003).GADS for dust is mainly based on Volz (1973); Levin et al. (1980) and Patterson and Gillette (1977), slightly modified in the infrared to take into account quartz absorption features in agreement with transmission measurements.Kaufman et al. (2001) suggested that the indices of refraction at solar wavelengths from Patterson and Gillette (1977) were excessively absorbing.In this sense we adopt refractive indices at solar wavelengths as an average of Patterson and Gillette (1977) and (Sinyuk et al., 2003).The latter are based upon Total Ozone Mapping Spectrometer (TOMS) retrievals and in situ Sun photometer measurements from the Aerosol Robotic Network (AERONET).

Simulations and evaluation
To evaluate the model predictions against observations we have simulated the dust distribution with the regional and global configurations for years 2006 and 2000, respectively.

Regional simulation over Northern Africa, Middle East and Europe for year 2006
We selected 2006 since it is the reference year for model evaluation and intercomparison within the Sand and Dust Storm Warning Advisory and Assessment System (SDS-WAS) for Northern Africa, Middle East and Europe (http://sds-was.aemet.es/).The SDS-WAS is a project under the umbrella of World Meteorological Organization (WMO) and intends to achieve comprehensive, coordinated and sustained observations and modeling capabilities of dust storms, in order to improve their monitoring state and increase the understanding of their formation processes.It shall enhance the ability of countries to deliver timely and quality sand and dust storm forecasts, observations, information and knowledge to users through an international partnership of research and operational communities.

Model set-up
The model domain consists of 385 × 281 grid points with a horizontal grid spacing of 0.25 • × 0.25 • and 40 vertical layers.The atmospheric model's fundamental time step was set to 40 s.In this simulation, dust advection and lateral diffusion are computed every 2 time steps, dust emission and vertical diffusion every 4 time steps, and convection and large scale precipitation (as well as dust convective mixing and wet scavenging) every 8 time steps.The dust distribution is simulated between 15 December 2005 and 31 December 2006 using 1 • × 1 • NCEP final analyses (FNL) as initial and 6-h boundary conditions.The first two weeks are discarded to remove dust spin-up effects.Meteorological initial conditions are reinitialized every 24 h with a spin-up of 12 h.In this contribution, simulations were carried with the the operational GFDL radiation scheme which does not allow feedback between dust and radiation.A tuning factor of 0.66 was used to minimize the error with respect to the observations.

Observational data
We compare the model dust AOD to daily AOD from AERONET Sun photometer stations (Fig. 3 and Table 3) and seasonal satellite aerosol distributions from OMI and MISR (Fig. 4).
From OMI we use the AI which is a measure of how much the wavelength dependence of backscattered ultraviolet (UV) radiation from an atmosphere containing aerosols differs from that of a pure molecular atmosphere with only Rayleigh scattering.In the UV, aerosol retrievals over deserts are facilitated by the intrinsically low surface reflectance at these wavelengths (Torres et al., 2002) and represent a valuable semi-quantitative product.Note that the AI depends upon height of the aerosol layer which makes comparison to model only qualitative.In the visible and near infrared (IR), deserts are highly reflective and accurate retrievals of Fig. 3. Location of the AERONET stations used in the regional evaluation.The map includes the acronyms of the stations defined in Table 3. AOD are difficult because single view multispectral satellite instruments are generally unable to separate the atmospheric and surface contributions to the measured radiances.Multiangle instruments like MISR (Meloni et al., 2004) make use of the directional properties of the surface to assist in the separation procedure.Studies have shown that MISR provides reliable AOD over desert surface and that although the MISR observation repeat time is only 3 or 4 visits per month over our region of interest, the major seasonal dust activity is captured at this temporal resolution (e.g., Martonchik et al., 2004;Kahn et al., 2010).However, one should keep in mind the current limitations of the satellite data when comparing to the modeled dust distribution.
AOD at 550 nm from AERONET stations was obtained from quality-assured data between 440 and 870 nm following the Ångström law.The typical uncertainty in the AOD measured by AERONET instruments ranges from 0.01 to 0.02 and is spectrally dependent with higher errors in the UV spectral range (Holben et al., 1998;Dubovik et al., 2000).Additionally, direct-sun AOD processing includes the Spectral Deconvolution Algorithm (SDA) described in O' Neill et al. (2003).This algorithm yields AOD of sub-micron aerosols (AODfine) and AOD of super-micron aerosols (AODcoarse) at a standard wavelength of 500 nm.The algorithm fundamentally depends on the assumption that the coarse mode Ångström exponent and its derivative are close to zero.The AODcoarse fundamentally describes the AOD of sea-salt and desert dust.Since sea-salt is related to low AOD values and mainly affects coastal stations, high AODcoarse values are related to mineral dust.These mode products are not available for all AERONET sites because at least three wavelength combinations including 490, 500 or 675 nm in addition to the standard 440 and 870 nm are needed.We compare the model at stations in the Sahel and West Africa (Fig. 5), the Middle East (Fig. 6), the Eastern Subtropical Atlantic, North Africa and Iberian Peninsula (Fig. 7) and the Mediterranean (Fig. 8).The location of the stations is shown in Fig. 3 and Table 3.The selection of the stations was based on the amount of the data during our study year and the availability of the AODcoarse product.We included 4 stations without AODcoarse (Agoufou, Banizoumbou, Cape Verde and Dakar in Fig. 5) because of their location in the Sahel and Western Africa and the predominance of dusty conditions.

Results and discussion
Figure 4 presents the geographic and seasonal distribution of the simulated dust AOD compared to the MISR AOD and the OMI AI.The data from satellites show two major dust sources, the Bodélé Basin and the Mali/Mauritania border source, specially in spring and summer.Other several important dust source areas can be identified as well, including the zone of Chotts in Algeria, several source regions in Libya, the Tokar Delta in Eastern Sudan and the south and east of the Arabian Peninsula.
Dust emission and its transport downwind varies seasonally.In fall and winter, the model reproduces the spatial extent of the average dust plume and the high AOD and AI values observed over the Bodélé, south of Niger, Nigeria, Benin, Ghana and the gulf of Guinea, when dust is mainly carried southwestward from the Bodélé and adjacent areas by the Northeasterly Harmattan winds.Strong dust events reaching maximum AOD values between 0.5 and 2.5 are observed and reproduced by the model in the stations located south of Saharan sources (Fig. 5, Niamey, Banizoumbou, Agoufou and Ilorin).These stations present large contributions of fine aerosols (with high Ångström Exponents) due to the well-known presence of fine biomass burning aerosols originating from the sub-Sahel zone in winter.In the rest of the domain, dust activity is relatively low during this period.Aerosol signatures over Libia, Sudan and the Arabian Peninsula are fairly well captured by the model.
In summer, the main dust patterns are generally well reproduced by the model from the west coast of Africa to the Arabian Peninsula.The comparison to the satellite distribution highlights the good agreement of the model in Eastern Sudan, the Tokar Delta, the Red and Arabian Seas.However, some important sources over Eastern Niger and the dust AOD at 550 nm.In the upper panels, observations at Ilorin and Niamey stations include AODtotal (green circles) and AODcoarse (red circles).Each panel includes the correlation (R), bias and root mean square error (RMSE) between the model and AODcoarse and AODtotal.In the middle and lower panels, observations at Agoufou, Banizoumbou, Capo Verde and Dakar include AODtotal (green circles) and Ångström exponent between 440 nm and 870 nm (blue circles).The top right of each panel includes R, bias and RMSE between the model and the AODtotal for Ångström exponents below 0.8 and AODtotal.0.0 0.5 1.0 1.5 2.0 2.5  Mali/Mauritania Border appear to be misrepresented mainly in summer but also throughout the year.In the end of spring and the beginning of summer the model overestimates the emission over the Bodélé leading to too high optical depths in Southern Niger, Northern Nigeria and Burkina Faso.Indeed, Fig. 5 shows that while the model tends to overestimate the AOD in the end of spring and summer in Baninzoumbou and Agoufou, it underestimates the AOD over Dakar and Cape Verde.Most probably, these regional differences are due to the overestimation of the Bodélé emissions and the underestimation of the Mali/Mauritania border emissions detected from the comparison to satellite derived data.Overall, the  daily correlations are good and range from 0.52-0.59 in Cape Verde, Dakar and Banizoumbou to 0.65-0.71 in Agoufou and Ilorin.

Dhabi
It is not straightforward to attribute the discrepancies to specific aspects of the model since the emission scheme depends on multiple surface, soil, and meteorological features and includes threshold processes and non-linear relationships.However it is clear from Fig. 1 that the Mali/Mauritania border source is mostly omitted by the topographic preferential source.A good candidate to improve this aspect of the model is the use of aeolian roughness derived from satellites as indicator for the location of preferential sources and/or in the drag partition scheme, which highlights this area as a important dust source.The problem is complex since model improvements in some regions might be accompanied by deterioration in some others.
In the eastern part of the domain, along the coastal or near coastal stations in the northeast of the United Arab Emirates (in Fig. 6 Dhabi, Dhadnah and Hamim) the model reproduces very well the daily variability of the AODcoarse with correlations of 0.75-0.77and correctly captures the seasonal distribution of the dust activity which peaks in summer.During the summer season, the southwest monsoon introduces a northwesterly flow over the Arabian Peninsula bringing extremely dry and dust-laden air from the Iraq and Southern Iran deserts (Liu et al., 2000).We can observe high total AOD and low AODcoarse, mainly in the end of fall and beginning of winter, associated with low dust activity and the influence of pollution from petroleum industry emissions in the region.In Solar Village, located in the middle of the Arabian Peninsula, near At Riyad and far from the Persian Gulf or other industrialized areas, the model underestimates some important dust storms in spring while it overestimates some episodes in summer leading to an overall correlation of 0.37.A more detailed study would be required to understand whether the model discrepancies in this site are mainly due to an inaccurate source prescription or to the inability of the model to reproduce the associated meteorology.We do not discard a measurement problem since occasionally Sun photometers miss to differentiate between high dust events and water clouds.
In the Canary Islands over the Eastern Subtropical Atlantic, at about 100 km west of the Moroccan coast, we find AERONET sites in Santa Cruz de Tenerife and La Laguna at sea level and Izana at 2391 m a.s.l.(Fig. 7).The background conditions at Izana (i.e. at periods without dust events) are associated to very low AOD values.High AOD (above 0.15) are associated to dust events which mainly are detected in summer.Santa Cruz de Tenerife and La Laguna sites are located in the city of Santa Cruz de Tenerife in the vicinity of the city harbor.The model captures most AODcoarse peaks with correlations ranging form 0.7 to 0.76.The overestimation of the dust activity in Izana is related to the calculation of the AOD from sea level in the model because of the steep topography of the island and the high altitude of the station not represented in the model.At higher latitudes, we find Blida in Algeria and El Arenosillo and Granada in the Iberian Peninsula (Fig. 6).These sites show a frequent background associated to low AOD values.High extinctions are associated with African dust events which are more frequent in spring and summer.The model reproduces the daily variability and the frequent events in the region with correlations around 0.7.
In the Central Mediterranean, the AOD at Lecce University and IMAA Potenza sites is highly affected by anthropogenic pollution (Fig. 8).Lecce University is also affected by fine particles originating from frequent summertime forest fires (Perrone et al., 2005).Saharan dust events are observed in the AODcoarse from spring to autumn and are captured by the model with correlations around 0.6-0.7.In the Central-Eastern Mediterranean, Thessaloniki is characterized by high pollution levels being strongly influenced by regional (Central and Eastern Europe) and local urban and industrial sources as well as by biomass burning (Gerasopoulos et al., 2003;Kazadzis et al., 2007;Gobbi et al., 2007).Several desert dust events are observed from spring to fall which are captured by the model resulting in a correlation of 0.71.
In the Eastern Mediterranean (Fig. 8), dust events are related to long-range transport from the Sahara, and, to a minor degree, from the Anatolian plateau, Negev desert and the Middle East (Dayan et al., 1991;Kubilay et al., 2000;Barnaba and Gobbi, 2004;Derimian et al., 2006;Basart et al., 2009).The aerosol climatology at Forth Crete site is significantly affected by the maritime environment with sea-salt aerosols constituting the background conditions.The model reproduces the daily variability of the AODcoarse at Forth Crete, Sede Boker and Nes Ziona with correlations of 0.82, 0.76 and 0.66, respectively.

Model set-up
For these experiments, the global domain was configured with an horizontal grid spacing of 1.4 • × 1 • and 24 vertical layers.Note that the global domain's expected forecast resolution is 0.47 • × 0.33 • and 64 vertical layers.The atmospheric model's fundamental time step was set to 180 s.As before, dust advection and lateral difusion is computed every 2 time steps, dust emission and vertical diffusion every 4 time steps, and convection and large scale precipitation (as well as dust convective mixing and wet scavenging) every 8 time steps.The dust distribution is simulated between 1 December 1999 and 31 December 2000 using 1 • × 1 • NCEP final analyses (FNL) as initial conditions.The first month is discarded to remove dust spin-up effects.Meteorological initial conditions are reinitialized every 24 h with a spin-up of 12 h.As in the regional configuration, simulations were carried out with the GFDL radiation scheme.A tuning factor of 2 was selected to minimize the error with respect to the observations.The difference in the tuning factor with respect to the regional simulation is mainly due to the following factors: differences in wind speed over sources due to the resolution of the model, the use of a different set of observations (type, number, distribution and temporal resolution) and the different year for each simulation.

Observational data
We evaluate monthly and annual means from a global model run for year 2000 against a variety of global observations making use of the tools developed at the LSCE within the AEROCOM project (Kinne et al., 2005;Textor et al., 2006) (http://nansen.ipsl.jussieu.fr/AEROCOM/).We use the AE-ROCOM dust benchmark data set used for global dust model evaluation and inter-comparison.A detailed description of the observations and other global dust model evaluations can be found in Huneeus et al. (2011).Here, we briefly describe the in-situ measurements of surface concentration, deposition and optical depth used in this contribution.
Monthly dust concentration measurements are used from 20 sites managed by the Rosenstiel School of Marine and Atmospheric Science from the University of Miami (Prospero et al., 1989;Prospero, 1996;Arimoto et al., 1995).Dust concentrations are derived from measured aluminium concentrations assuming an Al content of 8 % in soil dust (Prospero, 1999) or from the weights of filter samples ashed at 500 • C after extracting soluble components with water.The measurements were taken in the 1980s and 1990s with different measurement periods at each station.We also use monthly dust concentrations at Rukomechi, Zimbabwe (Maenhaut et al., 2000a;Nyanganyura et al., 2007) and Jabiru, Australia (Maenhaut et al., 2000b;Vanderzalm et al., 2003).Finally, we also use measurements from the year 2000 at Barbados station and at Miami.
The data of total deposition consists of 84 sites with yearly dust deposition fluxes not coincident with the model simulated year (Huneeus et al., 2011).We first use three compilations giving deposition fluxes over land: (1) those given in Ginoux et al. (2001) based partly upon measurements taken during the SEAREX campaign (Prospero et al., 1989) (2) those given in Mahowald et al. (2009) measuring iron and/or dust deposition and assuming a 3.5 % iron content and (3) deposition fluxes derived from ice core data (Mahowald et al., 1999).We also use a selection of deposition fluxes from sediment traps from the Dust Indicators and Records in Terrestrial and Marine Paleoenvironments (DIRTMAP) database (Tegen et al., 2002;Kohfeld and Harrison, 2001).
Deposition and surface concentration measurements are mostly not coincident with the simulated year.We follow Huneeus et al. (2011) and consider these datasets to approximate the present climatology of dust deposition.However, some of these measurements do not cover a long enough period to be climatological in a strict sense.
AOD from AERONET Sun photometers is used here as well.We concentrate on the monthly and annual means providing a comprehensive evaluation of the seasonal dust cycle and we exclude the evaluation of the frequency and intensity of dust events.The performance of the model to simulate individual dust events at a regional scale was analyzed in Sect.4.1.We use all available dusty stations with measurements for the year 2000 and a climatology constructed considering the multi-annual database 1996-2006.The selection method for dusty stations is detailed in Huneeus et al. (2011).

Results and discussion
We first analyze the model's ability to reproduce the observations from the 20 sites managed by the Rosenstiel School of Marine and Atmospheric Science from the University of Miami and at Rukomechi, Zimbabwe, and Jabirun, Australia.The model performance is assessed in terms of yearly averages and seasonal variability (Fig. 9).Following Huneeus et al. (2011) we group the stations according to their range of measured surface concentration in remote sites with measurements lower than 1 µg m −3 (orange in Fig. 9), stations under the influence of minor dust sources of the Southern Hemisphere or remote sites in the Northern Hemisphere (violet) and stations downwind of major dust sources (blue).The observations present a strong gradient between the three groups with the largest values in stations downwind of the main dust sources and the lowest ones in the remote sites.The overall correlation is high (0.87) in the upper range of AEROCOM models (Huneeus et al., 2011) mainly because the model reproduces very well the surface concentration over the sites located downwind of the major dust sources in Africa (Barbados (18), Miami (19) and Bermuda (21)) and Asia (Hedo (20) and Cheju ( 22)).The model also generally reproduces the gradient between the stations downwind of the main sources and those under the influence of minor dust sources but has difficulties in reproducing the gradient between the latter and the remote sites mainly due to an overestimation of the surface concentration in remote sites.This is reflected by the lower but still significant logarithmic correlation (0.73).The model strongly underestimates the observations in Antarctica (sites 8 and 9) and in Rukomechi, South Africa (site 17).
The seasonal variability of measured and simulated surface concentration is presented as Hovmoller-like diagrams also in Fig. 9.Each row corresponds to the monthly surface concentration of a particular station of the network.The stations are grouped as well as low, medium and high according to their surface concentration regime and we identify each group of stations with a colored bar on the left side of the figures.
In stations dominated by dust from major dust sources, the model successfully simulates the seasonal variability and magnitude of the surface concentration in the American stations of Barbados (18),Miami (19) and Bermuda (10).The simulated magnitude at these sites, while mostly overestimated in periods of maximum surface concentration, is within the variability of the observations suggesting that the model not only manages to simulate the dust transport across the Atlantic but also its latitudinal extent towards the north.Similarly, the model not only simulates the seasonal cycle of the dust surface concentration east of the Asian dust sources (Hedo (20) and Cheju ( 22)) but also the long-range transport of Asian dust into the Central North Pacific as revealed by the seasonal cycle in Midway Island (15).The model mostly underestimates in periods of maximum concentrations in the East China Sea (Hedo and Cheju) and mostly overestimate them at Midway Island.In general the simulated values are within the variability of the measurements at these three sites.Note that the model successfully reproduces the annual mean and seasonal variability for the stations downwind of the two major source regions (North Africa and Asia) at the same time, a common complex issue in global dust modeling.Finally the model does not manage to reproduce the seasonal cycle at Rukomechi ( 16) in South Africa presenting relatively constant values throughout most of the year mainly underestimating the surface concentration.
For stations labeled as low, the model overestimates the concentrations throughout most of the year except for Mawson (1) in the Antarctica were the concentration is underestimated for most of the months.Those stations (3,4,5,6,7) are located at tropical and subtropical latitudes far away from sources suggesting an overestimation of small particles, due to excessive emission, inaccurate vertical transport and/or deposition.The model presents the same dust regime in the geographically close stations of New Caledonia (2) and Norfolk Island (12) suggesting that both stations are in the southeast dust pathway from Australia documented in Mackie et al. (2008), whereas the observations attributes both stations different regimes.In this case, the difference between model and observations might be due to the climatological nature of the data and the episodic nature of the dust emissions in Australia.
For stations labeled as medium, comparison with Palmer (8) and King George (9) in the Antartica outline underestimation and/or inaccurate location of South American sources.Note that strong underestimation in these stations and Rukomechi ( 17) is a common feature in most AERO-COM dust models (Huneeus et al., 2011).The model reproduces the general features of the seasonal cycle at Cape Point (11) and Hawaii ( 14).In Mace Head (16) the model simulates the seasonal variability during the first half of the year but does not manage to do so in the second half presenting relatively constant concentrations.In Cape Grim (10) the simulated period of maximum surface concentrations is out of phase with the observed one whereas in Jabirun (13) the model presents too large monthly variability not seen in the measurements.
The model also successfully reproduces the seasonal cycle of surface concentration of the year 2000 at Barbados and Miami (Fig. 10).For most of the months the simulated concentration is within the measured variability except for the month of July where it is overestimated at both stations.The model also overestimates the surface concentration in May at Barbados anticipating the period of maximum concentration by one month.
Figure 11 shows the location of the dust deposition sites and the comparison of the model against the observations.At most of the stations the simulated deposition is within a factor 10 of the observations.The model mostly underestimates the total deposition.At individual regions, the model mostly overestimates the total deposition over the Southern Ocean, South Atlantic, Europe, the Indian Ocean and at ice core sites and mostly underestimates it in the Pacific Ocean and the North Atlantic.The normalized root mean square error (NRMS) and the correlation (0.84) lie within in the upper range performance of AEROCOM models (Huneeus et al., 2011).
We finally compare the model simulation to the climatological AOD constructed from the multi-annual database and to the data from year 2000.The AERONET stations are grouped according to their location into African stations (orange in Fig. 12), stations in the Middle East (pink) and American ones (blue).Stations not belonging to any of these groups are illustrated in black.The model is successful in reproducing the AOD gradient observed among these three groups for both data sets, climatology and year 2000.African stations present in general higher AOD than the Middle East ones and these in turn have larger values than the American Stations.The correlation is 0.88, within the upper range of AEROCOM models.The simulated AOD is mostly within a twofold of the observations.The bias is slightly negative mainly due to the significant underestimation in Kanpur (25).It overestimates the AOD in Africa and the Middle East and underestimates it in America.In terms of difference with respect to the observations, the smallest NRMS is seen in the Middle East (NRMS = 0.15) followed by African stations (NRMS = 0.20) and the largest one is seen in American stations (NRMS = 0.5).In the American stations, the influence of sea salt aerosol in the observations may be playing a role in the underestimation.
The AOD climatology shows a seasonality characterized by high AOD in Africa with maximum values from December to April in the most southern stations shifting progressively to July till September in the most Northern African stations (Fig. 13).The model reproduces this seasonal cycle with its latitudinal shift and the latitudinal gradient of AOD with larger values in the south and decreasing towards the north.The model shows a tendency to overestimate the AOD in periods of maximum AOD and underestimates it elsewhere.In the Middle East a yearly cycle with maximum AOD from May/June to September is observed.The model simulates this period of maximum AOD but with a delay of one month.Again the AOD is mostly overestimated during the period of maximum AOD.In contrast to what is seen in Africa and the Middle East, the AOD in American stations is underestimated at all stations and throughout the year.The model reproduces the seasonal cycle of larger AOD from June to September in the Caribbean associated to the trans-Atlantic dust transport, yet the model does not reproduce the seasonal cycle in Andros Island (22).In Surinam (17) the   (Huneeus et al., 2010).Bottom panels show the averaged AOD at 550 nm versus modeled one for the climatology (left) and for year 2000 (right).Root mean square error (RMS), bias, ratio of modeled and observed standard deviation (sigma) and correlation (R) are indicated in the scatter plot.Normalized mean bias and normalized root mean square error are given in parenthesis next to RMS and mean bias, respectively.Black continues line is the 1:1 line whereas the black dotted lines correspond to the 2:1 and 1:2 lines.1996-2006 (left) and for year 2000 (right).Names and locations for each selected station are given in (Huneeus et al., 2011).Bottom panels show the averaged AOD at 550 nm versus modeled one for the climatology (left) and for year 2000 (right).Root mean square error (RMS), bias, ratio of modeled and observed standard deviation (σ ) and correlation (R) are indicated in the scatter plot.Normalized mean bias and normalized root mean square error are given in parenthesis next to RMS and mean bias, respectively.Black continues line is the 1:1 line whereas the black dotted lines correspond to the 2:1 and 1:2 lines.Fig. 13.Left panels are from AERONET observations and right panels are from model output.At the top, AOD at 550 nm at dusty stations for the climatology, and at the bottom AOD at 550 nm at dusty stations for year 2000.Each row corresponds to the seasonal cycle at one of the stations.They have been grouped into African (AF, orange), Middle East (ME, violet) and American (AM, blue) stations and stations elsewhere in the world (OT, black).Each one of these groups is identified by a colored bar on the left side of the left panels.Stations are ordered from south to north within each group.The row for each station corresponds to the numbers presented in Fig. 12. Name and location of each station are given in (Huneeus et al., 2011).White color corresponds to months without measurements or months not complying with the selection criteria.66 Fig. 13.Left panels are from AERONET observations and right panels are from model output.At the top, AOD at 550 nm at dusty stations for the climatology, and at the bottom AOD at 550 nm at dusty stations for year 2000.Each row corresponds to the seasonal cycle at one of the stations.They have been grouped into African (AF, orange), Middle East (ME, violet) and American (AM, blue) stations and stations elsewhere in the world (OT, black).Each one of these groups is identified by a colored bar on the left side of the left panels.Stations are ordered from south to north within each group.The row for each station corresponds to the numbers presented in Fig. 12. Name and location of each station are given in (Huneeus et al., 2011).White color corresponds to months without measurements or months not complying with the selection criteria.model does not manage to reproduce the period from February to April with larger AOD values and limits it to the month of March.Finally, the model simulates the dust transport offshore in Western Africa throughout the year as illustrated by the AOD record in Capo Verde overestimating it from February to April.
Fewer stations are included in the analysis of data for the year 2000 since the number of available stations for this particular year is smaller.The model underestimates the average AOD in American stations.In the African stations it overestimates the AOD for the two southernmost stations and underestimates it in Dakar (3).The latter feature was already observed in the evaluation of the regional model.In terms of the seasonal variability the model mostly overestimates the AOD in Ouagadougou (1) and Banizoumbou (2) from February to May but underestimates it in Dakar (3) in September and October.In addition it also reproduces a period of high AOD in Surinam (17) but limits it to March and April.The model transports dust southwards in the winter months which is a feature most global models do not reproduce (Huneeus et al., 2011).Finally the model reproduces the year-round dust transport off Africa at Capo Verde (8) where in contrast to the climatology it slightly underestimates the AOD between February and April.The overestimation over Cape Verde between February and April when compared to the AOD climatology is explained by the strong year-to-year AOD variability in this region.The year-to-year variability can be very strong over the region in February and March as it is very sensitive to the phase of the North Atlantic Oscillation (NAO).During positive NAO winters the Azores anticyclone intensifies and the stronger easterlies over West Africa increase the dust load over the Atlantic Ocean when compared with negative NAO winters.In year 2000, the dust load was very high with NAO indexes for February and March of 4.37 and 0.54, respectively.The strong positive NAO in year 2000 explains the overestimation of the model when compared to the AOD observed climatology between 1996 and 2006.

Summary and conclusions
We present the NMMB/BSC-Dust, a new online multiscale dust model prepared for regional and global simulation domains.We compare a regional simulation covering Northern Africa, Middle East and Europe to daily AOD observations from the AERONET Sun photometer network and we show that the model reproduces significantly well the daily variability and seasonal spatial distribution of the dust in the regional domain.Correlations of the simulated against observed AOD are high in general and vary depending on the region.The model is particularly good in the Eastern Subtropical Atlantic (0.7-0.76), North Africa and Iberian Peninsula (0.7-0.71), the Central (0.59-0.72) and Eastern Mediterranean (0.66-0.82) and the Sahel (0.65-0.71).By comparison to satellite derived data, the model shows its ability to reproduce the dust spatial distribution in Eastern Sudan, the Tokar Delta and the Red and Arabian seas.However it misrepresents dust sources over Eastern Niger and the Mali/Mauritania border decreasing the AOD correlation over the west coast of Africa in Dakar and Cape Verde (0.52-0.59).A slight underestimation in Cape Verde for year 2000 is also observed in the global study.The topographic preferential source map used in the model does not reflect the Mali/Mauritania Border source and in our model leads to underestimation of the emission.Other choices of preferential sources or the use of satellite derived roughness lengths in the drag partition scheme will be the object of a forthcoming study where we expect to improve and further understand this aspect of the model.
The global model is compared to variables with a direct link to the estimation of the direct radiative effect, the dust impact on the biogeochemical cycle and air quality, i.e., AOD, deposition and surface concentration.The AERO-COM dust benchmark data set has already been presented in Huneeus et al. (2011) where it is used for a multi-model intercomparison of a total of 15 global aerosol models.Here we use four different dust deposition compilations of total annual flux, annual and monthly averaged dust surface concentration across the world and AOD at dusty sites from the AERONET network.
The annual correlations of the model with observations are high (0.87 for surface concentration, 0.84 for deposition and 0.88-0.89for AOD) and lie in the upper range of the AE-ROCOM model evaluation performance scores.The model reproduces the annual mean and the seasonal variability of the surface concentration in stations downwind of the major dust sources (North Africa and Asia) while it strongly underestimates the dust concentration at the Antarctica stations affected by South American sources and at the Rukomechi station affected by the Kalahari desert in South Africa.Note that the model includes a global tuning factor and that no regional tuning factors have been applied to correct these deficiencies.To improve the representation of the Kalahari desert, the source map will be revisited.In the case of South Amer-ica, we expect to improve the source estimation at higher model spatial resolution.
The tendency of the model to overestimate the very low background concentrations far away from sources points towards an overestimation of the smallest dust particles (clay) due to either inaccuracies in the size distribution of the emissions, vertical transport and/or removal.Cakmur et al. (2006) apply separate tuning factors for the emission of clay and silt globally to minimize the error of the model against a variety of global climatological observations.In the future, we plan to optimize our global emission estimates with this method.
The model simulates the annual mean dust deposition within a factor 10 with respect to observations which was also found in the global dust model intercomparison of Huneeus et al. (2011).The model mostly overestimates the total deposition over the Southern Ocean, South Atlantic, Europe, the Indian Ocean and at ice core sites and mostly underestimates in the Pacific Ocean and the North Atlantic.Dust removal by stratiform and convective rain is very sensitive to the prescribed dust solubility and other uncertain parameters in the model.Extensive evaluation and sensitivity tests of these processes in the model need to be performed.In this sense, the future inclusion of chemistry in the model is expected to improve the representation of dust solubility.As discussed in Huneeus et al. (2011), data issues cannot be discarded since deposition measurements considered in this study are not from year 2000 and do not represent the deposition climatology in a strict sense.
The model reproduces well the AOD gradient among the different dusty regions and the seasonal cycle of Northern African and Middle Eastern dust.We note that the model reproduces the southward displacement of the Saharan dust cloud during winter in contrast to most AEROCOM models.
The developments showed in this study represent a part of a larger effort towards the development of a unified chemical weather forecast system including other aerosol components and gas-phase chemistry (Jorba et al., 2011).
In November 2011, the NMMB/BSC-Dust has started to provide pre-operational dust forecasts at the Barcelona Supercomputing Center.

Appendix A Collection efficiencies A1 Rain
The capture efficiency of water droplets E l k (Slinn, 1984) incorporates the effects of directional interception (E int ), inertial impaction (E imp ) and Brownian diffusion (E BD ) where ρ p is the density of particle, C c is the Cunningham slip factor, k b is the Boltzmann constant and T is the absolute temperature.
In the NMMB microphysical scheme, the mean rain drop diameters vary from 0.05 mm to 0.45 mm.The terminal velocity or fall speed is calculated as a function of mean drop diameter and provided through a lookup table assuming exponential size distributions for raindrops.In order to increase the computational efficiency of the NMMB/BSC-Dust we have also created lookup tables for the capture efficiency of water droplets depending on fall speed of raindrops.

A2 Snow
The capture efficiency of snow E s k (Slinn, 1984) is parameterized according to where α = 2/3 and the characteristic capture length λ = 100 µm are set to typical values of rimed crystals (Gong et al., 1997).In order to increase the computational efficiency we have also created lookup tables for the capture efficiency of snow depending on fall speed of snow and rime factor.

Fig. 1 .
Fig. 1.Preferential sources in the model.The resolution of the map is 1.4 • × 1 • which is used for the global simulation in Sect. 4.

Fig. 2 .Fig. 2 .
Fig.2.Clay, silt, fine-medium sand and coarse sand mass fractions based on the soil textures of the hybrid STATSGO-FAO soil map.Table1provides the mass fractions of clay, silt, fine-medium sand and coarse sand for each soil texture class estimated from the textural triangle.

Fig. 5 .
Fig. 5. Daily AOD average of model vs. AERONET stations in the Sahel and Western Africa for 2006.Model averages (black lines) representdust AOD at 550 nm.In the upper panels, observations at Ilorin and Niamey stations include AODtotal (green circles) and AODcoarse (red circles).Each panel includes the correlation (R), bias and root mean square error (RMSE) between the model and AODcoarse and AODtotal.In the middle and lower panels, observations at Agoufou, Banizoumbou, Capo Verde and Dakar include AODtotal (green circles) and Ångström exponent between 440 nm and 870 nm (blue circles).The top right of each panel includes R, bias and RMSE between the model and the AODtotal for Ångström exponents below 0.8 and AODtotal.

Fig. 6 .
Fig. 6.Daily AOD average of model vs. AERONET stations in the Middle East (Dhabi, Dhadnah, Hamim and Solar Village) for 2006.Model averages (black lines) represent the dust AOD at 550 nm.Observations include AODtotal (green circles) and AODcoarse (red circles).Each panel includes the correlation (R), bias and root mean square error (RMSE) between the model and AODcoarse and AODtotal.

Fig. 7 .
Fig. 7. Daily AOD average of model vs. AERONET stations in the Eastern Subtropical Atlantic (La Laguna, Santa Cruz de Tenerife and Izana), North Africa (Blida) and Iberian Peninsula (El Arenosillo and Granada) for 2006.Model averages (black lines) represent the dust AOD at 550nm.Observations include AODtotal (green circles) and AODcoarse (red circles).Each panel includes the correlation (R), bias and root mean square error (RMSE) between the model and AODcoarse and AODtotal.

Fig. 8 .
Fig. 8. Daily AOD average of model vs. AERONET stations in the Central (Lecce University and IMMA Potenza) and Eastern Mediterranean (Forth Crete, Thessaloniki, Sede Boker and Nes Ziona) for 2006.Model averages (black lines) represent the dust AOD at 550 nm.Observations include AODtotal (green circles) and AODcoarse (red circles).Each panel includes the correlation (R), bias and root mean square error (RMSE) between the model and AODcoarse and AODtotal.

Fig. 9 . 58 Fig. 9 .
Fig. 9. Network of stations measuring surface concentration in the upper left panel.Names and locations for each selected station are given in Huneeus et al. (2010).In the upper right panel, yearly averaged measured surface concentration versus modeled one at each station in µg/m 3 .It includes the root mean square error (RMS), bias, ratio of modeled and observed standard deviation (sigma) and correlation (R).Normalized mean bias and normalized root mean square error are given in parenthesis next to RMS and mean bias, respectively.The logarithmic correlation is given in parenthesis next to R. Black continues line is the 1:1 line whereas the black dotted lines correspond to the 10:1 and 1:10 lines.Bottom panels show monthly averages of measured (left) and simulated (right) surface concentration.Each row corresponds to the seasonal cycle at one of the stations.White color corresponds to month without measurements.58

Fig. 10 .
Fig. 10.Measured (blue lines) and simulated (red lines) surface concentration in Barbados (upper panel) and Miami (bottom panel) for year 2000.Standard deviation bars in blue indicate the interannual variability of the records in those locations.Units are µg m −3 .

Fig. 11 .Fig. 11 .
Fig. 11.Measured yearly deposition fluxes versus modeled ones in the bottom panel; units are g/m 2 /yr.Location for each data point in the scatter plot is given in the left panel.Number and letters are coloured regionally for West/East Pacific (red/brown), North/Tropical/South Atlantic (orange/black/light-blue), Middle East/Asia/Europe (violet/purple/light green), Indian/Southern Ocean (dark green/dark blue) and pink ice core data in Greenland, South America and Antartica.Data from Ginoux et al. (2001)/Mahowald et al. (2009)/DIRTMAP/Mahowald et al. (1999) are indicated by letters/non-italic numbers/ italic numbers/lower-case letters.Root mean square error (RMS), bias, ratio of modeled and observed standard deviation (sigma) and correlation (R) are indicated for each model in the lower right part of the scatterplot.Normalized mean bias and normalized root mean square error are given in parenthesis next to RMS and mean bias, respectively.The correlation with respect to the logarithm of the model and of the observation is also given in parenthesis next to R. Black continues line is the 1:1 line whereas the black dotted lines correspond to the 10:1 and 1:10 lines.

Fig. 12 .
Fig. 12. Upper panels show the location of selected AERONET dusty sites based on the climatology build from the multi-annual database 1996-2006 (left) and for year 2000 (right).Names and locations for each selected station are given in(Huneeus et al., 2010).Bottom panels show the averaged AOD at 550 nm versus modeled one for the climatology (left) and for year 2000 (right).Root mean square error (RMS), bias, ratio of modeled and observed standard deviation (sigma) and correlation (R) are indicated in the scatter plot.Normalized mean bias and normalized root mean square error are given in parenthesis next to RMS and mean bias, respectively.Black continues line is the 1:1 line whereas the black dotted lines correspond to the 2:1 and 1:2 lines.

Fig. 12 .
Fig. 12. Upper panels show the location of selected AERONET dusty sites based on the climatology build from the multi-annual database1996-2006 (left)  and for year 2000 (right).Names and locations for each selected station are given in(Huneeus et al., 2011).Bottom panels show the averaged AOD at 550 nm versus modeled one for the climatology (left) and for year 2000 (right).Root mean square error (RMS), bias, ratio of modeled and observed standard deviation (σ ) and correlation (R) are indicated in the scatter plot.Normalized mean bias and normalized root mean square error are given in parenthesis next to RMS and mean bias, respectively.Black continues line is the 1:1 line whereas the black dotted lines correspond to the 2:1 and 1:2 lines. .

Table 2 .
Zhang et al. (2001)us of collectors A and parameters α d and γ in the dry deposition scheme dependent on NMMB model land use type.Values were reassigned to 27 WRF-USGS land use types from the 15 land use types used inZhang et al. (2001).
Zhang et al. (2001) efficiency for impaction which depends on the Stokes number St, which takes the form St = V g u * /gA for vegetated surfaces and St = V g u * 2 /µ otherwise; and E IN = 0.5(d p /A) 2 is the collection efficiency by interception.αd is an empirical parameter and A is the characteristic radius of collectors, both dependent on land use.The parameters A, α d and γ depending on 15 land use categories provided inZhang et al. (2001)were reassigned to the 27 NMMB model's WRF-USGS land use categories as shown in Table

Table 3 .
Names, acronyms and coordinates of the AERONET stations used in the regional evaluation.
• E) Lat( • N) Alt.(m) Pérez et al.:Dust modeling from meso to global scales -Part 1 Sc 1/3 + 0.16Re 1/2 Sc 1/2 (A4) where d k is the particle diameter, D is the raindrop diameter, µ a is the viscosity of air, µ w is the viscosity of water, ρ a is the density of air.The Reynolds number of raindrops Re, the Strokes parameter of the collected particles St, the Critical Strokes number St * and the Schmidt number for collected particles Sc are expressed as (D) is the terminal velocity of raindrops, v t (d k ) is the terminal velocity of particles.The characteristic relaxation time of the particle τ and its Brownian diffusivity are www.atmos-chem-phys.net/11/13001/2011/Atmos.Chem.Phys., 11, 13001-13027, 2011 13022 C.