Mineral dust cycle in the Multiscale Online Nonhydrostatic AtmospheRe CHemistry model (MONARCH) Version 2.0

We present the dust module in the Multiscale Online Non-hydrostatic AtmospheRe CHemistry model (MONARCH) Version 2.0, a chemical weather prediction system that can be used for regional and global modeling at a range of resolutions. The representations of dust processes in MONARCH were upgraded with a focus on dust emission (emission parameterizations, entrainment thresholds, considerations of soil moisture and surface cover), lower boundary conditions (roughness, potential dust sources), and dust–radiation interactions. MONARCH now allows modeling of global and regional mineral dust cycles 5 using fundamentally different paradigms, ranging from strongly simplified to physics-based parameterizations. We present a detailed description of these updates along with four global benchmark simulations, which use conceptually different dust emission parameterizations, and we evaluate the simulations against observations of dust optical depth. We determine key dust parameters, such as global annual emission/deposition flux, dust loading, dust optical depth, mass-extinction efficiency, singlescattering albedo, direct radiative effects. For dust particle diameters up to 20 μm, the total annual dust emission and deposition 10 fluxes obtained with our four experiments range between about 3,500 and 6,000 Tg, which largely depend upon differences in the emitted size distribution. Considering ellipsoidal particle shapes and dust refractive indices that account for size-resolved mineralogy, we estimate the global total (longwave and shortwave) dust direct radiative effect (DRE) at the surface to range between about −0.90 and −0.63 W m−2 and at the top of the atmosphere between −0.20 and −0.28 W m−2. Our evaluation demonstrates that MONARCH is able to reproduce key features of the spatio-temporal variability of the global dust cycle with 15 important and insightful differences between the different configurations.

under favorable atmospheric conditions and on long time scales, provided there is a sufficient supply of loose dust particles at Thompson Thompson et al. (2008) WSM6 Hong and Lim (2006) Radiation RRTMG Iacono et al. (2008) GFDL Fels and Schwarzkopf (1975) geometric standard deviation of 2 (Schulz et al., 1998;Zender et al., 2003).
Our new developments presented below have mostly focused on aspect (1): dust generation and uplift by surface wind and turbulence. In particular, we have implemented and tested a variety of dust emission and drag partition parameterizations, along with new datasets for dust source areas, source type (i.e. natural and anthropogenic), surface roughness, and vegetation.
Additional upgrades include the option to calculate dust extinction assuming non-spherical particle shape, as well as new 130 diagnostic capabilities (output of 3-dimensional single-scattering albedo and extinction, clear-sky aerosol optical depth (clearsky AOD) and AOD at satellite overpass times). In the following, we present the MONARCH dust module. We first describe the Dust emission scheme/Reference Abbreviation Approach Marticorena and Bergametti (1995) MB95 Dust emission based on saltation flux and soil texture Ginoux et al. (2001, modified) G01-U/G01-UST Dust emission based on a topographic dust source function Shao (2001) S01 Dust emission based on volume removal by saltation Shao (2004) S04 Dust emission based on volume removal by saltation (parameterized saltation bombardment efficiency) Shao et al. (2011) S11 Dust emission based on volume removal by saltation (reduced form) Kok et al. (2014b) K14 Dust emission based on brittle fragmentation by saltation treatment of dust emission, summarizing previous and detailing new developments. Then, we recapitulate the implementation of dust transport and deposition, and interactions with radiation.

