Parameterization of Dust Emissions in the Global Atmospheric Chemistry-climate Model Emac: Impact of Nudging and Soil Properties

Airborne desert dust influences radiative transfer , atmospheric chemistry and dynamics, as well as nutrient transport and deposition. It directly and indirectly affects climate on regional and global scales. Two versions of a pa-rameterization scheme to compute desert dust emissions are incorporated into the atmospheric chemistry general circulation model EMAC (ECHAM5/MESSy2.41 Atmospheric Chemistry). One uses a globally uniform soil particle size distribution, whereas the other explicitly accounts for different soil textures worldwide. We have tested these two versions and investigated the sensitivity to input parameters, using remote sensing data from the Aerosol Robotic Network (AERONET) and dust concentrations and deposition measurements from the AeroCom dust benchmark database (and others). The two versions are shown to produce similar atmospheric dust loads in the N-African region, while they deviate in the Asian, Middle Eastern and S-American regions. The dust outflow from Africa over the Atlantic Ocean is accurately simulated by both schemes, in magnitude, location and seasonality. Approximately 70 % of the modelled annual de-position data and 70–75 % of the modelled monthly aerosol optical depth (AOD) in the Atlantic Ocean stations lay in the range 0.5 to 2 times the observations for all simulations. The two versions have similar performance, even though the total annual source differs by ∼ 50 %, which underscores the importance of transport and deposition processes (being the same for both versions). Even though the explicit soil particle size distribution is considered more realistic, the simpler scheme appears to perform better in several locations. This paper discusses the differences between the two versions of the dust emission scheme, focusing on their limitations and strengths in describing the global dust cycle and suggests possible future improvements.