Dust emission and lower boundary conditions
135 Several different parameterizations of dust emission are available in MONARCH, which cover different paradigms and range from more simplified to more physics-based descriptions. To describe dust emission generated by saltation, MONARCH includes the parameterizations from Marticorena and Bergametti (1995) (MB95), Ginoux et al. (2001) with modifications detailed below (G01), Shao (2001Shao ( , 2004) (S01, S04), Shao et al. (2011, Eq. 34) (S11), and Kok et al. (2014b) (K14). Including such a large number of dust emission schemes in MONARCH gives us the unique opportunity to directly compare them in the exact 140 same framework (e.g. meteorology, land-surface conditions) of a global model. As by nature none of the schemes can predict dust emission perfectly, discrepancies between results obtained with the different schemes can in comparison with observations provide insights into aspects of parameterizing dust emission, which are particularly uncertain (or not) and which may need more attention in future research. The six available dust emission schemes in MONARCH are summarized in Tab. 2.

Dust emission flux 145
In saltation-based dust emission schemes, the vertical dust emission flux F depends on the horizontal flux of saltating soil particles or particle aggregates. In the MB95 scheme, F is directly proportional to the total streamwise saltation flux Q, In MONARCH, the vertical-to-horizontal-flux ratio α q can either depend on the clay content of the parent soil as originally proposed in Marticorena and Bergametti (1995), or on the soil texture as proposed in Tegen et al. (2002) and described in 150 Pérez et al. (2011). In the latter case, which is the default in MONARCH, α q is determined as a mass-weighted average of the vertical-to-horizontal-flux ratios of four soil particle size classes with mean diameters 2, 15, 160, and 710 µm. S is a globally variable dust source scaling function, which was not part of the original formulation in Marticorena and Bergametti (1995), but which is introduced here as it was found to lead to improved results (Pérez et al., 2011, see also Sec. 3.1.6). The dust emission flux resulting from Eq. (1) is a bulk flux, which we distribute across particle sizes using a predefined particle size-distribution 155 (Sec. 3.1.5).
The G01 dust emission scheme does not include an explicit formulation of Q. It seeks to avoid the need for detailed descriptions of soil characteristics and instead introduces a topography-based dust source function, S, representing the availability of sediment. The dust emission flux is originally obtained as 160 where C G01 is a dimensional factor, s bare is the bare soil fraction (see Sec. 3.1.4), u 10m is the 10 m wind speed, and u t is a threshold wind speed below which no dust emission occurs (Sec. 3.1.3). Note that s bare was included in S in the original formulation. We apply s bare to all schemes in MONARCH. The bulk dust emission flux F is distributed across particle size classes using predefined fractions s p (Sec. 3.1.5). To ease comparison with other schemes, we also implemented a modified version of the G01 scheme, which estimates F using friction velocity and threshold friction velocity, u * and u * t (G01-UST), 165 instead of u 10m and u t (G01-U), such that In both implementations, G01-U and G01-UST, we introduced additional modifications on, respectively, u t and u * t , and on s p as described in Secs. 3.1.3 and 3.1.5.
The S01 scheme is a physics-based dust emission scheme, which calculates size-resolved dust emission based on the soil 170 volume removed by impacting saltation particles and explicitly considers aggregate disintegration as a dust emission process in addition to saltation bombardment. The emission of dust particles of size d i by saltation particles of size d s is given by (Shao, 2001, Eq. 52) dust fraction follows as η ci = η f i − η mi . The γ-function (Eq. 5) and therein the parameter κ determine how easily a soil is disaggregated (Shao et al., 2011;Klose et al., 2019, see also Sec. 3.1.2). Here, we use κ = 1 globally. A spatially variable, for example soil-texture dependent, κ could be easily implemented, if future investigations support such a dependency. The 185 emission flux of dust particles with diameter d i , i.e. for all saltation particle sizes, is obtained as with d max = 20 µm in MONARCH.
The S04 scheme is a simplification of the S01 scheme in which the saltation bombardment efficiency, σ m = m Ω /m ps , is approximated as (Shao, 2004, Eq. 11) 190 σ m = 12u 2 * with soil plastic pressure P . The larger u * , the more soil mass is ejected by a saltation particle impact for a given soil. The dust emission flux is given by (Shao, 2004, Eq. 6) and F (d i ) follows from Eq. 4 assuming η f i ≈ η ci . Based on a detailed comparison with field measurements, a basic version 195 of the scheme (denoted here as S11) was suggested by Shao et al. (2011, Eq. 34), which makes use of the total (instead of size-resolved) saltation flux: F S11 (d i ) = s bare c y η mi (1 + σ m ) gQ u 2 * .
Shao et al. note, however, that this simplification may be specific to the experimental data set, which had a narrow soil PSD.
The K14 dust emission scheme uses the concept of the fragmentation of brittle material (Kolmogorov, 1941;Åström, 2006). 200 It is also a physics-based dust emission scheme that includes a dynamical dependency of soil erodibility on threshold friction velocity. Although the kinetic energy supplied by saltating particles is taken into account in the scheme, it does not include Q explicitly. The K14 dust emission flux is given as where f clay is the clay fraction (from STATSGO-FAO inventory, see Sec. 3.1.6), ρ a is air density, u * st = u * t ρ a /ρ a0 with ρ a0 = 1.225 kg m −3 is a standardized threshold friction velocity, C α is a constant coefficient, and C e is a u * st -dependent coefficient representing soil erodibility.

Saltation flux
For the schemes that contain explicit representations of the saltation flux (MB95, S01, S04, S11), the saltation flux of particles 210 with diameter d s , Q s , is calculated following Kawamura (1964) (same as White, 1979) as where c Q is a coefficient, u * t (d s ) the threshold friction velocity for particles with diameter d s , and u * the friction velocity for the bare surface. In Eq. (12), the saltation of particles of different sizes is treated independently. For a soil that consists of a mixture of different sized loose particles of sufficient availability, particle impacts can cause saltation in a wider size range 215 than it would be expected based on u * t (d s ) (Ungar and Haff, 1987;Martin and Kok, 2019). In the MB95 implementation, the total saltation flux Q is used and obtained as a weighted average taking into account the relative surface area of particles in four size classes (see Sec 3.1.1 for mean diameters) as a function of soil texture (Pérez et al., 2011, Eq. 2). The S11 scheme is also based on the total saltation flux, but takes a different approach and obtains Q by weighting Q s with the particle-size distribution estimated for airborne sediment, p s (d s ), as The γ-function (Eq. 5) determines how rapidly p s approaches p f with increasing u * , i.e. how easily soil aggregates are disintegrated (Shao et al., 2011;Klose et al., 2019). Both, p m and p f , are estimated for each soil texture class as a combination of up to four lognormal distributions. The coefficients used for those distributions are given in Tab. 3. PSDs are calculated with 225 60 size-bins distributed logarithmically using a quarter-ϕ scale (Krumbein, 1934(Krumbein, , 1938 with reference diameter 2000 µm. The S01 and S04 schemes directly use the spectral, i.e. size-resolved, saltation flux from Eq. (12) (cf. Sec. 3.1.1). The G01 and K14 schemes do not contain explicit formulations for saltation flux.

Threshold friction velocity and soil moisture correction
The implementation of the threshold friction velocity for ideal (dry) conditions, u * t0 (d i ), varies depending on the dust emission 230 scheme and its requirements. In the MB95 implementation, we use the relationship from Iversen and White (1982) for the four saltation size classes (cf. Sec. 3.1.2) as described in Pérez et al. (2011). The original parameterization of Ginoux et al. (2001) estimates dust emission based on 10 m wind speed instead of friction velocity (Sec. 3.1.1) and specifies a threshold wind speed for each dust size bin, which can typically be expected to be larger than for saltation-particle sizes (see Sec. 1). In combination with the relatively simple and constant distribution of soil particles across clay and silt particle-sizes (Ginoux et al., 2001), this 235 dust-size dependent threshold wind speed leads to a more variable particle-size distribution at emission. Here, we revise this implementation and specify the entrainment threshold for saltation in G01-UST as . Coefficients for minimally-dispersed particle-size distributions as assigned to the 12 USGS soil texture classes. Each PSD is composed of four lognormal distributions p1, p2, p3, and p4. Coefficients are taken from Klose (2014), 3 Different sample used as reference than in Table 6.1 of Klose (2014), but same underlying data set.
4 Sandy clay loam in Table 6.1 of Klose (2014) based on the theoretical expression for u * t0 (d i ) from Shao and Lu (2000). In the G01-U implementation, we use a fixed minimal threshold of u t d0 = 5 m s −1 , a typical 10 m-wind lower limit for dust emission under favorable land-surface conditions 240 (Kurosaki and Mikami, 2007;Pu et al., 2020). To obtain a more realistic PSD at emission in combination with a dust-particle size independent entrainment threshold, we replace the PSD described in Ginoux et al. (2001) with that of Kok (2011a) (see Sec. 3.1.5). The K14 scheme also makes use of a particle-size independent threshold friction velocity and in MONARCH, u * t d0 for K14 is obtained based on Iversen and White (1982) for 70 µm, a diameter in the size-range where u * t0 becomes minimal.
In the implementations of the S01, S04, and S11 dust emission schemes, u * t0 (d i ) is described as in Shao and Lu (2000) and 245 the minimum value (Eq. 15) is used in Eq. (5).
Models are known to underestimate the strong-wind tail of the wind speed distribution by different degrees depending on their resolution. This is particularly relevant for dust emission (Cakmur et al., 2004;Timmreck and Schulz, 2004;Cowie et al., 2015). If the frequency of occurrence of wind speeds or friction velocities above the threshold for particle entrainment is underestimated, dust emission will be underestimated, too. For coarse model resolutions (temporal or spatial), this underestimation 250 might be considerable in some regions, for example in areas with frequent moist convection or pronounced topography. In some models, this effect is mitigated by introducing sub-grid scale wind variability (e.g. Cakmur et al., 2004;Lunt and Valdes, Table 4. Summary of available options in MONARCH to account for soil moisture in the particle entrainment threshold.

255
As a result, dust emission is initiated more often and over larger areas.
When the soil is moist, the threshold friction velocity above which particles are lifted is higher than under dry conditions, because soil-water capillary forces increase the cohesion between the soil particles (Chepil, 1956;Zimon, 1982;Chen et al., 1996). This is implemented by first estimating the threshold friction velocity for dry conditions, u * tdry , and then applying a correction factor, f w > 1, to obtain the threshold friction velocity for the given (moist) conditions (McKenna Neuman and 260 Nickling, 1989;Fécan et al., 1999;McKenna Neuman, 2003;Cornelis et al., 2004a, b;Klose et al., 2014): In MONARCH, the soil moisture corrections from Belly (1964), Fécan et al. (1999), and Shao and Jung (unpublished manuscript, 2000; see Klose et al. 2014) are available in combination with all schemes. The options to account for the impact of soil moisture on dust emission in MONARCH are summarized in Tab. 4 and further detailed below.

265
The soil moisture correction from Belly (1964) is implemented as described in Ginoux et al. (2001): where θ is the volumetric soil moisture [m 3 m −3 ], f w Bwet is a large value prohibiting particle movement (here f w Bwet = 100), and c f1 is an optional calibration factor described below. This correction is used as the default for the G01 scheme. The soil moisture correction after Fécan et al. (1999) is implemented as with gravimetric soil moisture content w [%], gravimetric air-dry residual soil moisture content w r [%], and coefficients a = 1.21 and b = 0.68 (Fécan et al., 1999). w r is obtained based on Eq. (14) in Fécan et al. (1999). The conversion from volumetric soil moisture content θ to w is implemented as described by Zender et al. (2003, Eqs. 7-9): Here, ρ l is the density of water, ρ bd is the bulk density of dry soil, ρ pa is the average soil particle density (here ρ pa = 2500 kg m −3 ), θ sz is the volumetric soil moisture at saturation, and M sand is the sand fraction in the soil (Pérez et al., 2011, Tab. 1). The factor 100 converts soil moisture content from kg kg −1 into %. As the top-layer soil moisture in models is usually 280 obtained for a layer of several centimeters and is therefore typically higher than at the actual surface-atmosphere interface (which is relevant for dust emission), the soil moisture correction f w using the model's soil moisture is often too high and precludes dust emission. An optional calibration factor, c f1 or c f2 , can therefore be applied if needed. effect to reduce f w . We recommend using either c f1 or c f2 , but not both at the same time. Shao and Jung (2000, unpublished manuscript) and Klose et al. (2014) developed a soil moisture correction similar to that of Fécan et al. (1999), but based on the soil-water retention curve from Brooks and Corey (1964) rather than that from Gardner (1970). Including the optional coefficient c f1 , the correction is

290
where θ r is the volumetric air-dry residual soil moisture, h w is a function combining different constants, and ψ s is the saturation capillary pressure head . Eq. (23) is consistent with Eq. (19), as can be seen when setting Note that the volumetric soil moisture θ [m 3 m −3 ] is used in Eq. (23). θ s is the saturation (volumetric) soil moisture. The values for α, β, and θ r were obtained in Shao and Jung (2000) through fitting with observations and were published in Klose

Surface roughness, drag partition, and cover
Surface roughness through, e.g., vegetation, pebbles or rocks, absorbs momentum from the air flow and reduces the atmospheric momentum available for particle entrainment. We account for this drag partitioning using either the scheme of Raupach et al. Table 5. Summary of available options in MONARCH to account for surface roughness in particle entrainment.

Roughness correction/Reference
Description of input data Marticorena and Bergametti (1995); King et al. (2005) static roughness length (Prigent et al., 2012) and dynamic roughness length from monthly MODIS LAI (Myeni et al., 2015) Raupach et al. (1993) dynamic frontal area index from monthly vegetation cover (Guerschman et al., 2015, photosynthetic and non-photosynthetic vegetation) or AVHRR (Gutman and Ignatov, 2010, green vegetation) drag partition correction is applied to u * t , which is phenomenologically, but not physically, correct as discussed in Kok et al. (2014b). For use with all schemes in MONARCH, we apply the drag partition correction, f v < 1, on the friction velocity u * NMMB provided by the atmospheric model NMMB, such that the friction velocity acting on the erodible surface and used in Eq. (12) is In the parameterization of Raupach et al. (1993), the ratio f v between the friction velocity acting on the erodible surface and the total friction velocity supplied by the atmosphere is given as where τ is the total stress, τ s = τ s (mλ) is the maximum surface stress on the exposed area estimated from the average surface 310 shear stress on the exposed area, τ s , for a surface with lower roughness density using the constant m ≤ 1, σ v is the ratio of roughness-element basal to frontal area, and β R is the ratio of roughness-element to surface drag coefficients. Here we chose σ v = 1, β R = 200, and m = 0.5 (Shao et al., 2015). We estimate the frontal area index, λ, based on the vegetation cover fraction as (Shao et al., 1996)

315
where c λ is a coefficient. If the roughness elements are uniformly distributed and isotropically oriented, c λ = 1 Shao et al., 1996). As this is typically not the case, a value of c λ = 0.35 was proposed by Shao et al. (1996) based on measurements for stubble roughness. Stubble roughness can typically be associated with agricultural land use for which vegetation and its remains after the growing season are still relatively homogeneously distributed. An even smaller value for c λ , which leads to a weaker effect of vegetation cover in the drag partition correction, may be more appropriate for roughness 320 elements that are distributed heterogeneously, as it is typical in semi-arid regions. Here we choose c λ = 0.2. In MONARCH, λ can be estimated using Eq. (27) based on monthly satellite-based retrievals of photosynthetic and non-photosynthetic vegetation cover (PV and NPV) (Guerschman et al., 2009(Guerschman et al., , 2015, interpolated to the day of simulation (used as η in Eq. 27). Although NPV is intended to represent only vegetation components, it may also include some geological features, which is advantageous for our purposes. Monthly climatologies of the same data set (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017) and also of green vegetation cover estimated from 325 AVHRR (1985-1990, Gutman and Ignatov, 2010 are also available. Figure 1 (a-c) shows annual averages for 2012 of PV, NPV and f v R . With the parameter settings as described above, only areas in northern Africa, the Middle East, and western East Asia (Taklamakan desert) experience a low or moderate roughness correction. Areas in other parts of East Asia, Central Asia, Australia, as well as parts of North and South America and southern Africa show a stronger correction, but one which can still allow dust emission under strong wind conditions. Dust emission 330 from other areas is typically suppressed by a larger vegetation coverage using this drag partition parameterization and the given parameters. Variations in the parameters used for f v R will lead to changes in the roughness correction, particularly in areas with moderate vegetation coverage.
In the formulation from Marticorena and Bergametti (1995), f v is given by

335
where z 0 is the aerodynamic roughness length, z 0s is the smooth aerodynamic roughness length, and X is a parameter related to the distance downwind from an individual obstacle. As the surface becomes rougher (corresponding to increasing z 0 ), f v becomes smaller and the stress on the erodible surface decreases, reducing emission. The smooth roughness length z 0s is estimated as

340
where d c = 650·10 −4 cm is assumed to be the coarsest diameter of particles in the soil bed (Sherman, 1992;Pierre et al., 2014).
The aerodynamic roughness length z 0 is obtained globally in MONARCH as a combination of two different data sets. In arid regions, we use a static roughness, z 0stat , which is derived from satellite microwave backscatter (ASCAT) and visible/nearinfrared reflectances (PARASOL) (Prigent et al., 2012). In semi-arid regions, including natural vegetation and cultivated areas, we estimate a time-varying or "dynamic" roughness (z 0dyn ) based on the dimensions of green vegetation characterized using 345 the MODIS Leaf Area Index (LAI). The calculation of z 0dyn is based on empirical relationships from Marticorena et al. (2006): where h is the vegetation height and λ the roughness density (or frontal area index), defined as λ = n · a f , where n is the number density of roughness elements (number per unit area) having a frontal area a f . λ is calculated assuming patches of vegetation of diameter D η = 5 m, the number of which increases with the vegetation cover fraction η, n = η/(π · (D η /2) 2 ) 350 (Pierre et al., 2012). With a f = h · D η , it follows that scale with LAI: where h max is the maximum annual vegetation height and LAI max is the LAI above which dust emission is precluded. This approximation entails that η = 1 and h = h max for LAI = LAI max , decreasing linearly until η = 0 and h = 0 for LAI = 0. Due to the lack of data at global scale we currently assume h max = 0.4 m, a value obtained for the Sahel (Mougin et al., 1995;Pierre et al., 2012). We also set LAI max = 0.3 as in the Community Land Model (Mahowald et al., 2010;Kok et al., 2014a) although 360 this value should be further tested and constrained in future studies. We note that while η should scale with LAI at low fractional ground cover, the scaling may be weaker as leaves start overlapping, an aspect that is currently neglected in our simplified approach. In model grid cells, in which both z 0stat and z 0dyn are available, we use the larger value, z 0 = max (z 0stat , z 0dyn ).
The correction f v is smallest (i.e. roughness is largest) for roughness elements like stones or tall and closely spaced vegetation. Although Eq. (28) incorporates these dependencies, there is uncertainty related to characterizing the height and spacing 365 of roughness elements, particularly where they are of irregular size or spacing. For example, in regional studies, X has been set to 10 cm (Marticorena and Bergametti, 1995),  Marticorena and Bergametti (1995) or Pierre et al. (2014) instead, the resulting drag partition is substantially more restrictive.
Apart from the effect of vegetation or other roughness elements to absorb atmospheric momentum, they also directly prohibit particle entrainment from the area they cover. Similarly, areas covered by snow/ice (η snow ) or bedrock (η br ) preclude particle emission. We take this into account by scaling the obtained dust emission flux with s bare, in combination with the drag partition from Marticorena and Bergametti (1995) and with s bare, combination with the drag partition parameterization from Raupach et al. (1993). The area covered by vegetation is already accounted for in the latter, which determines the fraction of shear stress acting on the uncovered surface Webb et al., 2020). Alternatively, the bare soil fraction can be applied to the saltation flux. Accounting for s bare in either the dust emission flux or the saltation flux, but not both, assumes that saltation impacts eject dust close to their origin, i.e. saltation trajectories are short. This may not always be the case and saltating particles may also impact on the vegetated surface fraction 385 in a grid cell, where no emission occurs.

Particle-size distribution at emission
The particle-size distribution of emitted dust is key to quantifying the emitted dust mass, dust loading in the atmosphere, dust interactions with the energy and water cycles, along with more general impacts of dust upon climate. Whether or not the emitted dust PSD changes with the magnitude of atmospheric forces is still debated (e.g. Kok, 2011b;Shao et al., 2020). The 390 S01, S04, and S11 dust emission schemes estimate size-resolved dust emission fluxes, the PSDs of which vary with atmospheric forcing. In contrast, the K14 scheme assumes a PSD that is independent of wind speed. The G01 scheme originally distributed the estimated bulk dust emission flux across four particle-size classes (Sec. 3.1.3) and the MB95 scheme does not include assumptions of emitted dust particle sizes. For the latter two schemes, a pre-specified PSD is assigned to the estimated bulk dust emission flux that can be chosen to follow either D' Almeida (1987) or Kok (2011a). Figure  . The K14 PSD is shifted toward coarser particle sizes compared to the DA87 PSD, indicating that the DA87 PSD describes dust after more settling of coarse constituents. Both PSDs show a peak in the diameter range 4-8 µm. The mean PSD based on S04 is continuously increasing with particle size, however, the PSD corresponding to the 5th emission-weighted percentile of mean particle diameter with respect to annual emissions does also exhibit a peak around 8 µm, similar to the K14 PSD. In contrast, the S04-PSD belonging to the 95th 400 weighted percentile of mean particle diameter shows a somewhat steeper increase with particle diameter, and correspondingly a smaller fraction of small particles than the median S04 PSD and the DA87 and K14 PSDs. Differences in the PSD of dust at emission yield also differences in airborne dust PSD, which has important effects on the resulting dust optical depth and radiation interactions.  (Hsu et al., 2004;Ginoux et al., 2012, see Sec. 4.3.1). Note that this specification of potential dust source areas, is done in a binary sense and independent of any scaling of saltation or dust emission fluxes (see also next paragraph). This means that dust can be emitted if the topographic mask is non-zero, or the retrieved FoO(DOD>0.2) is 410 greater than a small value (here 0.025) (Fig. 3 top panel). Areas fully covered by vegetation, snow (obtained from reanalysis data used as boundary conditions), or bedrock (from STATSGO-FAO data) are excluded from potential dust sources as described in the final paragraph of Sec. 3.1.4, and a land-sea mask is applied.
In addition to the definition of areas from which dust emission is possible, a scaling of the calculated dust emission fluxes with the above-mentioned dust source functions is deployed in the MB95 and G01 schemes. The preferential source map from improve results compared to observations also for the MB95 implementation (Pérez et al., 2011). We have also added the option to apply the new FoO map as the preferential source function (Fig. 3 center panel). In this case, dust emission is enhanced in 420 areas with high FoO. The purpose of a source map scaling is to compensate for unrepresented processes and surface properties, which affect dust emission. The S01-S11, and K14 schemes are not scaled with any preferential source map. The physics of these schemes are assumed to account for spatial variations in the emitted dust mass and the retrieved FoO map is only used as a mask defining the areas from which dust emission is possible as described above.
An additional special feature of MONARCH is its ability to tag dust originating from natural and anthropogenic (agricul- MONARCH's tagging functionality can be adapted to track dust also from other predefined source origins (Kok et al., 2021b).

430
Vegetation in MONARCH is prescribed based on satellite data, using either an AVHRR monthly climatology of green vegetation cover fraction (Gutman and Ignatov, 2010), or monthly photosynthetic and non-photosynthetic vegetation cover based on MODIS and Landsat surface reflectance (Guerschman et al., 2015) either as a climatology or for the actual year of simulation. The two cover fraction data sets can be used consistently within MONARCH's meteorological and dust modules.
Additionally, monthly MODIS leaf-area-index (LAI) data (Myeni et al., 2015) is available for use in the dust module, for the actual year or as a climatology, by default in combination with the AVHRR climatological vegetation used for meteorology (aerodynamic roughness length and evaporative fluxes). Soil texture class information in MONARCH is obtained from the hybrid STATSGO-FAO data set at a resolution of 30 arc seconds (0.0083°) (Pérez et al., 2011). Additional soil information, such as on soil mineral content, is currently being implemented (Gonçalves Ageitos et al., 2021b, in prep.). To aggregate soil texture data to model resolution, MONARCH 440 utilizes a predominance approach, i.e. the predominant soil texture class in each grid cell is applied to the entire cell.

Dust transport and deposition
Dust transport and deposition in MONARCH has been thoroughly described in Pérez et al. (2011) and is only briefly summarized in this section. The numerical schemes for dust transport by advection and turbulent diffusion are the same as those of other scalars in the NMMB model. Horizontal advection is solved with the Adams-Bashforth scheme and vertical advection 445 with the Crank-Nicholson scheme. Lateral diffusion follows the Smagorinsky non-linear approach. Gravitational settling of dust is solved implicitly from top to bottom using a gravitational settling velocity based on the Stokes-Cunningham approximation. As the settling velocity increasingly deviates from Stokes settling for large particles (approximately >10 µm) and to correct for potential numerical diffusion (Ginoux, 2003) and other unaccounted phenomena (Stout et al., 1995;van der Does et al., 2018;Dey et al., 2019), we successively reduce the settling velocity using bin-wise tuning factors. By default we 450 use 1, 1, 1, 1, 0.5, 0.3, 0.2, 0.1 from bins 1 to 8. An explicit formulation is also now available in the model. Dry deposition through turbulent diffusion is based on Zhang et al. (2001), which accounts for Brownian diffusion, impaction, interception, and gravitational settling (Slinn, 1982). Wet deposition in MONARCH includes in-cloud and below-cloud scavenging from both stratiform (grid-scale) and convective (sub-grid scale) clouds. In-cloud scavenging from stratiform clouds is proportional to dust mass and solubility along with the conversion rate of cloud water to rain by autoconversion, accretion, and shedding 455 of accreted cloud water and to the conversion rate of cloud ice to precipitation through melting. Dust solubility is assumed to have intermediate values between purely hydrophobic and purely hydrophilic particles, with values decreasing with increasing particle size (Zakey et al., 2006). Below-cloud scavenging for rain and snow is based on Slinn (1984) and includes the effects of directional interception, inertial impaction and Brownian diffusion. For convective scavenging, the model follows the principles of the Betts-Miller-Janjic (BMJ) convective parameterization scheme developed by Betts (1986); Betts and Miller 460 (1986);Janjic (1994). In-cloud scavenging is proportional to dust mass and solubility along with the production of precipitation in the convective cloud. Below-cloud scavenging also follows Slinn (1984) assuming a raindrop diameter of 1 mm. BMJ is a convective adjustment scheme and therefore does not represent mass fluxes. Dust is vertically mixed by deep convection in analogy with the vertical adjustment of moisture (Pérez et al., 2011). Currently, dust particles do not affect cloud formation in MONARCH. Parameterizations representing the effect of dust particles on cloud formation, as they act as cloud condensation 465 and ice nuclei, are planned to be implemented in the future.

Radiation and optical properties
The model's radiation scheme is RRTMG (Iacono et al., 2000(Iacono et al., , 2008. MONARCH allows multiple options for setting the dust microphysical properties. In the longwave (LW), we assume refractive indices from the OPAC dataset (Hess et al., 1998) and spherical particle shape. In the shortwave (SW) we use mineralogy-based refractive indices and non-spherical shapes.
470 Table 6. Physical and optical particle properties available in MONARCH for eight particle-size bins: equivalent volume radius (rv), effective radius (re), density (ρp), real and imaginary parts of the refractive index (refREAL, refIMAG), mass-extinction efficiency (MEE, [m 2 g −1 ]), single-scattering albedo (SSA), and asymmetry factor (ASY). The optical properties are for a wavelength of 550 nm and MEE, SSA, and ASY are given assuming ellipsoidal (index ell) or spherical (index sph) particle shape. The diameter ranges of each bin are given in Sec. The multi-component Maxwell Garnett theory (Markel, 2016) (Balkanski et al., 2007) of dust refractive index (Gonçalves Ageitos et al., 2021a, in prep.). We account for the effects of the substantial dust asphericity (Huang et al., 2020) on dust optics by combining the probability distributions of particle shape obtained in Kok et al. (2017) based on laboratory measurements (e.g. Okada et al., 2001;Kandler et al., 2007) 480 with the dust single-scattering database of Meng et al. (2010). Table 6 summarizes key dust properties used in MONARCH.

Model performance and evaluation
A range of global model simulations were performed with MONARCH for one year (2012) to demonstrate MONARCH's dust modeling capabilities. We used different configurations in the runs covering different dust emission schemes. We evaluate the presented simulations against MODIS (Ginoux et al., 2012;Hsu et al., 2013) and Aerosol Robotic Network (AERONET, Giles

Experimental setup
The global model runs performed with MONARCH were conducted at a horizontal resolution of 1°latitude × 1.4°longitude with 48 vertical layers and a computational time step of 3 min. Turbulence, surface layer, dust emission, sedimentation and dry deposition routines were called every 4 computational times steps, moist convection, microphysics and wet scavenging routines 490 every 8 time steps, and short-and longwave radiation routines were called every 20 time steps. The runs were initialized using ERA Interim reanalysis data (Berrisford et al., 2011;Dee et al., 2011). The meteorological fields are re-initialized daily, whereas dust fields and soil moisture are transferred between the daily runs. We used one year of spinup for soil moisture and one month of spinup for the dust fields before the one-year simulation. A simple double-call mechanism computes the total (all size bins together) direct radiative effect (DRE), and a more complex multiple-call mechanism generates the DRE per size bin.

495
The DRE per bin depends on the vertical distribution of particles in a specific bin with respect to those in other bins. Hence, the sum of the DRE per bin does not exactly equal the total DRE, especially for locations with high dust loading. To minimize errors due to such non-linearities, the DRE per bin is calculated as the difference between the total DRE with all bins included (reference state) and the total DRE without the specific bin. Results were output three-hourly for the global runs.
Here we present results of global MONARCH simulations using the MB95, G01-UST, S04, and K14 dust emission schemes, 500 a set of well-known and frequently used parameterizations. In all runs, we scaled soil moisture using c f1 = 0.1 and applied the default soil moisture corrections listed in Tab. 4 (Sec. 3.1.3). We used the drag partition from Marticorena and Bergametti (1995) with X = 12, 255 cm (MacKinnon et al., 2004) for all runs presented here. The intention of using the same drag partition is to ease inter-comparison between the runs, and not to achieve the best possible results for each run. For the latter, different settings for each of the schemes may be more appropriate. The drag partition from Marticorena and Bergametti (1995) uses 505 MODIS LAI as input, while photosynthetic and non-photosynthetic vegetation cover fractions are used for MONARCH's meteorology. Dust emissions in both the MB95 and G01 schemes include a scaling with the topographic source mask from Ginoux et al. (2001) shown in Fig. 3 (bottom), whereas the S04 and K14 schemes do not receive any scaling accounting for preferential dust sources.
The dust fields of all model runs were calibrated using experiment-specific global calibration factors, which were obtained 510 by comparing monthly averages of modeled coarse DOD (size range 1.2-20 µm) for each experiment with the DOD obtained from MODIS (see Sec. 4.3.1 for more detail) and minimizing the overall error . This calibration only removes the general global bias for each run and does not affect the spatio-temporal variability of the dust emission, transport, deposition, and interactions. the S04 scheme produced substantially more dust emission and deposition than the other schemes. This is a result of the on average coarser particle-size distribution in the S04 scheme above 10 µm (Fig. 2) and also reflected in the shorter lifetime of dust aerosol obtained with the S04 experiment (Tab. 7). All experiments were calibrated so that their global DOD resembled that of MODIS. The coarser particles in the S04 experiment have only a small contribution to DOD, but constitute a large amount of the emitted and deposited dust mass.

Dust emission and deposition
pression. In comparison, deposition and dust loading are more intense in northwestern Africa and the Middle East in the S04 scheme, and more homogeneous in the K14 scheme. Figure A1 shows the percent contribution of dust emission and deposition at each location to their respective global and 540 annual totals to investigate differences between the four experiments independent from the overall flux magnitudes. The relative emission (deposition) confirms the differences highlighted before: The spatial patterns of the MB95 and G01-UST are similar due to the use of the preferential source function. In contrast, the S04 experiment produced relatively more dust in north-western Africa, while the K14 scheme generated relatively homogeneous patterns across northern Africa and the Middle East. Kok et al. (2021b) found that dust models often overestimate dust from North African dust sources, but underestimate dust South American dust source about 3-4% in two configurations; Middle Eastern and Central Asian sources about 28-33% in three configurations; East Asian sources > 9% in one configuration; western North African sources contributed 6% more 555 than eastern North African sources in one configuration. Differences across configurations demonstrate the benefit of having multiple dust emission options available in one model and will help to identify aspects of parameterizing dust emission that need particular attention in future research.

Dust optical depth
The global annual average of dust optical depth (DOD) is 0.034, 0.032, 0.041, and 0.035 in the MB95, G01-UST, S04, and 560 K14 runs. This results in an average mass-extinction efficiency MEE g of, respectively, 1.10, 1.11, 1.15, and 1.01 m 2 g −1 for the four runs considering dust up to 20 µm in diameter (Tab. 7), calculated from grid-based annual average DOD and dust load, and, correspondingly, 0.60, 0.57, 0.52, and 0.57 based on global annual average DOD and dust load ( MEE ). The discrepancy between MEE g and MEE is a result of the different emphasis put on locations with high and low dust loading in the two averaging methods: MEE g is calculated from DOD and dust load at each grid cell and then averaged, which puts equal weight 565 on grid cells with and without dust. In contrast, MEE is calculated from globally averaged DOD and dust load and hence focuses more on dusty locations. To provide a comprehensive, yet concise evaluation, we compare the DOD averaged across the four model runs with retrieved DOD from MODIS Deep Blue Sayer et al., 2013) and AERONET (Holben et al., 1998;O'Neill et al., 2003;Giles et al., 2019). Our objective here is to evaluate the overall behavior of MONARCH across dust emission schemes, rather than that of each individual scheme.

Comparison of modelled DOD with MODIS Deep Blue
We estimate DOD from MODIS using daily AOD and SSA at 550 nm, and Ångström exponent ( ison with the MODIS DOD, which discriminates coarse particles from the total AOD, we use modeled DOD in the size-range (1.2-20 µm) and refer to it as DOD coarse . We chose 1.2 µm as the lower diameter limit for DOD coarse as it coincides with that 590 established in the AERONET Inversion product (Dubovik and King, 2000). The size-range of DOD coarse therefore corresponds to the five coarsest bins in MONARCH. Differences in the modeled all-sky co-located and clear-sky DOD coarse underline the impact of the time, location, and number of missing values on the average DOD coarse . The clear-sky DOD coarse tends to be somewhat smaller compared to the all-sky co-located DOD coarse , indicating a discrepancy between modeled and observed clouds, in combination with differences in the underlying DOD. However, the spatial patterns between both model products are overall consistent and achieve comparable pattern correlations and root-mean-square-errors compared to MODIS. The reduced DOD coarse in northern Africa matches even better with the observed DOD, whereas the modeled clear-sky DOD coarse in the Arabian Peninsula is smaller than in the observations. Other areas show very similar results between the all-sky co-located and clear-sky model results.  frames from approximately April until November, and the G01 and S04 runs being intermediate.

Comparison of modelled DOD with AERONET
AERONET is a global network of ground-based solar photometer stations (Holben et al., 1998;O'Neill et al., 2003;Giles 650 et al., 2019). The primary parameter derived by AERONET (i.e. direct-sun) is the AOD in multiple spectral channels with uncertainties lower than 0.03. AOD data are computed for three data quality levels: level 1.0 (unscreened), level 1.5 (cloudscreened), and level 2.0 (cloud-screened and quality-assured). The products from inverting sky radiance measurements are the aerosol size distribution, single scattering albedo, refractive index, effective radius, and asymmetry factor. AERONET has very good coverage across the globe, albeit with lower station density in remote dust source regions, such as northern Africa, the 655 Middle East, central/western Asia, and Australia. Recently, sun-sky-lunar photometers extended the use of AERONET during nighttime (Barreto et al., 2013), allowing continuous aerosol monitoring.
Through AOD, AERONET gives information about the aerosol content and the mode-dominant type (i.e. fine or coarse modes) in the atmospheric column, but not the atmospheric dust burden. Almost pure mineral dust is difficult to find, except in  specific areas close to desert dust sources. Instead, dust is often mixed in variable percentages with other aerosols. To isolate 660 the atmospheric dust burden and estimate the DOD, two approaches are typically used.
The first approach aims to identify records in which the measured aerosol is dominated by mineral dust based on AE. AE is in general inversely related to the average size of the airborne particles and can be used to distinguish species with large particles like dust and sea salt. As a rule of thumb, a larger AE indicates smaller particle size. AE is typically in the range 0-4, where the upper limit corresponds to molecular extinction, and the lower limit corresponds to coarse-mode aerosols (sea-salt 665 and mineral dust), indicating no wavelength dependence of AOD (O'Neill et al., 2003). Since sea-salt is related to low AOD (< 0.03; Dubovik et al. 2002) and mainly affects coastal stations, large coarse-mode AOD values are mainly related to mineral dust. According to previous studies (Dubovik et al., 2002;Wang et al., 2004;Todd et al., 2007;Basart et al., 2009), AE values between 0.75 and 1.2 are associated with mixed aerosols (including dust). An AE lower than 0.2-0.3 is associated with a highly dominant coarse mode in the AERONET bi-modal size distribution (Schuster et al., 2006), which corresponds to almost 670 pure dust conditions over land. Here, we use AE < 0.3 to estimate AERONET DOD for comparison with the DOD (all sizes) obtained from MONARCH.
The second widely used methodology to estimate AERONET DOD is based on the Spectral Deconvolution Algorithm (SDA) retrievals (O'Neill et al., 2003). The SDA algorithm estimates fine (sub-micron) and coarse (super-micron) AOD at a standard wavelength of 500 nm (AOD fine and AOD coarse , respectively). Near dust source regions, DOD coarse ≈ AOD coarse . The advantage of this method is the availability of retrievals in regions where dust occurrence is sporadic and other aerosols are predominant, and where a more restrictive criterion, such as AE < threshold may filter out some dust intrusions (Cuevas et al., 2015). As DOD coarse from MONARCH, we use DOD in the diameter-range 1.2-20 µm.
For comparison with AERONET, we use bilinear interpolation to extract time series from the 3-hourly global model DOD and DOD coarse for the locations of AERONET measurements. We use 3-hourly averages of AERONET observations, such 680 that a comparison with the 3-hourly instantaneous MONARCH data assumes a statistical similarity between the temporally averaged AERONET DOD and DOD coarse and the spatially interpolated MONARCH DOD and DOD coarse . Figure 8 shows Taking into account the entire station list (Appendix C), the correlation is 0.61 with a root-mean-square error (rmse) of 695 0.31 for the total DOD and, 0.71 with an rmse of 0.14 for DOD coarse , based on the three-hourly MONARCH data (Fig. 9).
For total DOD, the correlation remains fairly constant when comparing daily and monthly instead of three-hourly values. The rmse decreases slightly with an increasing averaging period as then discrepancies for individual peaks become less relevant.
A similar behaviour is found for DOD coarse , but with a slightly more pronounced increase also of the correlation. The overall agreement between MONARCH and AERONET is also similar across experiments (Fig. C2), however, the MB95 and S04 700 schemes tend to overestimate events with large DOD, whereas the G01 and K14 show an underestimation of such situations.
As a result, the individual schemes show a slightly lower correlation and higher rmse than the experimental average for both DOD and DOD coarse .  the SFC and TOA for all particle sizes as well as from each bin (relative contribution with respect to the DRE for all sizes) are summarized in Fig. 11 and listed in Tab. D1.

Direct radiative effect
We note that the mineralogy-based set of refractive indices used in this work describes a more scattering dust in the shortwave with respect to other widely used prescriptions (e.g. Patterson et al., 1977;Hess et al., 1998). At the TOA, the shortwave DRE closely oscillates around zero over bright surfaces, such as in the Sahara and Saudi Arabia, and the total DRE does not exceed 720 about 10 Wm −2 in these regions. In contrast, Balkanski et al. (2007) for example obtained a total DRE at the TOA of up to 20 Wm −2 over the Sahara when using Patterson et al. (1977) (Volz (1973) in the longwave), and lower positive values in this region (in agreement with our values) when using the mineralogy-based refractive indices with a hematite content of 1.5% by volume. Note that Balkanski et al. (2007) found these refractive indices to be in a better agreement with AERONET retrievals (Dubovik et al., 2002), similar to what is found by Gonçalves Ageitos et al. (2021a, in prep.) for our refractive indices.

725
Moreover, Miller et al. (2006) obtained a global average total DRE at the TOA of −0.39 Wm −2 using refractive indices from Sinyuk et al. (2003) (Volz (1973) in the longwave). On the other hand, Miller et al. (2014) reported a value of 0.39 Wm −2 calculated using the dust distribution from Miller et al. (2006) and refractive indices from Patterson et al. (1977). Our negative value of −0.24 Wm −2 is therefore again more comparable with more scattering dust as described by the refractive indices of Sinyuk et al. (2003). The inclusion of particles larger than 20 µm would likely lead to a slightly more positive DRE.  (Rubin et al., 2016). This can be achieved using different meteorological fields as initial and boundary conditions in the meteorological driver of MONARCH (NMMB) for each forecast run in the ensemble, in addition to the dust perturbations. In the dust reanalysis currently in production at the BSC, we use 755 an ensemble based on stochastic perturbations of emission parameters, in conjunction with emission schemes with different physics and different meteorological initial and boundary conditions (Di Tomaso et al., 2021, in prep.).
Airborne dust is not a homogeneous entity, but a mixture of minerals, the relative amounts of which depend on the source region. Mineralogy affects a variety of dust-related impacts, e.g. interaction with radiation, atmospheric chemistry or nutrient supply to certain ecosystems. The capability to explicitly represent dust composition was recently added to MONARCH al-760 lowing the tagging of up to 12 different minerals. This new feature is currently used to assess the relevance of dust mineralogy for dust impacts and to provide insights for the near-term atmospheric and climate modeling communities (Gonçalves Ageitos et al., 2021b, in prep.).
The combination of different vegetation input data sets, drag partition approaches, and the source tagging capability, allows to represent the seasonal vegetation dynamics and provides an ideal basis to investigate the importance of dust from anthropogenic 765 (agricultural) sources, for which a key driver is the seasonal vegetation growth and decay. The benefit of online estimates within a modeling framework is that not only the emission, but also the transport, deposition, and effect of anthropogenic dust can be investigated (Klose et al., 2018).
Code availability. The MONARCH Version 2.0.0 source code used in this work is available at https://earth.bsc.es/gitlab/es/monarch.
Data availability. MONARCH output presented in this paper will be available via a EUDAT B2SHARE
Appendix A: Dust emission and deposition Figure A1 shows the percent contribution of dust emission and deposition at each location to their respective global and annual totals to visualize regional differences between the different experiments independent from the overall emission (deposition) magnitudes obtained.

Appendix C: Comparison with AERONET
The AERONET stations used for comparison with MONARCH and to obtain the global calibration factor are listed in Tab. C1 and shown in Fig. C1. They cover the main dust source regions around the globe. The intention of using only a subset of all stations is to increase confidence in that aerosol detected by AERONET photometers is predominantly dust. Appendix D: Size-dependent dust direct radiative effect Table D1 gives the relative contributions of each particle-size bin to the global average DRE for each MONARCH run.         List of symbols α coefficient in soil moisture correction from Klose et al. (2014) α q vertical-to-horizontal-flux ratio β coefficient in soil moisture correction from Klose et al. (2014) β R coefficient in drag partition correction from Raupach et al. (1993)  σ v ratio of roughness-element basal to frontal area σ pi η mi /η f i τ total stress τ s maximum surface stress on exposed area τ s average surface stress on exposed area a coefficient in soil moisture correction from Fécan et al. (1999) b coefficient in soil moisture correction from Fécan et al. (1999)  X parameter in drag partition from Marticorena and Bergametti (1995) 905 z 0 aerodynamic roughness length z 0dyn dynamic aerodynamic roughness length z 0stat static aerodynamic roughness length z 0s smooth aerodynamic roughness length Boucher, O., Randall, D., Artaxo, P., Bretherton, C., Feingold, G., Forster, P., Kerminen, V.-M., Kondo, Y., Liao, H., Lohmann, U., Rasch, P., Satheesh, S. K., Sherwood, S., Stevens, B., and Zhang, X. Y.: Clouds and Aerosols, in: Climate Change 2013: The Physical Science Basis.