Introduction
Desert dust is an important atmospheric constituent, given its potential to affect air quality, nutrient deposition and climate (Sokilik and Toon, 1996;Tegen and Lacis, 1996;Ramanathan et al., 2001;Mahowald et al., 2005;Forster et al., 2007).Global models, being independent of lateral boundary constraints and less dependent on initial conditions than limited area models, are useful tools to study the large-scale transport dynamics, the physico-chemical behaviour and deposition distributions of dust aerosols.Furthermore, anthropogenic influences, including interactions between pollutant gases and aerosols with dust particles can be analysed, and their role in atmospheric chemistry and climate change simulated (Kallos et al., 2007;Astitha et al., 2010).The uncertainties associated with these processes and research questions are indicated by the large variety of global models using different emission parameterization schemes, input parameters and representations of aerosol removal processes (Ginoux et al., 2001;Tegen et al., 2002;Zender et al., 2003;Stier et al., 2005;Pringle et al., 2010).In many cases models are tuned Published by Copernicus Publications on behalf of the European Geosciences Union.
(either the dust fluxes or the threshold friction velocity) towards available observations, which are often performed at large distances from the main dust source areas (Ginoux et al., 2001;Zender et al., 2003;Prigent et al., 2005;Heinold et al., 2008;Li et al., 2008;Ridley et al., 2012;Perez et al., 2011).An important effort to compare the dust distributions from different global models, the Aerosol Comparisons between Observations and Models (AeroCom) project, has revealed that models apply a wide range of global dust sources and sinks (514 to 5999 Tg yr −1 ) and simulate very different atmospheric burdens (8.2 to 54 Tg) (Textor et al., 2006(Textor et al., , 2007;;Prospero et al., 2010;Huneeus et al., 2011), indicative of the uncertainties involved.
Mineral dust particles enter the lower atmosphere primarily through a mechanism called saltation bombardment, which is strongly dependent on the meteorological conditions near the surface, as well as the soil texture and particle size classification (Shao et al., 1993;Alfaro and Gomes, 2001;Grini et al., 2002).The emissions of dust particles have important regional and global consequences, whereas the phenomenon itself is of episodic nature driven by processes on small spatial and temporal scales.Therefore, global models need to be based on a number of assumptions to simplify and generalize the dust emission schemes.Several methods have been proposed to estimate the dust flux into the atmosphere, some more detailed than others (Marticorena and Bergametti, 1995;Marticorena et al., 1997;Shao, 2004;Shao et al., 1993Shao et al., , 2011;;Alfaro and Gomes, 2001;Nickovic et al., 2001;Tegen et al., 2002;Zender et al., 2003;Balkanski et al., 2004).In most cases a process-based dust emission scheme is pursued that explicitly takes into account the soil surface characteristics.One of the limitations in a global setup involves the availability of input parameters for which measurement data are lacking.While this applies to both global and regional models, detailed datasets have been collected of soil characteristics for specific desert areas (Callot et al., 2000;Laurent et al., 2005Laurent et al., , 2006Laurent et al., , 2008)), and regional models are typically more sophisticated in representing dust emissions as their resolution is closer to the small-scale dynamical processes involved.Further, the difficulty of directly measuring the dust emission fluxes in the source areas hinders the evaluation of model results.
The goal of our study is to develop and test two versions of a dust emission parameterization scheme, implemented into the atmospheric chemistry -general circulation model EMAC (Jöckel et al., 2005(Jöckel et al., , 2006(Jöckel et al., , 2010;;Roeckner et al., 2006), and investigate the dependence of the dust distribution on the soil properties and the emitted particle size distribution.The two versions of the dust emission scheme are primarily different in the explicit representation of soil particle size distributions and also in the emitted size distribution at the source.In the first version, the particle size distribution at the source is globally uniform, whereas in the second, it is explicitly accounted for based on the Zobler geographical soil texture classification and four soil populations.
We address (i) the physical processes that lead to the injection of dust particles into the atmosphere and (ii) the role of the input parameters in representing the spatial heterogeneity of dust emissions.One advantage of using EMAC is the direct coupling to meteorological calculations at each time step (∼ 10 min), which is expected to realistically represent the grid-scale temporal variability, e.g., compared to off-line calculations with a chemistry-transport model based on 3 or 6-hourly meteorological analyses.We thus combine the meteorological variables with specific input fields of soil properties without using a-priori (or preferential) source distributions of dust or pre-calculated tables of input variables.The first objective is to implement a dust emission scheme that makes use of the air temperature, humidity, density and friction velocity, thus directly linking into the meteorological calculations, and simulate the threshold friction velocity and dust fluxes online.The second objective is to investigate the role of the soil particle size distribution and texture.The spatial heterogeneity of the dust emissions is related to the uniqueness of the regions the particles originate from.This is reflected in the soil texture and size distribution, the local meteorological conditions, the topographical characteristics and the general synoptic conditions.
The paper is organized in five sections.The next section describes the EMAC model and the configuration used for the simulations.In Sect.3, the dust production scheme is presented and analysed considering the two implementations.The observational data, collected to compare with the model results, is discussed in Sect.4, which includes in-situ measurements from the AeroCom dust benchmark dataset (Huneeus et al., 2011) and remote sensing data.In Sect. 5 the effects of nudging the model to meteorological analyses are discussed, and an extensive model evaluation with subsections on concentrations, deposition fluxes and aerosol optical depth is presented.A discussion about the differences between the two versions of the parameterization is also included.

EMAC model configuration
The EMAC model combines the ECHAM5 general circulation model (Roeckner et al., 2006) and the Modular Earth Submodel System (Jöckel et al., 2005(Jöckel et al., , 2006(Jöckel et al., , 2010)).The MESSy system (version 2.41) is modular and all sub-models follow strict coding standards to allow easy implementation within other models.It also provides the option of running the model with multiple representations of processes to systematically test and improve the results.The model was run with a spectral resolution of T106 (∼ 1.1 • × 1.1 • ) and 31 vertical levels up to 10 hPa for the year 2000.The reason for selecting this year is the availability of dust measurement data (for Miami, Barbados, Haifa), important for the evaluation procedure and the recently published work on global modelling of dust for the same year (Huneeus et al., 2011;Perez et al., 2011).Gläser et al. (2012) used four horizontal resolutions, with T106 being the highest, and concluded that important differences occur in the dust emissions, which are not well described in the coarser resolution setups.The model output is recorded every 5 h and the output fields represent averages of the previous 5 h.The simulations performed include a simplified sulphur chemistry scheme allowing the production of sulphuric acid and particulate sulphate, which play an important role in transforming the dust particles from hydrophobic into hydrophyllic, thus affecting their ability to interact with clouds and be removed by precipitation.More details on the sulphur chemistry mechanism and the set of chemical reactions can be found in Gläser et al. (2012).The aerosol microphysics sub-model is M7 (Vignati et al., 2004;Stier et al., 2005), which describes the aerosol size distribution according to 7 lognormal modal size fractions; 4 soluble and 3 insoluble, encompassing the nucleation, Aitken, accumulation and coarse modes.
The emissions of gases and aerosols are treated using online and offline routines (Pozzer et al., 2009) and the processes taken into account in the simulations include advection, convection, deposition (wet and dry), simple sulphur chemistry and radiation.The set of sub-models activated in the EMAC version are described in Table 1, including references to the original work.The following emission submodels are used: (a) the fixed (offline) emissions include sulphur dioxide from anthropogenic, volcanic sources and biomass burning; nitrogen oxide from anthropogenic, aircraft, biogenic sources and biomass burning; black carbon and organic carbon from wildfires; dimethyl sulphide from terrestrial sources; formaldehyde, formic acid and methanol from anthropogenic sources and biomass burning.(b) The online emissions include dimethyl sulphide from water bodies, nitrogen oxide from soil biogenic sources, organic and black carbon (bio-fuel, fossil fuel, and secondary species), dust and sea salt.The present work extends this scheme with an online representation of mineral dust sources.The aerosols are not radiatively coupled to the model in this work; the full aerosol chemistry and thermodynamics scheme is not implemented in the simulations for computational efficiency as the focus is to investigate and evaluate the desert dust production, without the complexity of climate and chemistry feedbacks, which will be the subject of planned future work.
The modelled AOD is calculated at 550 nm using concentrations of dust and sea salt particles, biomass burning products (black carbon and organic carbon) and sulphate aerosols.More specifically, the aerosol optical properties are calculated with the EMAC submodel AEROPT (Table 1).It is based on the scheme by Lauer et al. (2007) and makes use of predefined lognormal modes (i.e. the mode width σ and the mode mean radius have to be taken into account), for which lookup tables with the extinction coefficient, the single scattering albedo and the asymmetry factor for the shortwave and extinction coefficient for the longwave part of the spectrum are created (Pozzer et al., 2012).The consid-ered species are organic carbon, black carbon, dust, sea salt, water-soluble compounds (WASO, i.e. all other water soluble inorganic ions, e.g.: NH + 4 , SO 2− 4 , HSO − 4 , NO − 3 ) and aerosol water (H 2 O).The refractive indices for dust are taken from Hess et al. (1998).More details on the method of the AOD calculation can be found in Pozzer et al. (2012).The anthropogenic fraction of the aerosol species is less relevant, since the subject of this work is to investigate the dust production and transport, and we focus on areas primarily influenced by desert dust outbreaks.
The model performance, simulating the global dust distribution, was also tested using the nudging option, in contrast to the free-running mode.As pointed out by Timmreck and Schulz (2004), significant differences can occur between nudged and free-running simulations, in particular with respect to the mean geographical distribution and seasonal variability of mineral dust.For the nudged simulations the European Centre for Medium-Range Weather Forecasts (ECMWF) 40-yr re-analysis (ERA-40) data for the year 2000 has been used.The prognostic variables nudged towards the observations (i.e., the re-analyses) are temperature, vorticity, divergence and surface pressure, and the nudging weights are chosen such that the boundary layer and the upper troposphere/lower stratosphere are not directly influenced (Lelieveld et al., 2007).A discussion about the effects of nudging on the simulated dust concentrations and deposition is included in Sect. 5.

Dust emission parameterization
The previous dust emission scheme in the EMAC model was based on the work of Balkanski et al. (2004) as discussed in Stier et al. (2005).It includes three pre-calculated tables of clay content, emission source strength factors (kg s 2 m −2 ) and threshold wind friction velocities (m s −1 ), representing the entire globe.The modelled temperature and precipitation fields are used to adjust the soil wetness and thus limit the dust production in wetted soil areas.The vertical dust flux (only coarse dust particles were emitted) was calculated with the use of the diagnosed wind speed at 10 m, the threshold velocity and the emission source strength factor from the precalculated tables.The results of the Balkanski dust emission scheme are also discussed in Gläser et al. (2012), comparing it with the Tegen et al. (2002) scheme recently implemented in the EMAC model.
The methodology followed in this work is based on previous dust emission schemes for regional (Perez et al., 2006;Spyrou et al., 2010;Laurent et al., 2008Laurent et al., , 2010;;Marticorena et al. 1997) and global modelling systems (Zender et al., 2003;Tegen et al., 2002).Two versions of a dust emission parameterization have been included into EMAC, presented in detail in the following sub-sections.

Dust sources and input parameters
Dust particles are injected into the atmosphere through saltation and sandblasting (Marticorena and Bergametti, 1995;Marticorena et al., 1997;Grini et al., 2002;Zender et al., 2003).The saltation process is initiated when the drag near the surface exceeds the gravitational inertia of the sand-size particles (diameter > 60-80 µm) moving them downwind horizontally.With this movement the large particles disaggregate and release smaller size silt and clay particles (sandblasting) (Grini et al., 2002;Alfaro and Gomes, 2001).The height of the saltation layer is of the order of 1m, which underscores that the dust emissions take place on small spatial scales.The direct emission of small and coarse dust particles, referred to as aerodynamic entrainment (Shao, 2004;Shao et al., 2011), is negligible because of cohesive and gravitational forces, respectively, which bind the particles to the soil; this mechanism is not included here considering the negligible contribution compared to saltation bombardment.A third mechanism of dust entrainment into the atmosphere is disintegration (self-abrasion) of large aggregates or fragments (Shao et al., 2011), which is considered difficult to model on a global scale due to the lack of input data that characterize the aggregates in the soil.The dust particles are considered as mobilized in the atmosphere when the wind friction velocity, a proxy of the surface drag properties, exceeds a threshold value.This threshold value depends on the soil size distribution and soil texture classification.Details on the calculation of the threshold friction velocity are given in Sect.3.2.
For the new implementation two formulations of the online dust production are tested, from here on referred to as DU1 and DU2.DU1 utilizes a homogeneous global soil size distribution of dust particles and DU2 uses an explicit geographical representation instead.Another difference between  (Olson, 1992) as described in Table 2.This database provides global information on the location of deserts (sand or clay), semi-deserts (steppe, shrub, sparse grass) and flat desert playas (Fig. 1, left plot).The set of input parameters also includes the clay fraction of the soils (Scholes and Brown de Colstoun, 2011) as shown in Fig. 1 (right plot), the rooting depth (Schenk and Jackson, 2009) and the monthly vegetation area index (sum of leaf and stem area index) as discussed in Zender et al. (2003).The schemes use the online meteorological fields from the EMAC model: temperature, pressure, relative humidity, soil moisture and the surface friction velocity.Each relevant parameter will be discussed in the following subsections.

Threshold friction velocity
A central part of the dust production scheme is the calculation of the threshold friction velocity (u * thr ), above which the emission of dust particles into the air is considered possible.It is based on an empirical relationship derived by Marticorena and Bergametti (1995) and analyzed in Marticorena et al. (1997) over smooth surfaces based on the proposed formulations of Iversen and White (1982).The relationship utilizes the friction Reynolds number B, which depends on the soil particle size D p , in Eqs.(1-3).This relationship indicates that the minimum threshold friction velocity occurs for soil particle diameters around 60-70 µm (Fig. 2a).
For 0.03 < B < 10 : (1) where K = where u * ts is the threshold friction velocity over smooth surfaces, D p is the diameter of the soil particle, v is the kinematic viscosity of air, ρ p is the particle density, ρ α is the air density and g the gravitational acceleration.The above Eqs.(1-3) are implemented in the first version of the dust emission scheme (DU1), with the iterative Eq. (2) being used in the second time step whereas the analytical solution ( 3) is used at the start of the calculation.The difference with the DU2 scheme is that in DU1 the D p is constant and equal to an optimal diameter for saltation (60 µm), when in DU2 the soil size distribution is included as an input field, and the threshold wind friction velocity is calculated for particle diameters in the range of 0.1 to 1000 µm.Two corrections are imposed in the calculated u * ts which are related to the drag partition scheme near the surface (Marticorena et al., 1997) and the soil moisture (Fecan et al., 1999).Since the initial computation of the u * ts is an empirical relation for smooth surfaces, a correction is imposed depending on the surface roughness length (z o ) and the local roughness length of the uncovered surface (z os ).The empirical relation is shown in Eq. ( 4) and it is valid for small values of aeolian roughness (z o < 1 cm) (Darmenova et al., 2009).An example of how the correction factor (1/f drag is the multiplication factor for u * ts ) changes with the surface roughness length z o and local roughness length z os is given in Fig. 2b.As shown in the graph, the correction factor is higher for higher z o and lower z os , which leads to higher values of u * thr (Eq.5).The reason for this is that a higher threshold friction velocity must be assigned for a surface with more obstacles.When the soil includes an increased number of non-erodible elements (solid obstacles, i.e. rocks, pebbles, vegetation), which translates into higher values of z o , the threshold friction velocity increases causing a decrease in the emission of dust particles.When z os increases (i.e., smoother surfaces) the correction factor decreases (1/f drag ), giving smaller values of u * thr .
For both dust emission schemes DU1 and DU2, globally uniform values have been used for z o and z os (0.01 and 0.00333 cm, respectively, as in Zender et al., 2003), given the lack of data for these parameters in the global scale.In future implementations, we plan to substitute these fixed values with more appropriate ones.
The final correction in the calculation of the threshold friction velocity is the soil moisture adjustment as proposed by Fecan et al. (1999) and applied by Marticorena et al. (1997) and also Zender et al. (2003), Laurent et al. (2005Laurent et al. ( , 2006Laurent et al. ( , 2008)), among others.The principle behind this correction is that the threshold friction velocity must increase in wetted soils and this is accomplished by relating the residual soil moisture to the clay content of the soil.The residual soil moisture w is calculated with the empirical Eq. ( 6), and its physical meaning is given as the ratio of the mass of water to the mass of dry soil: Besides the residual soil moisture, the soil moisture correction scheme requires the gravimetric soil moisture w, which is calculated from the modelled soil moisture w s (m) by dividing with the rooting depth (Hillel, 1980).The final correction of the threshold friction velocity is calculated based on Eq. ( 7), depending on whether the soil is dry or wet.
For w < w (dry soil) : (7) For w > w (wet soil) : Both emission schemes DU1 and DU2 use the above soil moisture correction in the formulation of Eqs. ( 6) and ( 7), since it is not dependent on the soil size distribution.The sequence described above provides the threshold friction velocity u * t that enters the calculation of horizontal and vertical fluxes, as discussed in the following subsection.

Horizontal and vertical flux calculations
The final step is the calculation of horizontal (H ) and vertical (V ) fluxes of the dust particles entering the atmosphere.Following Kawamura (1964) and White (1979) as implemented in Marticorena et al. (1997) the horizontal flux is calculated with Eqs. ( 8) or ( 9), depending on whether or not the soil size distribution is accounted for: where c = 1 (Darmenova et al., 2009), ρ air is the air density, u * is the friction velocity, u * t is the threshold friction velocity, S rel is the relative surface area covered from particles with diameter D p (assuming particles of spherical shape).
In DU1 the horizontal flux is calculated using Eq. ( 8) (Zender et al., 2003, Spyrou et al., 2010), since there is no size distribution assigned to the soil particles.In DU2, Eq. ( 9) is implemented to estimate the horizontal flux per soil particle of size D p and the total horizontal flux per soil source i, H tot,i .
The vertical flux V is defined as the mass of dust particles emitted from unit area per unit time (kg m −2 s −1 ).V is considered proportional to the horizontal flux H and is calculated by the relationship V = a • H established from combined measurements of H and V by Gillette (1979), accounting for all dust particles with diameters less than 20 µm.The parameter a is the sandblasting efficiency that depends on the clay content of the soil.Marticorena and Bergametti (1995) established an empirical relationship between the ratio V /H and the soil clay content, which is adopted in this work: a = 10 0.134 (% clay)−6 (10) Equation ( 10) is valid for clay percentages up to 20 %.For higher clay contents the sandblasting efficiency is set constant to the values proposed by Tegen et al. (2002).Specifically: The final equation that provides the dust vertical flux V for both DU1 and DU2 is described in the following sub-section, in which we discuss the particle size distribution, which is stated to be of central importance by Kok (2011).

Soil and transport size distributions
The main difference between the two implementations involves the soil size distribution and the emitted size distribution.In the DU1 version of the scheme we adopt an optimum size for saltation of D o = 60 µm at which the threshold friction velocity has a minimum (Marticorena and Bergametti, 1995;Spyrou et al., 2010).It assumes that all erodible regions contain particles of size D o so that saltation is initiated when the friction velocity exceeds the threshold at D o .Thus, the threshold friction velocity and the horizontal flux H are calculated without a dependency on soil particle size.To estimate the dust vertical flux V distributed over the different size modes, the particle sizes at the sources are assumed to follow a tri-modal distribution based on the "background" modes suggested by D'Almeida (1987) and used in Zender et al. (2003).This source size distribution is considered globally uniform and the parameters of the distribution are as follows: mass median diameters D v = 0.832, 4.82, 19.38 and geometric standard deviation σ g = 2.10, 1.9, and 1.6, followed by the mass fraction of each mode: m i = 0.036, 0.957, and 0.007, respectively.The size distribution of the soil elements is different from the size distribution of the particles set in motion into the atmosphere.Once the particles are emitted into the air, the particle size distribution is adjusted to the 8 size bin modes from 0.2 to 20 µm in diameter after Perez et al. (2006).Because the mass in the source modes is log-normally distributed, the mass fraction overlap M ij of each source mode i, carried in each transport bin j , is calculated with the use of the standard error function (Schulz et al., 1998;Zender et al., 2003): where D min and D max are the minimum and maximum diameters of each transport bin j , and D v is the mass median diameter of the source i.
Finalizing the formulation of the dust emissions scheme DU1, the vertical flux (kg m −2 s −1 ) for each transport bin j , is calculated using: where F 1 is an empirical conversion factor (10 −4 for DU1), B is the fraction of bare soil exposed in a grid cell (depends on the percentage of each Olson biome, the fraction of land covered by water or snow and the fraction of ground covered by vegetation; the vegetation area index is used to calculate the monthly fraction of ground covered by vegetation), α is the sandblasting efficiency, H is the horizontal flux and 3 i=1 m i M i,j is the mass fraction of each transport bin.The empirical factor F was chosen following Zender et al. (2003) where a similar tuning factor was applied (7 × 10 −4 ).The choice of the particular F for each formulation was based on providing reasonable global emissions of dust.
In the second scheme DU2 an explicit size distribution of the soil is introduced, based on the Zobler soil type categorization, thus representing each grid cell by fractions of defined particle sizes (7 types from coarse to fine organized in 4 modes as shown in Tegen et al. (2002), and shown in Table 3).The mass fraction mf i of each size population is also listed in Table 3.As in DU1, the transport size distribution is adjusted to 8-size bins from 0.2 to 20 µm in diameter.The horizontal dust flux is calculated with Eq. ( 9) for each diameter of the soil particles (ranging from 0.1 to 1000 µm).Consequently, the vertical flux V for every transport bin is estimated as follows: Table 3. Zobler soil texture classification and mass fraction for each soil size population after Tegen et al. (2002).For all particle sizes the geometric standard deviation σ g is 2.0.where The parameters in Eq. ( 13) are the same as described in Eq. ( 12), with the exception of the empirical conversion factor (F 2 = 10 −3 for DU2) and mf i , which is the mass fraction of each of the four source sizes in Table 3.The vertical flux in DU2 includes the soil size distribution dependence of the horizontal flux H (Eq. 9) and the assumption that the breakage of the soil aggregates that result in the atmospheric suspension of dust particles (sandblasting) is related to the soil size distribution.The emission scheme of Marticorena and Bergametti (1995) cannot represent the sandblasting procedure explicitly (Grini and Zender, 2004;Grini et al., 2002) and the assignment of the emitted size distribution is investigated by attributing this to the soil properties whereas in DU1 the d'Almeida particle size distribution is imposed uniformly.
For consistency with the aerosol module in the EMAC model, the 8 transport size bins are grouped into the accumulation and coarse insoluble modes used for all aerosol physical and chemical processes (Aitken mode particles are not produced by the dust scheme).The aerosol module used in this work is M7 (Vignati et al., 2004), using 7 size and solubility modes for the aerosol distribution.The freshly emitted dust particles are assumed to be initially insoluble, with a geometric standard deviation and mass median diameter σ g = 2, mmd = 3.5 µm for the coarse mode and σ g = 1.59, mmd = 0.7 µm for the accumulation mode (Cheng et al., 2008).The dust particle size distribution used in this work has not been assessed except indirectly (by the mass concentration and AOD at distant locations) and will be the subject of future work.The differences between the two versions of the dust emission scheme (DU1 and DU2) are summarized in Table 4.

Observational data
For the comparison with model results, available observational datasets of monthly and annual dust concentrations and deposition, aerosol optical depth from sun-photometer data include those of the Aerosol Robotic Network (AERONET), and additionally from the MODIS-Terra (v5.1 -Level 3 product), MISR-Terra (v.31 -Level 3 product) satellite instruments and Deep Blue algorithm (which are included in the Supplement).The concentration and deposition datasets are taken from the AeroCom benchmark dataset presented by Huneeus et al. (2011) (N. Huneeus, personal communication, 2011).More specifically, the comparison is performed for monthly and annually averaged dust concentrations based on measurements at 24 stations, of which 22 are managed by the Rosenstiel School of Marine and Atmospheric Science, University of Miami (Prospero et al., 1989;Prospero, 1996;Arimoto et al., 1995).Two additional stations, Jabirun and Ruchomechi, are added as in Huneeus et al. (2011).All measurement stations are located downwind though remote from the main dust source areas, and the measurement period covers the 1980s and 1990s, while each station has been active during different periods.
The 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 as described in Huneeus et al. (2011).Figure 3a shows the locations of the 24 stations; the names and coordinates are given in Table S1 of the Supplement.From this dataset we also use the annual average from each station to evaluate the model calculated annual and seasonal dust distributions.These measurements are multi-annual and not for the simulation year 2000, which will be taken into account in the model evaluation in Sect.5.1.Furthermore, monthly dust concentration measurements for 2000 are available for Miami and Barbados (J.Prospero, personal communication, 2010) and for Tel Shikmona (Haifa, Israel) (B. Herut, personal communication, 2012).
The deposition of dust particles directly relates to their distribution and mass load in the atmosphere.Dust particles are mostly deposited through precipitation scavenging and sedimentation and to a lesser extent through dry deposition.Annual deposition data from 84 sites (names and locations are listed in Table S2 of the Supplement) are used to evaluate the modelled dust deposition.This dataset is a collection  from different sources, the measurement period being different for the stations; for more details we refer to the indicated publications.The 84 stations are shown in Fig. 3b, colourcoded according to region.The regions are: Europe, Asia, Nand S-Atlantic, S-and W-Pacific, S-Ocean and Indian Ocean (Table S2, Supplement).Further, the monthly dust masses deposited in three locations in Florida are used (Florida Atmospheric Mercury Study-FAMS network, Prospero et al., 2010).The deposition measurements for each site have been averaged by month over three years (1994)(1995)(1996), providing the total and the wet deposition fraction.More details about these datasets can be found in Huneeus et al. (2011).
The AOD observed by the AERONET stations (level 2 data) is used to compare with the model output for the year 2000 on a daily, monthly and annual basis.Since the AERONET Level 2AOD is not given at 500 nm or 550 nm for all the selected stations, the measured AOD at 550 nm is obtained by an interpolation method using the measurements at 440 and 870 nm (de Meij and Lelieveld, 2011).The AERONET stations provide AOD data during daytime and clear sky conditions, and the daily average of the available measurements is used for the evaluation.
The selection of the AERONET stations included in this evaluation presupposes that they are located in areas influenced mainly by dust sources, thus avoiding highly polluted regions.To select the appropriate stations for the year 2000, two criteria were applied to the daily AOD of each available station.The first criterion is that the AOD 550 nm exceeds 0.2 and the second that the Angstrom exponent (AE 500−870 nm ) is less than 1.2 to select the coarser particles associated with desert dust transport (Kinne et al., 2003;Dubovic et al., 2000).The station is considered "dusty" if at least 20 % of the available data satisfies the above criteria.This methodology has provided 19 available stations for the simulation year 2000.A list is presented in Table S3 of the Supplement and the location of each station is shown in Fig. 3d.For some stations the availability of the AERONET aerosol optical depth has been limited and in these cases monthly AOD values from MODIS-Terra, MISR and MODIS Deep Blue (Level 3) have been included to evaluate the model (included in the Supplement).Of course, there are limitations to the use of level-3 data, which are averages from level-2 pixels on a 1 × 1 degree resolution.The sampling of actual retrievals is highly non-uniform in space and time, even at the resolution of these products (MODIS 1 × 1 degree and MISR 0.5 × 0.5 degree) (Kahn et al., 2009).The use of the satellite data is complementary to that of the AERONET data, since it is recommended to cautiously compare the monthly level 3 data products with the AERONET data (Leptoukh, 2011;Levy et al., 2009).

Effects of nudging to meteorological analyses
As discussed by Timmreck and Schulz (2004), nudging a simulation with observed, i.e., analysed meteorological data compared to a free-running general circulation model, can substantially affect the simulated global dust budget and the seasonal variability of concentration distributions.In this work, the ERA40 reanalysis data has been assimilated into EMAC for the nudged simulation, and as a second option the free-running model was used by applying multi-annual mean sea surface temperatures and ice coverage as boundary conditions.The latter runs were included to study possible artefacts by the nudging and differences relevant for climate change applications.The simulations are denoted as DU1 and DU2 for the free-running model and DU1 ERA40 and DU2 ERA40 for the nudged model.At first glance the differences in the annual emissions from running the model in these two modes are hardly discernible.The annual emissions from DU1 ERA40 and DU2 ERA40 are presented in Fig. 4. The geographical distribution is largely the same, except for some areas where the emission flux appears to be lower in the nudged simulation (N-Africa, Arabian Peninsula).By calculating the global emission fluxes the differences become more apparent (Table 5).The differences in the seasonal emissions between nudged and free-running mode and among the DU1 and DU2 formulations appear in the Supplement (Figs.S1 and S2).
The nudged simulations (DU1 ERA40 and DU2 ERA40) produce less dust globally compared to the free-running model, in agreement with the results of Timmreck and Schulz (2004) for the ECHAM4 climate model that the nudging somewhat reduces the wind speed and the dust emissions accordingly.For DU1, the reduction of the annual emissions is ∼ 20 % and for DU2 it is ∼ 24 %.The atmospheric lifetime of airborne dust is slightly longer in the nudged version of the model.The annual deposition reduces in line with the emissions between 16 and 25 %, and the annual atmospheric loads are also reduced compared to the free-running model (−14 % for DU1 and −19 % for DU2).The parameterization of the dust emissions is such that even small differences in the friction velocity can have substantial effects (Eqs.8 and 9 include a dependency of the horizontal flux on the friction velocity to the power of three).The same is evident from the regional analysis of the dust budgets in Table 6.
The nudged simulations produce lower emissions compared to the free running model in all areas except Australia, because the wind speeds are less affected there by the nudging (Fig. S1 in the Supplement demonstrates the seasonal difference in the emissions for DU1 and DU1 ERA40).
In spite of the effects of the nudging on the dust emission strength, its importance for a direct comparison with observations is illustrated in Fig. 5, showing model results and measurement data at the station Miami (courtesy J. Prospero).This station is affected by dust transport from the Sahara, predominantly in summer, whereas during winter the prevailing westerly winds minimize the transported transatlantic Saharan dust.The lower panel shows the multi-year mean simulation where the seasonal cycle is captured by the free running model.The modelled dust concentration is within the variability of the measurements with r = 0.96, mean bias = −0.64,RMSE = 2.8 for DU1 and r = 0.95, mean bias = −0.59 and RMSE = 3.03 for DU2.The upper panel shows the comparison for the year 2000, indicating a bimodal seasonal profile, thus deviating from the multiannual mean.Again the model captures the seasonality of the observations, while underestimating the concentrations in the months September to November.During these months, the transport of desert dust from N-Africa is weak (as shown by satellite images and sun-photometer data).The differences between model and observations are related to the simulation of the removal processes, the inadequate northward transport of Saharan dust during these months or the existence of local sources that are not included in the model.Nevertheless, the correlation between modelled and observed concentrations is r = 0.91 for DU1 ERA40 and r = 0.90 for DU2 ERA40 with the regression line close to a 1 : 1 relationship.The mean bias is −2.8 and the RMSE is 3.7 for DU1 ERA40 and the mean bias is −3.1 and the RMSE is 3.8 for DU2 ERA40.

Dust concentrations
Figure 6 presents a comparison of the annual mean dust concentrations calculated for the year 2000 and the observed multi-annual means at 24 stations.The model appears to simulate the spatial variability within a factor of 10 with both emission parameterization schemes DU1 and DU2.The colours of the symbols in Fig. 6 correspond to the station locations shown in Fig. 3a, to help distinguish the geographical areas.The same colour coding is used throughout this section.January is excluded from the comparison to avoid differences due to the model initialization.The comparison for stations in the Atlantic region (Barbados, Izana, Bermuda, Miami and Mace Head shown in green) shows a relationship between modelled and measured values close to 1 : 1 for all simulations (DU1, DU2, DU1 ERA40 and DU2 ERA40).This demonstrates the ability of the model to simulate the transport of dust from the main source areas in N-Africa.The simulations for the Asian stations (Jeju and Hedo in pink) are also in agreement with the observed annual means, with the highest correlation for the DU2 and DU2 ERA40 model versions.By also considering Midway Island (no. 12 in Fig. 3a, orange square with an x-symbol in Fig. 6) as a sources that are not represented in global models.Finally, the dust observed in Australian stations (yellow), i.e., Cape Grim and Jabirun, is underestimated by all simulations, possibly related to the underestimation of the Australian dust sources.
The statistical analysis, comparing the annual average dust concentrations from multi-annual observations to the model results for the year 2000 (Table 7) shows correlations in the range 0.84-0.88 and a relatively low bias, especially for the nudged simulations.Also, the normalized mean bias is low for the nudged simulations, especially for DU1 ERA40 (1.5 %).The root mean square error (RMSE) is lowest for the nudged simulations, for which the slope of the linear regression is closest to one.The difference between the simulations is generally small, though the nudged runs slightly outperform the other model simulations.
Extending the comparison for these 24 stations to the monthly average dust concentrations, the general picture is similar (Fig. 7).Again, January has been excluded.The correlations of the simulation results for the Atlantic Ocean stations (green) with the observations in all simulations show   Dust concentration measurements for the year 2000 were kindly provided for 3 stations (Barbados, Miami, USA, and Tel Shikmona near Haifa, Israel) by J. Prospero (personal communication, 2010) andB. Herut (personal communication, 2012), respectively.The availability of these comprehensive datasets, together with AERONET data, has motivated our selection of the year 2000 for the simulations.Figure 8 compares the observations to the results of the nudged simulations.The two parameterizations yield similar results for Barbados and Miami (Fig. 8, upper and middle plot, respectively), with the Barbados concentrations being overestimated by the model and the Miami ones being closer to the observations during most months.In Barbados the model overestimates concentrations during the dusty summer season by a factor of 1.5 to 3. In Miami, the simulated concentrations agree with the observed values within a factor of 0.7 to 1.02 for this season, though the model underestimates concentrations in the transition seasons by a factor of 10 to 50, especially in September-October.For the Tel Shikmona station, the DU1 ERA40 scheme appears to perform better DU2 ERA40, though the differences are not large (DU1 ERA40 values are 0.6 to 1.1 times the DU2 ERA40).During two days in April (11-12) the model results are substantially higher than the observations by a factor of 2.5, i.e., with both versions of the parameterization scheme, which is related to the meteorological conditions in that period.The modelled wind speed at 10 m is similar to the measured at the WMO station in Haifa (http://www.ncdc.noaa.gov/oa/mpp/).The daily average modelled value is 2.6 m s −1 for April 11 and 5.3 m s −1 for 12 April, and the observations are 3.8 m s −1 for 11 April and 5.9 m s −1 for 12 April.The high modelled dust concentrations result from the underestimation of the precipitation during these days which underestimated the wet deposition fluxes of the dust particles.The daily precipitation rate from the TRMM online visualization and Analy-sis System (http://disc.sci.gsfc.nasa.gov/precipitation/tovas/) is 3.4 mm day −1 for 11 April and 0.02 mm day −1 for 12 April for Haifa, whereas in the model the precipitation is almost zero for both days at that location.For all three stations the correlation coefficients are rather high (0.73-0.91) and the spread in the scatter plots is low.

Dust deposition
Annual deposition measurements of dust are available for 84 stations, and monthly bulk and wet deposition dust data for three locations in Florida, as mentioned in Sect. 4. The 84 stations are shown in Fig. 3b, and are typically located downwind and partly remote from the main dust source areas.This implies that to a large degree we test the transport and deposition qualities of the model, and to a lesser extent the emission schemes, depending on the distance from the sources.The data represent multi-annual averages that unfortunately do not match the simulation year, so that the evaluation is rather qualitative.The comparison with the different simulations is shown in the scatter plots of Fig. 9, which also includes correlation coefficients, mean biases, RMSE and normalized mean biases (NMB).Again, the colour of the data points in the scatter plots relates to the location of the stations shown in Fig. 3b.
The measured deposition over Atlantic Ocean stations (green) is reproduced by the model in all four simulations within a factor of 0.5 to 2 (68-72 % of the modelled data are within 0.5 and 2 times the observations), in line with the above-described evaluation, indicating that the dust outflow from N-Africa over the Atlantic Ocean is reasonably simulated.Dust deposition in E-and W-Pacific locations (red and orange) tends to be underestimated by the model, though the DU2 improves the simulation for the E-Pacific region (the linear regression becomes y = 0.70x + 0.05, with r = 0.69, whereas for the other simulations is very low (not shown in the plots)).The modelled deposition over the S-Ocean (grey) is strongly overestimated in some locations reaching up to 100 times the observations, while indicating higher correlations with observations in the DU1 ERA40 simulation (y = 1.25x + 0.04, r = 0.95, not shown in the statistics in Fig. 9).The modelled deposition over the Indian Ocean (black) and at European stations (blue) agrees with the observations within a factor of 10, except one value in Lake Kinneret (Israel).In two locations, the Taklimakan desert in central Asia (purple) and Lake Kinneret in Israel (blue), the dust deposition is systematically underestimated by the model, with no significant differences between the simulations.The underestimation of some E-Pacific stations (red) is also persistent in all simulations, with small deviations among them.The statistical analysis of the four simulations shows that the mean bias, NMB and RMSE are lowest and the correlation highest for the DU2 simulation.Measurements of bulk and wet deposition at three locations in Florida during three consecutive years (1994)(1995)(1996) have been used to additionally evaluate the dust deposition simulations and also to assess the contribution of wet deposition.The stations are Lake Barco, Tamiami Trail and Little Crawl Key, as in Huneeus et al. (2011) (Prospero et al., 2010).The free-running simulation results are compared with the measurements in Fig. 10.The differences between the two versions of our parameterization are small, likely because deposition fluxes are largely governed by transport and rainout processes, being the same in the two versions.At Lake Barco, the model captures the seasonality and the magnitude of the deposited dust, in contrast to some of the results in the AeroCom study (Huneeus et al., 2011).For the Tamiami Trail station a time shift of one month appears in the maximum deposition flux (the model maximum is in July, not in June as indicated by the measurements), and the model overestimates the bulk deposition flux in July and August.For the southernmost station, Little Crawl Key, the model captures the seasonality though also overestimates the deposition maximum in July, being the result of a too strong wet deposition flux.Since we are comparing different model and measurement periods for those stations, we cannot expect conclusive (dis)agreement among the values.Nevertheless, the model captures the seasonality of the measurements with only one month shift in one of the three stations (Tamiami Trail).Also the predominance of wet deposition is described by the model, though the maximum wet deposited amount at the two stations located further south was overestimated by almost a factor of 2. The scatter plots on the right also list the linear regressions and correlation coefficients.

Aerosol optical depth
The modelled AOD, which is dominated by dust at the selected stations, analysed at 550 nm wavelength, has been evaluated on a daily, monthly and annual basis.Although AERONET provides measurements of AOD at high temporal resolution of the order of minutes, the evaluation in this paper is based mainly on the daily and monthly averages, focussing on the seasonal dust cycle and not on specific dust events.Nevertheless, the daily averages provide a rather detailed view of the desert dust distribution.The monthly AOD is estimated from the daily values for each AERONET station, and the modelled AOD is calculated for the same days as those for which AERONET data is available.Whenever the MODIS (v5.1),MISR (v31) or MODIS Deep Blue (v.5.31)AOD is used, the monthly modelled values are calculated from all days of the month.For this evaluation the nudged simulations for DU1 and DU2 are used.The MISR and MODIS related comparison for each AERONET station is included in the Supplement as it is complementary to the comparison with the sunphotometer data.
The comparison between daily measured and modelled AOD at the 19 AERONET stations for the year 2000 is shown in Fig. 11.The mean biases are small for both versions of the emission scheme (−0.008 for DU1, and 0.014 for DU2) and the average and standard deviation are close to the measured values (Table 9).The simulations do not show statistical significant differences between the two versions of the dust scheme, although at individual locations the AOD from one version can be twice the AOD from the second version.In addition, we consider each station individually and in groups according to the location (colours according to Fig. 3d).For the African stations (no. 4, 6, 7, 9 and 16 in red), the model reproduces the daily measured AODs, e.g., in the outflow region of the dust (6: Cape Verde,  7: Dakar in Table 10).At the other three stations located in Central Africa, the model tends to underestimate the AOD.This can be attributed to a possible underestimation of both biomass burning and dust emissions in this area, the latter from the Bodélé Depression, for example, which is a major dust source region in N-Africa.In contrast to some other models, we have not tuned preferential dust sources such as the Bodélé Depression because we favour process consistency within EMAC, with the risk of under-representing such pronounced source regions (Todd et al., 2008).The daily AODs over the S-American and W-Atlantic stations (no. 2,5,14 and 19) are represented by the model with correlation coefficients in the range of 0.48 to 0.70 for DU1 ERA40 and 0.46 to 0.70 for DU2 ERA40 (Table 10).The model underestimates the AOD over the station Suriname (no. 19) in northern S-America, and to a greater extent with the DU2 ERA40 simulation, possibly associated with the under-representation of anthropogenic aerosols in the model, notably biomass burning aerosols, which may influence this station.This is supported by the fine/coarse AOD from the AERONET database (not included) which shows that the coarse mode AOD dominates the total AOD during January to May and in July.During these months the model values are close to the observed ones (Fig. S3, station 19).
The discrepancies are more pronounced during June, and September to December when the fine and coarse AOD almost equally contribute to the total AOD.The overestimation of the AOD for the month of June can be due to overestimation of the dust and the sea salt flux as this site is near the coast in Suriname.The model results for the stations in the Middle East (no. 11,15 and 17) compared to the daily measured AODs have correlation coefficients in the range 0.60-0.69,while the model overestimates some peak values in April and May by a factor of 3 to 5, over Nes Ziona and Sede Boker in Israel.At the Arabian stations (no. 3 and 18) the model performance is poorer with low correlation coefficients (0.36 to 0.39) and slope of the linear regression line 0.31 to 0.60.For the AOD over the Asian station (no.1), located in S-Korea, the model underestimates the AODs in January to March, November and December by a factor of 3 to 10, while performing better in April, May, September and October with the modelled values 0.7 to 1.2 times the observations (other months are absent from the observational dataset).
One station is located on an island in the Indian Ocean (no.12, Kaashidho at Maldives) where the model underestimates the AOD compared to the AERONET data in the period January to March, whereas the model performance is much better for the months of April to October with the exception of July (Fig. S3).This is likely related to pollution outflow from the Indian subcontinent during the dry season (Lawrence and Lelieveld, 2010), being underrepresented in the model.Studying the MODIS small mode fraction images (http://disc.sci.gsfc.nasa.gov/giovanni/overview/index.html)shows that during autumn and winter the small particles dominate the AOD (caused by anthropogenic inorganic aerosols [SO = 4 , NO − 3 and NH + 4 ]), while during the summer a mixture of fine and coarse particles is observed.This station is expected to be influenced also by sea salt particles as it is located in an island, thus the contribution in the coarse mode AOD is from both dust and sea salt particles (such contribution cannot be quantitatively derived from the AOD evaluation).Finally, for the stations located in S-Europe (no.8, 10 and 13; two in Italy and one in Spain) we obtain high correlations with the measurements for the Italian stations and a lower correlation with the Spanish data (Table 10).The AODs over the Italian stations (Lampedusa and IMC Oristano) are typically overestimated by the model (both emission schemes, though better with DU1 ERA40), while in southern Spain (El Arenosillo) the modelled AODs appear to be too low (Statistical analysis in Table 10).Figure 11 summarizes these results into two global maps, showing the correlation coefficients for each of the AERONET stations, both for DU1 ERA40 and DU2 ERA40.This indicates that both versions of the emission scheme perform similarly in the Middle East, the Arabian Peninsula, S-Europe and some stations in S-America while DU1 ERA40 generally achieves the best agreement with AERONET data.
The monthly average AODs are calculated by only accounting for the days for which AERONET data are available, i.e., the same days from model results and observations.Figure 12 shows the scatter plots for the two versions of the dust emission scheme.The upper panel shows the linear and the lower panel the logarithmic relationships, colour coded by the location of each station (Fig. 3d).The monthly time-series for each station in comparison to AERONET and satellite data are presented in Fig. S3 in the Supplement.A seasonal analysis of the monthly values is included in the Supplement (Fig. S4); the data are segregated per trimester on a seasonal basis (winter: December to February, spring: March to May, summer: June to August, fall: September to November).The majority of the modelled values are within 0.5 to 2 times the observations during spring, summer and fall, whereas for the winter months the correlation is quite poor.The underestimation of the monthly AOD at the station Illorin (Nigeria) prevents a closer to 1 : 1 regression for the winter months.Again this may be related to the underrepresentation of anthropogenic aerosols in this location.Another reason for the poor correlation during winter is the inclusion of January, which is considered a spin-up month.
By grouping the stations according to their location (Fig. 3d), the modelled monthly AODs compared to the AERONET measurements (Fig. 12, lower panels) are presented for the S-American-Atlantic stations (purple), the Asian station in Anmyon (cyan, S-Korea), the Middle East (green) and the S-European stations (grey).The AOD at the African stations (red) is reproduced by both versions of the model, except for Illorin, as mentioned above.Poorer agreement is obtained for the Arabian stations (blue) and the Indian Ocean station (white, Maldives).The majority of the points, though, lay within the 1 : 2 and 2 : 1 lines (Fig. 12); 70 % of the DU1 ERA40 and 67.5 % of the DU2 ERA40 monthly modelled AOD values are within the range 0.5 to 2 times the observations.Similarly, the annual AOD evaluation in Fig. 13 is coloured by the location of the stations.It should be emphasized, however, that many of the AERONET measurement time series are incomplete; hence the "monthly" and "annual" data should be interpreted with care.Similar agreement with the AOD measurements is obtained for the two versions of our dust emission scheme for the African, Middle Eastern, S-European and S-American-Atlantic stations.Differences are more pronounced for the Arabian stations and Lampedusa (Italy), indicating better agreement using the DU1 scheme for Lampedusa and the Arabian stations.
The column aerosol mass burden (µg cm −2 ) and the AOD from the model and that derived from the MODIS-Terra satellite measurements (http://disc.sci.gsfc.nasa.gov/giovanni/) are shown for the month June 2000 in Fig. 14.We selected this month because of the generally intense dust activity in the large, arid areas such as N-Africa and the ensuing dust transport over the Atlantic Ocean.The images of the mass burden only provide a qualitative evaluation, since this MODIS satellite data product is not validated to the same extent as the AOD.The outflow of desert dust towards the Atlantic Ocean appears to be well-described using both emission schemes, matching both the spatial distribution pattern and the magnitude (Fig. 14).The column mass and the AOD over the Indian Ocean and India is also in good agreement with the satellite-retrieved data, with DU2 ERA40 being closer to the observed column burden than DU1 ERA40.
The N-African deserts are excluded from the MODIS data of the columnar aerosol mass but they are present in the AOD maps (Fig. 14f).The AOD is overestimated by the model in the Mauritania and West Africa region and also over the Arabian Peninsula with both schemes, with DU2 producing more dust in those areas than DU1.Further, the DU2 scheme seems to overestimate the dust burden from the Asian deserts.Finally, the N-American arid areas are represented by the DU1 scheme whereas DU2 again overestimates the dust load.Both schemes overestimate the amount of dust from S-America.Notwithstanding indications in the literature and from the MODIS data that S-America is only a weak dust source, especially compared to N-Africa (Prospero et al., 2002), both versions of the emission scheme generate significant dust plumes, while DU1 generates less dust than DU2.This is possibly related to the coarse resolution of the model, which smoothens the pronounced terrain, leading to too high friction velocities over the Patagonian desert.Again, this is also a consequence of applying one consistent parameterization throughout the global domain, whereas some models apply regionally tuned emission fluxes (Cakmur et al., 2006;Miller et al., 2006).The same pattern is indicated in the AOD plots shown in Fig. S5 from the MODIS-Terra and Deep Blue instruments compared to the model results.The AOD is calculated as an average between March to November 2000, based on the availability of the MODIS data.We should note here that the modeled AOD was averaged over the entire period and there is no available information for the missing data from MODIS, thus the comparison between model and satellite retrievals must be done with caution.Nevertheless, the results from the AOD comparison indicate that both versions of the dust emission parameterization adequately describe many of the sources, and that the EMAC model reproduces the dust transport over the Atlantic Ocean and the Mediterranean region realistically.The DU2 scheme, which includes an explicit soil particle size distribution, appears to perform better at locations like Anmyon, Dakar, Bahrain and Erdemli.Both schemes underestimate atmospheric dust at Illorin, possibly caused by too low emissions expected that our model somewhat underestimates the AOD as compared to remote sensing observations.

Discussion and conclusion
In this paper, two formulations of a process-based dust emission scheme have been introduced in the atmospheric chemistry general circulation model EMAC; one in which the soil properties are assumed globally uniform and the other in which soil properties are assigned from a gridded database.The emission scheme follows to a large degree the work of Marticorena et al. (1997) as well as other published work.
The effect of assigning a soil size distribution in each grid cell compared to a uniform pattern is investigated and also the effect of nudging the simulations towards observed meteorological data.The simulations provide encouraging results both for the spatial and temporal distributions of dust particles on a global scale.The new formulation includes more realistic spatial and temporal heterogeneity of the emitted dust, and the online coupling with the meteorology and the soil characteristics is more consistent.Furthermore, this formulation allows for the emission of dust particles of variable sizes (accumulation and coarse), which are linked to the aerosol module of the EMAC model, and can be expanded in the future.
The two versions of the dust emission scheme are primarily different in the explicit representation of soil particle size distributions and also in the emitted size distribution at the source.Whereas in DU1 the (tri-modal) particle size distribution at the source is globally uniform, in DU2 it is explicitly accounted for based on the Zobler geographical soil texture classification and four soil populations, as listed in Table 3 (Tegen et al., 2002).This influences the threshold friction velocity, which triggers the dust mobilization and hence influences the emission fluxes.It should be stressed that the Zobler soil texture classification has been derived using wet sedimentation measurement techniques, which break the soil aggregates.This increases the number of free clay particles, thus underestimating the number of large size aggregates (Laurent et al., 2010;Kok, 2011).Furthermore, Laurent et al. (2006) mention that there is no direct relationship between the soil grain size distribution and the soil texture in the Northeast Asian deserts.An advanced technique, based on dry sedimentation, has been used by Chatenet et al. (1996), followed by Laurent et al. (2006Laurent et al. ( , 2008Laurent et al. ( , 2010) ) though the measurements are limited to N-African, Arabian and NE Asian deserts and cannot be applied within a global framework at this stage.Nevertheless, making use of these data is planned as a next step in our work.
The model calculated atmospheric dust budgets based on the DU1 and DU2 dust emission schemes differ substantially (Table 5).We have to note here that the calculated budget is simply a model diagnostic tool and not an indication of quality.These values cannot be directly compared to measured/observed global quantities that would allow a proper characterization of the simulated budgets.In the free-running model the global source with DU2 is about 1000 Tg yr −1 stronger than with DU1, and in the nudged simulations the difference is about 713 Tg yr −1 .It appears that the results of the nudged DU1 ERA40 simulation are closest to the median source of the AeroCom exercise, i.e., 1123 Tg yr −1 (Huneeus et al., 2011).The DU2 scheme produces stronger emissions than DU1, mostly due to differences in the Asian and S-American deserts and to a lesser extent in N-Africa (Table 6).The stronger sources by DU2 primarily increase the dust loads over the source regions, and it would be useful to have access to additional measurement data there.For the N-African deserts the total dust emissions by the two versions do not deviate much except for DU2 ERA40, which seems relatively low (460 Tg yr −1 ).The two versions also produce similar emissions in the Middle East and Australia.For S-Africa, Asia, N-and S-America, on the other hand, the differences can be a factor of two to three.This is a result of the substantially different soil particle size distributions and emitted size distributions.The annual cycles by the DU1 and DU2 schemes are quite similar, also because the seasonality is predominantly determined by the meteorology rather than the soil classification (Fig. S2 in the Supplement shows the difference in the seasonal emissions between the two versions).The differences are mostly regional and related to the threshold friction velocity (i.e., the particle size distribution, soil moisture, drag partition correction).Over Africa, the geographical patterns of DU1 and DU2 can differ also, the latter emitting less dust from the Sahara, Mauritania and the Bodélé Depression, and more in Libya and Algeria.This is a direct effect of the size distribution assigned to the soils in DU2, because in parts of the Sahara, the Bodélé Depression and Mauritania deserts typically have relatively coarse particles, while in Libya and Algeria more medium size particles are found (according to the Zobler classification).
Our evaluation of the concentrations, deposition fluxes and the AOD does not provide conclusive evidence about quality differences between the two versions of our dust scheme.Even though the explicit soil particle size distribution in DU2 is considered more realistic, the simpler DU1 scheme appears to perform better in several locations.The general conclusion from our evaluation is that DU1 performs slightly better in reproducing the remotely sensed AOD for the year 2000 (higher correlation coefficient and slope and lower bias in the daily and monthly data), especially in the vicinity of the sources.DU2 leads to more realistic results in simulating the deposition fluxes in remote locations (higher correlation coefficient and lower bias).However, the two versions of the dust emission scheme do not show statistically significant differences in their performance, when compared to the dust concentration and deposition measurements.Future work, in which we aim to account for all classes of aerosols simultaneously and improve the representation of chemical "ageing" by dust particles in the atmosphere, will provide additional information to help evaluate and further improve our dust emission parameterization.

Fig. 1 .
Fig. 1.Left plot: Olson global ecosystem biomes.The indices on the colour bar correspond to the values in Table 2. Right plot: clay content of the soil (%).Both plots at 1 • × 1 • resolution.
Fig. 2. (a) Dependence of the threshold friction velocity on the diameter of the soil particle.The blue dots correspond to the correction of the threshold friction velocity from the drag partition scheme when f drag is constant.(b) Dependence of the drag correction parameter on the aerodynamic roughness length z o (left), and the smooth roughness length z os (right).

Fig. 3 .
Fig. 3. Stations used for the model evaluation: (a) for dust concentrations, (b) for dust deposition.The names of the stations that correspond to each number are given in the supplement.(c) The black boxes denote the areas for the calculation of the regional emissions in Table 6.(d) Location of the 19 AERONET stations used for the model evaluation of AOD.

Fig. 6 .
Fig. 6.Annual mean dust concentrations from the simulation of the year 2000 compared to measured multi-annual means at 24 stations.The colours correspond to the location of each station, as shown in Figure 3a.E-Pacific = red, W-Pacific = orange, S-Africa = blue, Atlantic = green, Australia = yellow, Asia = pink, S-Ocean = grey (the orange square with the x symbol is for the Midway Island station).The dotted lines denote the 1 : 10 to 10 : 1 range.

Fig. 7 .
Fig. 7. Comparison of modelled monthly with measured multi-annual dust concentrations at 24 stations (January excluded).The colours correspond to the location of each station, as shown in Fig. 3a.E-Pacific = red, W-Pacific = orange, S-Africa = blue, Atlantic = green, Australia = yellow, Asia = pink, S-Ocean = grey.The dotted lines denote the 1 : 10 to 10 : 1 range.

Fig. 8 .
Fig. 8.Comparison of modelled and measured dust concentrations for the year 2000 for stations in Barbados (upper panel), Miami (middle) and Tel Skikmona, Haifa (lower panel).Measurements for the year 2000 are indicated in blue; the multi-annual average concentrations are shown by the blue dotted line.The model results from DU1 ERA40 are shown in black and from DU2 ERA40 in green (nudged simulations).The scatter plots on the right also list the linear regressions and correlation coefficients.

Fig. 9 .
Fig. 9. Comparison of modelled and measured annual dust deposition (g m −2 yr −1 ) from 84 stations.The stations are shown in Figure 3b and the colour of the dots corresponds to the location of each station.E-Pacific = red, W-Pacific = orange, S-Africa = blue, Atlantic = green, Australia = yellow, Asia = pink, S-Ocean = grey.The dotted lines denote the 1 : 10 to 10 : 1 range.

Fig. 11 .
Fig. 11.Scatter plot of the modelled versus measured daily AOD550nm from the 19 AERONET stations for the DU1 ERA40 and the DU2 ERA40 simulations.The dotted lines correspond to 1 : 2 and 2 : 1 relationship.The two lower panels show the correlation coefficients for each of the 19 stations for the DU1 ERA40 (left) and DU2 ERA40 (right) simulations (daily average).

Fig. 13 .
Fig.13.Scatter plots of the modelled versus AERONET measured annual AOD550nm for the DU1 ERA40 (upper plot) and for the DU2 ERA40 (lower plot) simulations.The colour of the dots corresponds to the location of the stations: Africa = red, Indian Ocean = white, Europe = grey, Middle East = green, S-America = purple, Arabian Peninsula = blue, Asia = cyan.The dotted lines denote the 1 : 2 to 2 : 1 range.

Table 1 .
EMAC sub-models used in this study.

Table 2 .
Olson ecosystem biomes selected for the dust emission scheme.

Table 4 .
Characteristics of the two versions of the dust emission scheme.

Table 7 .
Statistics of the annual average dust concentrations.January has been excluded.
Head, Bermuda and Miami, where the modelled dust concentrations are significantly too low.This also applies to July in Mace Head, and to September and November in Bermuda.The same applies to the monthly mean model results at the Asian stations (pink), in line with the annual averages.The −3 ).The comparison for the S-Ocean stations indicates best agreement between the DU1 ERA40 results and the observations, though an overestimation is evident for all simulations.The statistical analysis shows that the comparison is generally best for the nudged simulations with www.atmos-chem-phys.net/12/11057/2012/Atmos.Chem.Phys., 12, 11057-11083, 2012

Table 8 .
Statistics of the monthly average dust concentrations.January has been excluded.

Table 9 .
Statistics of the daily average AOD for 19 AERONET stations.