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 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, single-scattering albedo, and direct radiative effects. For dust-particle diameters up to 20 µm, the total annual dust emission and deposition fluxes obtained with our four experiments range between about 3500 and 6000 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⁻² and at the top of the atmosphere between −0.20 and −0.28 W m⁻². Our evaluation demonstrates that MONARCH is able to reproduce key features of the spatiotemporal variability of the global dust cycle with important and insightful differences between the different configurations.

1 Introduction

The Multiscale Online Non-hydrostatic AtmospheRe CHemistry model (MONARCH) is a chemical weather modeling system that can be used at multiple spatial scales, ranging from regional scales at single-digit kilometer resolutions with explicit convection to coarse-resolution global scales with parameterized convection (Pérez et al., 2011; Badia et al., 2017). MONARCH is continuously developed at the Barcelona Supercomputing Center (BSC) with a focus on mineral dust and other aerosols (Pérez et al., 2011; Haustein et al., 2012; Spada et al., 2013; Spada, 2015), atmospheric chemistry (Jorba et al., 2012; Badia and Jorba, 2015; Badia et al., 2017), emissions (HERMES, Guevara et al., 2019), data assimilation (Di Tomaso et al., 2017), workflow management (Manubens-Gil et al., 2016), evaluation (Binietoglou et al., 2015; Ansmann et al., 2017), and operational forecasting (Basart et al., 2019; Xian et al., 2019). Daily dust forecasts using MONARCH are produced at the BSC and made available through the Barcelona Dust Forecast Center (a WMO regional specialized meteorological center with activity specialization on atmospheric sand and dust forecasts; https://dust.aemet.es/, last access: 8 October 2021), the WMO SDS-WAS Regional Center for Northern Africa–Middle East–Europe (NA-ME-E) (https://sds-was.aemet.es/, last access: 8 October 2021), and the International Cooperative for Aerosol Prediction (ICAP) (Sessions et al., 2015; Xian et al., 2019). Here, we present recent developments on the representation of mineral dust processes in MONARCH.

Mineral soil dust is the most abundant aerosol type in terms of global mass, competing only with sea salt (Textor et al., 2006). Global dust emissions are estimated to range between 3300 and 9000 Tg yr⁻¹ for particles smaller than 20 µm in diameter (Kok et al., 2021a). Soil dust is mainly emitted from arid and semi-arid regions, e.g., in Africa, the Middle East, central and northeastern Asia, India, Australia, Patagonia, and the southwestern United States, but it can in principle be emitted from any uncovered dry soil surface under windy conditions, e.g., from agricultural fields (Ginoux et al., 2012).

Mineral dust is emitted as soon as the forces that act to retain the soil particles at the surface (gravity and interparticle cohesion) are overcome either by atmospheric lifting forces generated by wind and turbulence (aerodynamic entrainment) or by the force generated by other impacting particles, i.e., sand grains or particle aggregates (saltation bombardment/aggregate disintegration) (Shao, 2008). Typically, soil particles in the silt and clay particle-size range (diameter <63 µm; Udden, 1914; Wentworth, 1922) are considered dust, whereas larger particles are referred to as sand. Soil particles in the size range 70–100 µm can typically be lifted most easily (Iversen and White, 1982; Shao and Lu, 2000). For larger diameters, the particle weight is the predominant inhibitor. For smaller particles, interparticle cohesion becomes more significant but likely exhibits stochastic behavior leaving a fraction of particles with substantially below-average cohesive forces (Shao and Klose, 2016). Particles with diameters of around 70 µm and larger are mainly transported in ballistic trajectories along the surface (saltation). The limit at which saltation is initiated, i.e., when the particle retarding forces are exceeded by the aerodynamic lifting forces, is expressed as a threshold friction velocity, u_*t [ms⁻¹]. Saltation bombardment is typically most efficient at generating dust emission (e.g., Shao et al., 1993; Houser and Nickling, 2001). Aerodynamic dust entrainment is typically less efficient because of on-average higher cohesive forces for dust-sized particles compared to sand (or saltation) particles but can be significant under favorable atmospheric conditions and on long timescales, provided there is a sufficient supply of loose dust particles at the surface (e.g., Loosmore and Hunt, 2000; Macpherson et al., 2008; Chkhetiani et al., 2012; Klose and Shao, 2013; Li et al., 2014; Zhang et al., 2016).

Once airborne, mineral dust particles interact with short- and longwave radiation through scattering and absorption (Boucher et al., 2013; Miller et al., 2014), which has important direct effects on the Earth's energy balance (Kok et al., 2018; Li et al., 2020). Dust particles are known to be efficient ice nuclei and can also act as cloud condensation nuclei (DeMott et al., 2003; Karydis et al., 2011; Cziczo et al., 2013; Kiselev et al., 2017). Nutrients transported with dust can create ecosystem responses due to, e.g., carbon uptake and storage (Jickells et al., 2005; Rizzolo et al., 2017; Kanakidou et al., 2018). Dust can cause respiratory and cardiovascular diseases (Meng and Lu, 2007), and can contribute to other ailments like meningitis (Pérez García-Pando et al., 2014) and valley fever (Tong et al., 2017). To study and quantify dust and its impacts, models that include advanced dust representations, such as MONARCH, are key tools.

Ferrier et al. (2002)Thompson et al. (2008)Hong and Lim (2006)Iacono et al. (2008)Fels and Schwarzkopf (1975)Janjic (1994, 1996b)Ek et al. (2003)Vukovic et al. (2010)Janjic (1996a, 2001)Hong and Pan (1996)Betts (1986); Betts and Miller (1986); Janjic (1994, 2000)Han and Pan (2011)

Table 1 Available physics schemes in MONARCH.

Existing dust emission parameterizations range from formulations that are strongly simplified (e.g., Ginoux et al., 2001) to those that aim to represent the physics of the emission processes (Shao, 2004; Klose et al., 2014; Kok et al., 2014b). The more simplified dust emission schemes are typically constrained by “preferential” source scaling functions and are commonly used in global but also in regional models. Such constraints have significantly improved the skill of models by approximately locating and enhancing dust emissions from prolific large-scale natural sources. However, these schemes are not very sensitive to changes in parameters known to affect dust emission (e.g., soil texture, soil moisture, land-surface properties), which at the same time can make models insensitive to changes in climate. In contrast, physics-based dust emission parameterizations are very sensitive to such changes but need more detailed input. This detailed input has traditionally been difficult to observe and/or estimate, in particular globally, and errors in the description of, for example, surface properties, translate non-linearly into errors in emitted and transported dust. How such errors compare with those arising when neglecting dust emission sensitivities entirely remains a subject of research and discussion. Modeling the dust cycle for past (e.g., mid-Holocene, Last Glacial Maximum, rapid transitions) or future periods requires balancing a sensitive parameterization with input data that are more uncertain than for the present day. All in all, a clear benefit of physics-based schemes with detailed sensitivities is that input data sets can easily be updated as more data become available and hence future improvements are more likely, in particular for climate applications.

In this work, we introduce recent advancements in the treatment of mineral dust in MONARCH. The model now has diverse available model configurations, in particular to estimate dust emission, which makes MONARCH unique among state-of-the-art models, and which makes it suitable for a variety of applications that range from process studies to operational forecasting and climate research. In the following sections, we briefly present the MONARCH modeling system and subsequently focus on the mineral dust cycle. We then demonstrate and evaluate MONARCH's dust modeling capabilities based on four annual global model runs.

Marticorena and Bergametti (1995)Ginoux et al. (2001, modified)Shao (2001)Shao (2004)Shao et al. (2011)Kok et al. (2014b)

Table 2 Summary of the six available dust emission schemes in MONARCH.

2 The MONARCH model

MONARCH (previously known as NMMB/BSC-CTM) consists of advanced chemistry and aerosol packages coupled online with the Non-hydrostatic Multiscale Model on the B-grid (NMMB) (Janjic et al., 2001; Janjic and Gall, 2012), whose nonhydrostatic dynamical core allows running both global and regional simulations with embedded telescoping nests. The global model works on a latitude–longitude grid with polar filtering and the regional model on a rotated longitude–latitude grid. In both cases, the Arakawa B grid and the hybrid pressure–sigma coordinate are used in the horizontal and vertical directions, respectively. The numerical schemes follow the principles described in Janjic (1977, 1979, 1984, 2003). The NMMB can be configured with a combination of different physics schemes (see Table 1). The configuration commonly used in production runs and in this work is as follows. Turbulence in the planetary boundary layer and the free troposphere is resolved using the Mellor–Yamada–Janjic (MYJ) level 2.5 turbulence closure scheme (Janjic, 2001). The surface layer scheme combines Monin–Obukhov similarity theory (Monin and Obukhov, 1954) with a viscous sublayer introduced over land and water (Zilitinkevich, 1965; Janjic, 1994, 1996b). The shortwave and longwave radiation fluxes are computed using the Rapid Radiative Transfer Model for general circulation models (GCMs) (RRTM-G) radiation package (Iacono et al., 2008). The model includes the Ferrier scheme for grid-scale cloud microphysics (Ferrier et al., 2002), and the Betts–Miller–Janjic convective adjustment scheme (Betts, 1986; Betts and Miller, 1986; Janjic, 1994, 2000). The Unified NCEP/NCAR/AFWA Noah (Ek et al., 2003) land surface model is used for the computation of heat and moisture surface fluxes.

The gas-phase chemistry in MONARCH solves the Carbon Bond 2005 chemical mechanism (CB05; Yarwood et al., 2005) extended with toluene and chlorine chemistry. The CB05 is well formulated for urban to remote tropospheric conditions, and it considers 51 chemical species and solves 156 reactions. The photolysis rates are computed with the Fast-J scheme (Wild et al., 2000) considering the physics of each model layer (e.g., aerosols, clouds, absorbers such as ozone). The aerosol module in MONARCH describes the life cycle of dust, sea salt, black carbon, organic matter (both primary and secondary), sulfate, and nitrate aerosols (Spada, 2015). While a sectional approach is used for dust and sea salt, a bulk description of the other aerosol species is currently adopted (Spada, 2015). A simplified gas–aqueous–aerosol mechanism has been introduced in the module to account for the sulfur chemistry, the production of secondary nitrate–ammonium aerosol is solved using the thermodynamic equilibrium model EQSAM (Metzger et al., 2002), and a two-product scheme is used for the formation of secondary organic aerosols from gas-phase precursors. Different meteorology-driven emissions are computed online in MONARCH (i.e., mineral dust, sea salt, and biogenic gas species). Sea salt emissions can be calculated with a wide range of available source functions (Spada et al., 2013), while the biogenic emissions are estimated with the Model of Emissions of Gases and Aerosols from Nature (MEGAN) v2.04 model (Guenther et al., 2006).

In addition to single model runs, MONARCH can be run in an ensemble mode for data assimilation applications, where the ensemble of model states is used to derive a flow-dependent background error covariance at the assimilation time, which evolves during the model forecast. The background error covariance is used to express model uncertainty within the data assimilation framework. Model uncertainty, together with observational uncertainty, is a key ingredient in the optimal integration of model simulations and observations for the production of an analysis that best represents the atmospheric state. The MONARCH ensemble is coupled with the local ensemble transform Kalman filter (LETKF) scheme (Hunt et al., 2007; Miyoshi and Yamane, 2007; Schutgens et al., 2010) for the estimation of dust analyses (Di Tomaso et al., 2017), as well as for the generation of dust reanalyses currently in production at the BSC (Di Tomaso et al., 2021).

3 The mineral dust cycle in MONARCH

The dust module in MONARCH (previously known as NMMB/BSC-Dust), initially described by Pérez et al. (2011), solves the mass balance equation for dust taking into account the following processes:

1. dust generation and uplift by surface wind and turbulence,

2. horizontal and vertical advection,

3. horizontal diffusion and vertical transport by turbulence and convection,

4. dry deposition and gravitational settling, and

5. wet removal by convective and stratiform clouds.

The dust size distribution is represented with eight bins ranging up to 20 µm in diameter: 0.2–0.36, 0.36–0.6, 0.6–1.2, 1.2–2.0, 2.0–3.6, 3.6–6.0, 6.0–12.0, and 12.0–20.0 µm. The effective and volume radii of each bin in the radiative and sedimentation schemes, respectively (see Sect. 3.3, Table 6) are time invariant and based on a lognormal distribution with mass median diameter of 2.524 µm and 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 data sets 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 diagnostic capabilities (output of three-dimensional single-scattering albedo and extinction, clear-sky aerosol optical depth (clear-sky AOD), and AOD at satellite overpass times). In the following, we present the MONARCH dust module. We first describe the treatment of dust emission, summarizing previous and detailing new developments. Then, we recapitulate the implementation of dust transport and deposition, and interactions with radiation.

3.1 Dust emission and lower boundary conditions

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 (2001, 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 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 Table 2.

Table 3 Coefficients for minimally dispersed particle-size distributions as assigned to the 12 USGS soil texture classes. Each particle-size distribution (PSD) is composed of four lognormal distributions (p₁, p₂, p₃, and p₄). Coefficients are taken from Klose (2014), Table 3 (and references therein), unless otherwise indicated.

¹ Coefficients for samples from Site D in Klose et al. (2019) (PSD_LEM).
² Coefficients from Table 3 of Klose et al. (2019) (PSD_LEM).
³ Different sample used as reference than in Table 6.1 of Klose (2014) but same underlying data set.
⁴ Sandy clay loam in Table 6.1 of Klose (2014).

Belly (1964)Ginoux et al. (2001)Fécan et al. (1999)Zender et al. (2003)Tegen et al. (2002)Fécan et al. (1999)Zender et al. (2003)Kok et al. (2014b)Klose et al. (2014)

Table 4 Summary of available options in MONARCH to account for soil moisture in the particle entrainment threshold.

Marticorena and Bergametti (1995); King et al. (2005)(Prigent et al., 2012)(Myeni et al., 2015)Raupach et al. (1993)(Guerschman et al., 2015, photosynthetic and non-photosynthetic vegetation)(Gutman and Ignatov, 2010, green vegetation)

Table 5 Summary of available options in MONARCH to account for surface roughness in particle entrainment.

3.1.1 Dust emission flux

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,

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }{F}_{\text{MB95}}=S{\mathit{\alpha }}_{q}Q.\end{array}$$

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 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 Sect. 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 (Sect. 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

$$\begin{array}{}\text{(2)}& F_{\text{G01-U}}\left(d_i\right)=\left\{\begin{array}{ll}{s}_{\text{bare}}{C}_{\text{G01}}S{s}_{\mathrm{p}}{u}_{\text{10m}}^{\mathrm{2}}\left(u_{\text{10m}}-u_{\mathrm{t}}\right)& u_{\text{10m}}>u_{\mathrm{t}}\\ \mathrm{0}& \text{otherwise}\end{array}\right.,\end{array}$$

where C_G01 is a dimensional factor, s_bare is the bare soil fraction (see Sect. 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 (Sect. 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 (Sect. 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), instead of u_10m and u_t (G01-U), such that

$$\begin{array}{}\text{(3)}& F_{\text{G01-UST}}\left(d_i\right)=\left\{\begin{array}{ll}{s}_{\text{bare}}{C}_{\text{G01}}S{s}_{\mathrm{p}}{u}_{*}^{\mathrm{2}}\left(u_{*}-u_{*\mathrm{t}}\right)& u_{*}>u_{*\mathrm{t}}\\ \mathrm{0}& \text{otherwise}\end{array}\right..\end{array}$$

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

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{\stackrel{\mathrm{̃}}{F}}_{\text{S01}}\left(d_i,d_{\mathrm{s}}\right)& =s_{\text{bare}}{c}_y\left[\left(\mathrm{1}-\mathit{\gamma }\right)+\mathit{\gamma }{\mathit{\sigma }}_{{\mathrm{p}}_i}\right]\times {\displaystyle \frac{g{Q}_{\mathrm{s}}}{u_{*}^{\mathrm{2}}{m}_{{\mathrm{p}}_{\mathrm{s}}}}}\\ & \times \left({\mathit{\rho }}_{\mathrm{b}}{\mathit{\eta }}_{\text{fi}}\mathrm{\Omega }+{\mathit{\eta }}_{\text{ci}}{m}_{{\mathrm{p}}_{\mathrm{s}}}\right),\end{array}\end{array}$$

with

$$\begin{array}{}\text{(5)}& \mathit{\gamma }=\exp \left[-\mathit{\kappa }{\left(u_{*}-u_{*\mathrm{t}}\right)}^{\mathrm{3}}\right],\end{array}$$

where​​​​​​​ c_y and κ are coefficients, η_fi, η_ci, and η_mi are the total, aggregated, and free dust fractions at diameter d_i, ${\mathit{\sigma }}_{{\mathrm{p}}_i}={\mathit{\eta }}_{\text{mi}}/{\mathit{\eta }}_{\text{fi}}$, Q_s the saltation flux of particles with diameter d_s, $m_{{\mathrm{p}}_{\mathrm{s}}}=m_{\mathrm{p}}\left(d_{\mathrm{s}}\right)$ the mass of a spherical particle with diameter d_s assuming a density of ρ_ps=2650 kg m⁻³, ρ_b≈1000 kg m⁻³ the soil bulk density, g gravitational acceleration, and Ω the soil volume removed by a saltating particle of size d_s (Lu and Shao, 1999, Eq. 8). The removed soil mass is given by m_Ω=ρ_bΩ. The bare soil fraction s_bare was added here for implementation in MONARCH. The dust fractions η_fi, η_ci, and η_mi can be estimated from the minimally and fully dispersed particle-size distributions (PSDs), p_m(d_i) and p_f(d_i), as ${\mathit{\eta }}_{\text{mi}}={\int }_{d_{i\mathrm{0}}}^{d_{i\mathrm{1}}}{p}_{\mathrm{m}}\left(d_i\right)\mathit{\delta }{d}_i$ and ${\mathit{\eta }}_{\text{fi}}={\int }_{d_{i\mathrm{0}}}^{d_{i\mathrm{1}}}{p}_{\mathrm{f}}\left(d_i\right)\mathit{\delta }{d}_i$, where d_i0 and d_i1 are the lower and upper limits of the particle-size bin corresponding to d_i. The aggregated dust fraction follows as ${\mathit{\eta }}_{\text{ci}}={\mathit{\eta }}_{\text{fi}}-{\mathit{\eta }}_{\text{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 Sect. 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 emission flux of dust particles with diameter d_i, i.e., for all saltation particle sizes, is obtained as

$$\begin{array}{}\text{(6)}& {\displaystyle }{\displaystyle }F\left(d_i\right)=\sum _{ds=d_{\text{max}}}^{\mathrm{\infty }}\stackrel{\mathrm{̃}}{F}\left(d_i,d_{\mathrm{s}}\right),\end{array}$$

with d_max=20 µm in MONARCH.

The S04 scheme is a simplification of the S01 scheme in which the saltation bombardment efficiency, ${\mathit{\sigma }}_m=m_{\mathrm{\Omega }}/m_{{\mathrm{p}}_{\mathrm{s}}}$, is approximated as (Shao, 2004, Eq. 11)

$$\begin{array}{}\text{(7)}& {\displaystyle }{\displaystyle }{\mathit{\sigma }}_m=\mathrm{12}{u}_{*}^{\mathrm{2}}{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{b}}}{P}}\left(\mathrm{1}+\mathrm{14}{u}_{*}\sqrt{{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{b}}}{P}}}\right),\end{array}$$

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)

$$\begin{array}{}\text{(8)}& {\stackrel{\mathrm{̃}}{F}}_{\text{S04}}\left(d_i,d_{\mathrm{s}}\right)=s_{\text{bare}}{c}_y{\mathit{\eta }}_{\text{fi}}\left[\left(\mathrm{1}-\mathit{\gamma }\right)+\mathit{\gamma }{\mathit{\sigma }}_{{\mathrm{p}}_i}\right]\left(\mathrm{1}+{\mathit{\sigma }}_m\right){\displaystyle \frac{g{Q}_{\mathrm{s}}}{u_{*}^{\mathrm{2}}}},\end{array}$$

and F(d_i) follows from Eq. (4) assuming η_fi≈η_ci. Based on a detailed comparison with field measurements, a basic version 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:

$$\begin{array}{}\text{(9)}& {\displaystyle }{\displaystyle }{F}_{\text{S11}}\left(d_i\right)=s_{\text{bare}}{c}_y{\mathit{\eta }}_{\text{mi}}\left(\mathrm{1}+{\mathit{\sigma }}_m\right){\displaystyle \frac{gQ}{u_{*}^{\mathrm{2}}}}.\end{array}$$

Shao et al. (2011) 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). 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

$$\begin{array}{}\text{(10)}& {\displaystyle }\begin{array}{rl}{F}_{\text{K14}}\left(d_i\right)& =C_{\mathrm{e}}{s}_{\text{bare}}{f}_{\text{clay}}{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{a}}\left(u_{*}^{\mathrm{2}}-u_{*\mathrm{t}}^{\mathrm{2}}\right)}{u_{*\text{st}}}}{\left({\displaystyle \frac{u_{*}}{u_{*\mathrm{t}}}}\right)}^{C_{\mathit{\alpha }}{\mathit{\psi }}_{*}}\\ & \text{for}{u}_{*}>u_{*\mathrm{t}}\end{array}\text{(11)}& {\displaystyle }{\mathit{\psi }}_{*}={\displaystyle \frac{u_{*\text{st}}-u_{*\text{st0}}}{u_{*\text{st0}}}},\end{array}$$

where f_clay is the clay fraction (from STATSGO-FAO inventory; see Sect. 3.1.6), ρ_a is air density, $u_{*\text{st}}=u_{*\mathrm{t}}\sqrt{{\mathit{\rho }}_{\mathrm{a}}/{\mathit{\rho }}_{\text{a0}}}$ with ρ_a0=1.225 kg m⁻³ is a standardized threshold friction velocity, C_α is a constant coefficient, and C_e is a u_*st-dependent coefficient representing soil erodibility.

3.1.2 Saltation flux

For the schemes that contain explicit representations of the saltation flux (MB95, S01, S04, S11), the saltation flux of particles with diameter d_s, Q_s, is calculated following Kawamura (1964) (same as White, 1979) as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{Q}_{\mathrm{s}}\left(d_{\mathrm{s}}\right)=c_Q{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{a}}}{g}}{u}_{*}^{\mathrm{3}}\left(\mathrm{1}+{\displaystyle \frac{u_{*\mathrm{t}}\left(d_{\mathrm{s}}\right)}{u_{*}}}\right)\left(\mathrm{1}-{\displaystyle \frac{u_{*\mathrm{t}}{\left(d_{\mathrm{s}}\right)}^{\mathrm{2}}}{u_{*}^{\mathrm{2}}}}\right)\\ \text{(12)}& {\displaystyle }& {\displaystyle }\text{for }{u}_{*}>u_{*\mathrm{t}}\left(d_{\mathrm{s}}\right),\end{array}$$

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

$$\begin{array}{}\text{(13)}& {\displaystyle }& {\displaystyle }Q=\int Q_{\mathrm{s}}\left(d\right)p_{\mathrm{s}}\left(d\right)\mathit{\delta }d,\text{ with}\text{(14)}& {\displaystyle }& {\displaystyle }{p}_{\mathrm{s}}\left(d_{\mathrm{s}}\right)=\mathit{\gamma }{p}_{\mathrm{m}}\left(d_{\mathrm{s}}\right)+\left(\mathrm{1}-\mathit{\gamma }\right)p_{\mathrm{f}}\left(d_{\mathrm{s}}\right).\end{array}$$

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 Table 3. PSDs are calculated with 60 size bins distributed logarithmically using a quarter-φ scale (Krumbein, 1934, 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. Sect. 3.1.1). The G01 and K14 schemes do not contain explicit formulations for saltation flux.

3.1.3 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 scheme and its requirements. In the MB95 implementation, we use the relationship from Iversen and White (1982) for the four saltation size classes (cf. Sect. 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 (Sect. 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 Sect. 1). In combination with the relatively simple and constant distribution of soil particles across clay and silt particle sizes (Ginoux et al., 2001), this 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

$$\begin{array}{}\text{(15)}& {\displaystyle }{\displaystyle }{u}_{*{\mathrm{t}}_{\text{d0}}}=\underset{d_i}{min}\left[u_{*\text{t0}}\left(d_i\right)\right],\end{array}$$

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_{{\mathrm{t}}_{\text{d0}}}=\mathrm{5}$ m s⁻¹, a typical 10 m wind lower limit for dust emission under favorable land-surface conditions (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 Sect. 3.1.5). The K14 scheme also makes use of a particle-size-independent threshold friction velocity and in MONARCH, $u_{*{\mathrm{t}}_{\text{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 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 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 subgrid-scale wind variability (e.g., Cakmur et al., 2004; Lunt and Valdes, 2002). In our model, we included an optional constant scaling parameter, c_thr≤1, such that the final threshold friction velocity for dry conditions, $u_{*{\mathrm{t}}_{\text{dry}}}$ is

$$\begin{array}{}\text{(16)}& {\displaystyle }{\displaystyle }{u}_{*{\mathrm{t}}_{\text{dry}}}\left(d_i\right)=c_{\text{thr}}\cdot u_{*\text{t0}}\left(d_i\right).\end{array}$$

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_{*{\mathrm{t}}_{\text{dry}}}$, and then applying a correction factor, f_w>1, to obtain the threshold friction velocity for the given (moist) conditions (McKenna Neuman and Nickling, 1989; Fécan et al., 1999; McKenna Neuman, 2003; Cornelis et al., 2004a, b; Klose et al., 2014):

$$\begin{array}{}\text{(17)}& {\displaystyle }{\displaystyle }{u}_{*\mathrm{t}}\left(d_i\right)=u_{*{\mathrm{t}}_{\text{dry}}}\left(d_i\right)\cdot f_w.\end{array}$$

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 Table 4 and further detailed below.

The soil moisture correction from Belly (1964) is implemented as described in Ginoux et al. (2001):

$$\begin{array}{}\text{(18)}& f_{w_{\mathrm{B}}}=\left\{\begin{array}{ll}\mathrm{1.2}+\mathrm{0.2}{\log }_{\mathrm{10}}\left(max\left(\mathrm{0.001},c_{{\mathrm{f}}_{\mathrm{1}}}\cdot \mathit{\theta }\right)\right)& \mathit{\theta }<\mathrm{0.5}\\ f_{w_{\text{B\hspace{0.17em}wet}}}& \text{otherwise}\end{array}\right.,\end{array}$$

where θ is the volumetric soil moisture [m³ m⁻³], $f_{w_{\text{B\hspace{0.17em}wet}}}$ is a large value prohibiting particle movement (here $f_{w_{\text{B\hspace{0.17em}wet}}}=\mathrm{100}$), and $c_{{\mathrm{f}}_{\mathrm{1}}}$ 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

$$\begin{array}{}\text{(19)}& {\displaystyle }{\displaystyle }{f}_{w_{F_w}}=\sqrt{\mathrm{1}+a{\left(c_{{\mathrm{f}}_{\mathrm{1}}}\cdot w-c_{{\mathrm{f}}_{\mathrm{2}}}\cdot w_{\mathrm{r}}\right)}^b}{c}_{{\mathrm{f}}_{\mathrm{1}}}w>c_{{\mathrm{f}}_{\mathrm{2}}}{w}_{\mathrm{r}},\end{array}$$

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

$$\begin{array}{}\text{(20)}& {\displaystyle }& {\displaystyle }w=\mathrm{100}\cdot \mathit{\theta }{\displaystyle \frac{{\mathit{\rho }}_l}{{\mathit{\rho }}_{\text{bd}}}},\text{(21)}& {\displaystyle }& {\displaystyle }{\mathit{\rho }}_{\text{bd}}={\mathit{\rho }}_{\text{pa}}\left(\mathrm{1}-{\mathit{\theta }}_{s_z}\right),\text{ and}\text{(22)}& {\displaystyle }& {\displaystyle }{\mathit{\theta }}_{s_z}=\mathrm{0.489}-\mathrm{0.126}{M}_{\text{sand}}.\end{array}$$

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⁻³), ${\mathit{\theta }}_{s_z}$ is the volumetric soil moisture at saturation, and M_sand is the sand fraction in the soil (Pérez et al., 2011, Table 1). The factor of 100 converts soil moisture content from kg kg⁻¹ into %. As the top-layer soil moisture in models is usually 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_{{\mathrm{f}}_{\mathrm{1}}}$ or $c_{{\mathrm{f}}_{\mathrm{2}}}$, can therefore be applied if needed. The coefficient $c_{{\mathrm{f}}_{\mathrm{1}}}\le \mathrm{1}$ directly reduces the soil moisture in Eq. (19) (e.g., Shao et al., 2010). Soil moisture remains unmodified outside of Eq. (19). Alternatively, the coefficient $c_{{\mathrm{f}}_{\mathrm{2}}}\ge \mathrm{1}$ (Zender et al., 2003) instead increases the air-dry soil moisture. Both coefficients have the effect to reduce f_w. We recommend using either $c_{{\mathrm{f}}_{\mathrm{1}}}$ or $c_{{\mathrm{f}}_{\mathrm{2}}}$, 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_{{\mathrm{f}}_{\mathrm{1}}}$, the correction is

$$\begin{array}{}\text{(23)}& {\displaystyle }{\displaystyle }{f}_{w_{\mathrm{K}}}=\sqrt{\mathrm{1}+{\displaystyle \frac{h_w}{{\mathit{\psi }}_{\mathrm{s}}}}{\left({\displaystyle \frac{c_{{\mathrm{f}}_{\mathrm{1}}}\mathit{\theta }-{\mathit{\theta }}_{\mathrm{r}}}{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}}}\right)}^{\mathit{\beta }}}{c}_{{\mathrm{f}}_{\mathrm{1}}}\mathit{\theta }>{\mathit{\theta }}_{\mathrm{r}},\end{array}$$

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 (Klose et al., 2014). Equation (23) is consistent with Eq. (19), as can be seen when setting

$$\begin{array}{}\text{(24)}& {\displaystyle }{\displaystyle }\mathit{\alpha }={\displaystyle \frac{h_w}{{\mathit{\psi }}_{\mathrm{s}}}}{\left({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}\right)}^{-\mathit{\beta }}.\end{array}$$

Note that the volumetric soil moisture θ [m³ m⁻³] 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 et al. (2014, Table 1). The optional tuning constant $c_{{\mathrm{f}}_{\mathrm{2}}}$ was not implemented in Eq. (23) for simplicity, as this would require modifying α.

3.1.4 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. (1993) or that of Marticorena and Bergametti (1995) with a correction published in King et al. (2005) (Table 5). Typically the 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

$$\begin{array}{}\text{(25)}& {\displaystyle }{\displaystyle }{u}_{*}=f_v\cdot u_{*\text{NMMB}}.\end{array}$$

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

$$\begin{array}{}\text{(26)}& {\displaystyle }{\displaystyle }{f}_{v_{\mathrm{R}}}={\left({\displaystyle \frac{{\mathit{\tau }}_{\mathrm{s}}^{\prime \prime }}{\mathit{\tau }}}\right)}^{\mathrm{1}/\mathrm{2}}={\left[{\displaystyle \frac{\mathrm{1}}{\left(\mathrm{1}-m{\mathit{\sigma }}_v\mathit{\lambda }\right)\left(\mathrm{1}+m{\mathit{\beta }}_{\mathrm{R}}\mathit{\lambda }\right)}}\right]}^{\mathrm{1}/\mathrm{2}},\end{array}$$

where τ is the total stress, ${\mathit{\tau }}_{\mathrm{s}}^{\prime \prime }={\mathit{\tau }}_{\mathrm{s}}^{\prime }\left(m\mathit{\lambda }\right)$ is the maximum surface stress on the exposed area estimated from the average surface shear stress on the exposed area, ${\mathit{\tau }}_{\mathrm{s}}^{\prime }$, for a surface with lower roughness density using the constant m≤1, σ_v is the ratio of the 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)

$$\begin{array}{}\text{(27)}& {\displaystyle }{\displaystyle }\mathit{\lambda }=-c_{\mathit{\lambda }}\ln \left(\mathrm{1}-\mathit{\eta }\right),\end{array}$$

where c_λ is a coefficient. If the roughness elements are uniformly distributed and isotropically oriented, c_λ=1 (Raupach et al., 1993; 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 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, 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–2017) and also of green vegetation cover estimated from the Advanced Very High Resolution Radiometer (AVHRR) (1985–1990, Gutman and Ignatov, 2010) are also available.

Figure 1a–c show annual averages for 2012 of PV, NPV, and $f_{v_{\mathrm{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, and 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 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_{\mathrm{R}}}$ will lead to changes in the roughness correction, particularly in areas with moderate vegetation coverage.

Figure 1 Comparison of roughness input and drag partition approaches: panels (a, b) show, respectively, annual averages for 2012 of photosynthetic and non-photosynthetic vegetation cover fractions (Guerschman et al., 2015), which we use to obtain the roughness correction $f_{v_{\mathrm{R}}}$ (label DPR) based on Raupach et al. (1993) (parameters c_λ=0.2, β=200, m=0.5, σ=1) shown in panel (c); panel (d) displays the annual average MODIS leaf area index and (e) static aerodynamic roughness length (Prigent et al., 2012), which we utilize for the roughness correction $f_{v_{\mathrm{M}}}$ (label DPM) after Marticorena and Bergametti (1995) (parameters $z_{\text{0s}}=(\mathrm{2}\times \mathrm{650}\times {\mathrm{10}}^{-\mathrm{4}}/\mathrm{30})$ cm, X=12 255 MacKinnon et al., 2004).

Figure 2 Normalized particle-size distributions (PSDs) based on D'Almeida (1987) (DA87, turquoise), Shao (2004) (S04, gray), and Kok (2011a) (as in K14, coral). The DA87 and K14 PSDs are invariant, while the S04 PSD varies in time and space. Shown for the latter are the PSDs corresponding to the 50th percentile (median; solid gray line) of the emission-weighted average diameters per model grid cell of annually accumulated dust emissions, framed by the PSDs belonging to the 5th and 95th percentiles (gray shading). The S04 PSDs were obtained from the S04 experiment presented in Sect. 4.1.

In the formulation from Marticorena and Bergametti (1995), f_v is given by

$$\begin{array}{}\text{(28)}& {\displaystyle }{\displaystyle }{f}_{v_{\mathrm{M}}}=\mathrm{1}-{\displaystyle \frac{\ln \left(z_{\mathrm{0}}/z_{\text{0s}}\right)}{\ln \left(\mathrm{0.7}{\left(X/z_{\text{0s}}\right)}^{\mathrm{0.8}}\right)}},\end{array}$$

where z₀ 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₀), f_v becomes smaller and the stress on the erodible surface decreases, reducing emission. The smooth roughness length z_0s is estimated as

$$\begin{array}{}\text{(29)}& {\displaystyle }{\displaystyle }{z}_{\text{0s}}={\displaystyle \frac{\mathrm{2}{d}_{\mathrm{c}}}{\mathrm{30}}},\end{array}$$

where $d_{\mathrm{c}}=\mathrm{650}\times {\mathrm{10}}^{-\mathrm{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₀ 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–near-infrared 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 the Moderate Resolution Imaging Spectroradiometer (MODIS) leaf area index (LAI). The calculation of z_0dyn is based on empirical relationships from Marticorena et al. (2006):

$$\begin{array}{}\text{(30)}& {\displaystyle }{\displaystyle }{z}_{\mathrm{0}\text{dyn}}=\left\{\begin{array}{ll}h\times {\mathrm{10}}^{\mathrm{1.3}\log \mathit{\lambda }+\mathrm{0.66}}& \mathit{\lambda }<\mathrm{0.041}\\ h\times {\mathrm{10}}^{-\mathrm{1.16}}& \mathit{\lambda }\ge \mathrm{0.041}\end{array}\right.,\end{array}$$

where h is the vegetation height and λ the roughness density (or frontal area index), defined as $\mathit{\lambda }=n\cdot a_{\mathrm{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=\mathit{\eta }/(\mathit{\pi }\cdot (D_{\mathit{\eta }}/\mathrm{2}{)}^{\mathrm{2}})$ (Pierre et al., 2012). With $a_{\mathrm{f}}=h\cdot D_{\mathit{\eta }}$, it follows that

$$\begin{array}{}\text{(31)}& {\displaystyle }{\displaystyle }\mathit{\lambda }=\mathrm{4}\mathit{\eta }{\displaystyle \frac{h}{D_{\mathit{\eta }}\mathit{\pi }}}.\end{array}$$

In Eq. (30), the dynamic roughness length increases with the characteristic height and density of the roughness elements. The influence of density is assumed to saturate above a sufficiently large value. In this implementation, η and h are assumed to scale with LAI:

$$\begin{array}{}\text{(32)}& {\displaystyle }{\displaystyle }h=h_{\text{max}}{\displaystyle \frac{\text{LAI}}{{\text{LAI}}_{\text{max}}}}.\end{array}$$

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 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_{\mathrm{0}}=max\left(z_{\mathrm{0}\text{stat}},z_{\mathrm{0}\text{dyn}}\right)$.

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 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), 40 cm (Sahel; Pierre et al., 2014), and 12 255 cm to extend its use to rougher vegetated surfaces (US; MacKinnon et al., 2004). The assumption of vegetation patches of 5 m in diameter was suggested as optimal for the Sahel (Pierre et al., 2012) but may be inadequate for other semi-arid regions. We note that this value can easily be modified in a static or dynamic way as soon as more detailed information becomes available. MONARCH uses maps of monthly LAI (actual year or climatology) and interpolates the monthly values to the day of simulation for each grid cell. The static roughness length and annual averages of the dynamic roughness length and the resulting drag partition correction, $f_{v_{\mathrm{M}}}$, using X from MacKinnon et al. (2004) are shown in Fig. 1d–f. Compared to $f_{v_{\mathrm{R}}}$, the correction $f_{v_{\mathrm{M}}}$ tends to be weaker with values typically above 0.35. Areas with low corrections generally coincide with those in $f_{v_{\mathrm{R}}}$, but $f_{v_{\mathrm{M}}}$ is smaller (weaker correction) in the Taklamakan and Gobi deserts and east/southeast of the Caspian Sea. When specifying X according to 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_{\text{bare,M}}=(\mathrm{1}-\mathit{\eta })\times (\mathrm{1}-{\mathit{\eta }}_{\text{snow}})\times (\mathrm{1}-{\mathit{\eta }}_{\text{br}})$ in combination with the drag partition from Marticorena and Bergametti (1995) and with $s_{\text{bare,R}}=(\mathrm{1}-{\mathit{\eta }}_{\text{snow}})\times (\mathrm{1}-{\mathit{\eta }}_{\text{br}})$ in 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 (Raupach et al., 1993; 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 in a grid cell where no emission occurs.

3.1.5 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 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 (Sect. 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 2 compares the PSDs based on D'Almeida (1987) (DA87), Shao (2004) (S04), and Kok (2011a) (as in K14). 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 fifth 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 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.

Figure 3 (a) Binary potential dust source areas defined based on $\text{FoO}(\text{DOD}>\mathrm{0.2})>\mathrm{0.025}$; light line patterns indicate anthropogenic dust sources using the method from Ginoux et al. (2012) considering cropland and pasture based on Klein Goldewijk et al. (2017); (b) binary dust source overlaid with FoO(DOD >0.2); (c) binary dust sources overlaid with the topographic source scaling function from Ginoux et al. (2001) without vegetation mask.

Table 6 Physical and optical particle properties available in MONARCH for eight particle-size bins: equivalent volume radius (r_v), effective radius (r_e), density (ρ_p), real and imaginary parts of the refractive index (ref_REAL, ref_IMAG), mass-extinction efficiency (MEE, [m² g⁻¹]), 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 Sect. 3.

Figure 4 Total annual dust emission (left), dust deposition (center), and annual average column dust load for 2012 using the configurations described in Sect. 4.1. Dust deposition includes gravitational settling, turbulent diffusion, and wet deposition from convective and non-convective precipitation. Shown are results for dust-particle diameters up to 20 µm.

Table 7 Statistical dust parameters of four global model simulations using the dust emission schemes MB95, G01-UST, S04, and K14 with the configurations described in Sect. 4.2. 〈MEE_g〉 and 〈MEE〉 are annual global averages of MEE (all sizes). 〈MEE_g〉 is calculated as the ratio of annual average grid-based DOD and dust load, whereas 〈MEE〉 is calculated from annual global average DOD and dust load. Parameters are for dust with particle diameters up to 20 µm.

3.1.6 Dust sources and lower boundary conditions for emission

In MONARCH, areas from which dust emission is possible are described using a map obtained from the climatological (for the years 2003–2015) frequency of occurrence (FoO) of the MODIS Deep Blue dust optical depth (DOD) greater than 0.2 (Hsu et al., 2004; Ginoux et al., 2012, see Sect. 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 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 Sect. 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 Ginoux et al. (2001) describes the sources as a function of topography. In practice, the topographic source term (S in Eqs. 2 and 3) enhances dust emission from enclosed basins in arid regions where soil particles have accumulated after fluvial erosion of the surrounding highlands (Fig. 3, bottom panel). Such scaling is part of the design of the G01 scheme and was found to 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 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 (agricultural) sources. For this purpose, the MODIS FoO-based map is linked with fractions of anthropogenic land use, following the approach described in Ginoux et al. (2012), but using an updated land-use data set (Klein Goldewijk et al., 2017). When considering cropland and pasture as anthropogenic dust sources, the main anthropogenic source regions are in the Sahel, India, China, and the United States (Fig. 3, top panel). Besides tagging natural and anthropogenic dust sources (Klose et al., 2018), MONARCH's tagging functionality can be adapted to track dust also from other predefined source origins (Kok et al., 2021b).

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 LAI data (Myeni et al., 2015) are 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). To aggregate soil texture data to model resolution, MONARCH utilizes a predominance approach; i.e., the predominant soil texture class in each grid cell is applied to the entire cell.

3.2 Dust transport and deposition

Dust transport and deposition in MONARCH have been thoroughly described in Pérez et al. (2011) and are 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 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 use 1, 1, 1, 1, 0.5, 0.3, 0.2, and 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 (subgrid-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 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 (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 and ice nuclei, are planned to be implemented in the future.

3.3 Radiation and optical properties

The model's radiation scheme is RRTM-G (Iacono et al., 2000, 2008). MONARCH allows multiple options for setting the dust microphysical properties. In the longwave (LW), we assume refractive indices from the Optical Properties of Aerosols and Clouds (OPAC) data set (Hess et al., 1998) and spherical particle shape. In the shortwave (SW), we use mineralogy-based refractive indices and non-spherical shapes. The multi-component Maxwell Garnett theory (Markel, 2016) is used to calculate refractive indices of internal mixtures of eight minerals (Gonçalves Ageitos et al., 2021a), whose size-resolved proportions are estimated based on the mineralogical atlas from Claquin et al. (1999) combined with the brittle fragmentation theory of Kok (2011a). The single-mineral refractive indices are taken from Scanza et al. (2015). We obtain size- and wavelength-dependent real and imaginary indices for each of the 28 soil types in the atlas and we take the median values. Note that the dependence of our refractive indices upon size is due to changes in mineralogy with size. Our median imaginary indices compare better than OPAC values (too absorbing) with recent chamber-based retrievals (Di Biagio et al., 2019), in situ aircraft measurements (Denjean et al., 2016), and ground-based remote sensing (Balkanski et al., 2007) of dust refractive index (Gonçalves Ageitos et al., 2021a). 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) with the dust single-scattering database of Meng et al. (2010). Table 6 summarizes key dust properties used in MONARCH.

4 Model performance and evaluation

A range of global model simulations was performed with MONARCH for 1 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 et al., 2019) products in terms of DOD.

4.1 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 four computational times steps; moist convection, microphysics, and wet scavenging routines were called every eight 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 1 year of spinup for soil moisture and 1 month of spinup for the dust fields before the 1-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. 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 3-hourly for the global runs.

Here we present results of global MONARCH simulations using the MB95, G01-UST, S04, and K14 dust emission schemes, a set of well-known and frequently used parameterizations. In all runs, we scaled soil moisture using $c_{{\mathrm{f}}_{\mathrm{1}}}=\mathrm{0.1}$ and applied the default soil moisture corrections listed in Table 4 (Sect. 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 intercomparison 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 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 by comparing monthly averages of modeled coarse DOD (size range 1.2–20 µm) for each experiment with the DOD obtained from MODIS (see Sect. 4.3.1 for more details) and minimizing the overall error (Cakmur et al., 2006). This calibration only removes the general global bias for each run and does not affect the spatiotemporal variability of the dust emission, transport, deposition, and interactions.

4.2 Dust emission and deposition

The total mass of dust emitted globally during 2012 was 3489, 3627, 5994, and 3739 Tg, respectively, for the MB95, G01-UST, S04, and K14 dust emission schemes. Correspondingly, the total dust deposition (dry/wet) obtained with the four schemes was 3442 (2435/1007), 3541 (2131/1410), 5893 (3929/1964), and 3664 (2215/1449) Tg. Dry dust deposition here includes both gravitational settling and turbulent diffusion. Wet deposition is due to convective and non-convective precipitation. The globally integrated annual average column dust load for the four configurations resulted as 29.0, 29.1, 40.6, and 31.4 Tg. Figure 4 shows the global spatial distribution of the total annual dust emission and deposition, as well as average column dust load for the four model runs. Values are summarized in Table 7.

The similarity in global dust emission between the MB95 and G01-UST schemes is a result of the scaling with the topographic source mask. Nevertheless, differences in the magnitude of dust emission are evident, in particular in the Middle East, central Asia, and Australia. Neither the S04 nor K14 scheme uses a preferential source function besides the binary treatment explained in Sect. 3.1.6. Hence, dust emissions are independent of this source function and differences to other experiments are more pronounced. Compared with the MB95 and G01-UST runs, for example, the Bodélé Depression in Chad does not stand out as much compared to the runs using the topographic source mask. Dust emissions in Asia extend over a larger area in the S04 and K14 runs and tend to be larger in North and South America. The S04 run shows decreased dust emissions in the eastern Sahara, whereas north African and Middle Eastern dust emissions are relatively homogeneous in the K14 run. Overall, 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 (Table 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.

Consistent with the differences in dust emission between the four runs, the annual total dust deposition and annual average dust load are similar in the MB95 and G01-UST runs, with pronounced individual source regions such as the Bodélé Depression. 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 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 northwestern 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 from Asian and other dust sources. MONARCH contributed to this study with a preliminary model version, which did not yet include all developments presented here. Our new results (Table A1) show that, compared to the results from Kok et al. (2021b), the contribution of North African dust sources is overestimated in the four presented configurations (range 55 %–71 % compared to 50 %), which suggests that we may underestimate the contribution from other source regions. This is most likely a result of larger vegetation/roughness element coverage in, for example, North America, which acts to suppress dust emission to a large extent in the coarse global model grid we used. However, besides the overestimation of North African sources, our model results agree with many of the aspects highlighted by Kok et al. (2021b): our global dust emission fluxes are in the range 3489–5994 Tg; Asian dust sources contributed to about 37 %–42 % to global emissions in three of four configurations; South American dust sources 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 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.

4.3 Dust optical depth

The global annual average of DOD is 0.034, 0.032, 0.041, and 0.035 in the MB95, G01-UST, S04, and K14 runs. This results in an average mass-extinction efficiency 〈MEE_g〉 of, respectively, 1.10, 1.11, 1.15, and 1.01 m² g⁻¹ for the four runs considering dust up to 20 µm in diameter (Table 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 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 (Hsu et al., 2013; 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.

Figure 5 Seasonally averaged MODIS Deep Blue DOD (left), MONARCH all-sky DOD_coarse at satellite overpass times co-located with MODIS DOD (middle), and clear-sky DOD_coarse at approximate satellite overpass times derived from 3-hourly model output from MONARCH (right). The model results were obtained for DOD_coarse averaged across the four model experiments. The seasonal averages were calculated with respect to the number of valid values per grid cell in the respective products. Area-weighted root mean square error (RMSE) and uncentered Pearson product-moment coefficient of linear correlation (r) between model and observations are indicated in the respective panels.

Figure 6 Seasonally averaged FoO of DOD >0.2, normalized by the number of days per season for MODIS Deep Blue (left), MONARCH all-sky DOD_coarse at satellite overpass time co-located with MODIS DOD (middle), and MONARCH clear-sky DOD_coarse at approximated satellite overpass times derived from 3-hourly model output (right). The FoO was calculated with respect to the number of days in the season. Area-weighted RMSE and uncentered Pearson product-moment coefficient of linear correlation (r) are indicated in the respective panels.

Figure 7 Globally averaged monthly global DOD_coarse (a) and FoO of DOD >0.2 (b) for MODIS (dashed green line) and MONARCH (DOD_coarse all-sky co-located with MODIS observations) (green solid line). The shading indicates the range of DOD_coarse across the four MONARCH experiments, which are also shown. Pearson correlation coefficients (r) and RMSE are given in the figure for both the experimental average as well as the individual runs.

Figure 8 Comparison of 3-hourly DOD between MONARCH (average (turquoise line) and standard deviation (shading)) and AERONET direct-Sun V3 level 2.0 for selected stations covering Cabo Verde and the Canary Islands (Capo Verde, Santa Cruz de Tenerife), the Sahara and Sahel (Ouarzazate, Tamanrasset, Cinzana, Banizoumbou), the Middle East (Eilat, Solar Village, Masdar Institute), Asia (Karachi, Issyk-Kul, Dalanzadgad), Europe (Granada), southern Africa (Henties Bay), Australia (Birdsville), and North and South America including the Caribbean (Railroad Valley, CASLEO, Ragged Point). The direct-Sun DOD is filtered for dust aerosol using Ångström exponent (AE) <0.3 (filled circles). Records which do not meet the AE criteria are less likely to be associated with dust and are shown as open circles. The Pearson correlation coefficients (corr) and RMSE are given in each panel.

4.3.1 Comparison of modeled DOD with MODIS Deep Blue

We estimate DOD from MODIS using daily AOD and SSA at 550 nm, and Ångström exponent (AE) of the Deep Blue Collection 6 Level 2 MODIS products (Hsu et al., 2013; Sayer et al., 2013) from the Aqua platform at 0.1^∘ resolution (Ginoux et al., 2010). As in Ginoux et al. (2012) and Pu and Ginoux (2018), DOD is estimated from AOD using a continuous function of AE (Anderson et al., 2005). Pu and Ginoux (2018) estimated an error of $\pm (\mathrm{0.08}+\mathrm{0.52}\text{DOD})$ for the DOD derived from MODIS Aqua.

To enable a direct comparison between MODIS satellite observations and MONARCH results independent of the model output frequency, MONARCH internally diagnoses the all-sky DOD for a given satellite overpass time for each day. The sampling of the satellite overpass time follows Quaas (2011, https://www.yumpu.com/en/document/read/19547757/modis-simulator-software-euclipse, last access: 8 October 2021) and is done based on a longitude-based local time (LLT). We assume 13:30 LLT as the overpass time of MODIS Aqua. Actual overpass times vary and may deviate slightly from this nominal time. The same diagnostic is also available for MODIS Terra (nominal overpass time 10:30 LLT). Other polar satellite overpasses can be implemented easily. For model evaluation, the MONARCH-modeled satellite DOD is additionally co-located in space and time with the satellite observations, i.e., grid cells for which the MODIS data contain missing values because of clouds are filtered from the MONARCH data for each day. MONARCH also estimates the DOD under clear-sky conditions (i.e., without clouds) based on the modeled cloud fraction and a coin-flipping method. The clear-sky DOD is currently diagnosed at the model output times (in contrast to the satellite overpass times available for the all-sky DOD). For that reason, we apply a post-processing and sample the clear-sky DOD for the output time closest to the satellite overpass time (subsequently termed “approximate satellite overpass time”). For comparison 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 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.

Figure 5 shows seasonal averages of MODIS DOD (left) and modeled global all-sky co-located DOD_coarse at satellite overpass times (center) and clear-sky DOD_coarse at approximate satellite overpass times (right) averaged across the four global MONARCH runs. The spatial patterns of observed DOD in northern Africa and the Middle East are well represented in MONARCH throughout the year. Distinct features are high DOD in the Bodélé area (somewhat overestimated mainly in MAM and SON) with elevated levels also toward the south/southwest in MAM and toward the west/northwest in JJA; and increased AOD along the eastern coast of the Arabian Peninsula in MAM and in its southern part in JJA. The spatiotemporal evolution of DOD in central Asia also agrees well between MONARCH and MODIS, with relatively low values in SON and DJF and increased DOD in particular in the Thar and Registan deserts in MAM and JJA. The DOD north of the Aral Sea is underestimated in MONARCH compared to MODIS in JJA. Likewise, DOD in the Taklamakan Desert is lower in MONARCH compared to MODIS in DJF and particularly in MAM. This may be related to the pronounced topographic features in the area, which are difficult to resemble at coarse model resolution. The DOD in Australia is relatively low throughout the year in both MODIS and MONARCH, but with areas of slightly increased DOD in northeastern Queensland in DJF, which are underrepresented in MONARCH. Similarly, a somewhat higher DOD in southern Africa in JJA and SON, in South America in SON and DJF, and in North America in MAM, JJA, and SON are underestimated in MONARCH. Global pattern correlations between MODIS and the MONARCH average range between 0.69 and 0.85 for all-sky co-located DOD_coarse (area-weighted RMSE between 0.09 and 0.14) with stronger correlations in MAM and JJA than in other seasons. Note that the algorithm used to derive DOD from MODIS AOD cannot perfectly discriminate dust from other aerosols. This may lead to an overestimation of DOD in areas in which dust is not the dominant aerosol type and where other aerosols are present. Seasonal averages of modeled all-sky co-located DOD_coarse for each individual model run are shown in Fig. B1.

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

Figure 6 shows the FoO of DOD >0.2, normalized by the number of days in each season, again for MODIS DOD (left), as well as modeled all-sky co-located (center) and clear-sky (right) DOD_coarse. The spatial patterns of the observed FoO are very well captured by the MONARCH runs, in particular for the all-sky co-located FoO (pattern correlations between 0.81 and 0.92), for key dust sources in northern Africa (e.g., Erg of Bilma/Bodélé Depression, Grand Erg Oriental/Erg Chech, El Djouf desert), the Middle East (e.g., Rub' al Khali and Nefud deserts), and central and eastern Asia (e.g., Registan/Thar and Karakum deserts, Taklamakan, and Gobi deserts). Discrepancies in the all-sky co-located FoO magnitude between the model ensemble and MODIS depend on the season: the FoO is slightly underestimated in deserts east of the Caspian Sea in MAM and JJA, in South America and South Africa in JJA and SON, in Australia in SON and DJF, and in North America throughout the year. FoO values are slightly overestimated in the Arabian Peninsula in JJA and SON, and eastern north Africa in MAM and JJA. The modeled seasonal all-sky co-located FoO for each individual model run is shown in Fig. B2.

The FoO obtained from the modeled clear-sky DOD_coarse is generally larger than that obtained from the all-sky co-located DOD_coarse. Due to the normalization with the number of days in the season for calculation of the FoO, differences in the frequency and location of clouds in MONARCH and MODIS directly impact the resulting FoO. Over dust source regions, MONARCH produces considerably fewer cloud pixels in its clear-sky product and hence a larger number of valid data values than are in the MODIS observations (not shown). As a result, the clear-sky FoO is based on a larger number of valid (high-DOD_coarse) values and is therefore larger than the MODIS and MONARCH all-sky co-located FoOs. With 0.77–0.88, clear-sky FoO pattern correlations are slightly weaker than those obtained for the all-sky co-located FoO.

Figure 7 shows global averages of monthly DOD and FoO from MODIS and the all-sky co-located DOD_coarse from MONARCH. Globally, the DOD obtained with MODIS is reproduced well with MONARCH in all four experiments (Pearson correlation coefficient of 0.98 for the experimental average; between 0.86 and 0.97 for the individual runs). The DOD range across the four experiments is relatively similar throughout the year with the spread being the largest during the Northern Hemisphere peak dust season in March. The MB95 experiment contributes the largest DOD in January until March and the lowest during much of the remaining year, whereas the K14 experiment shows opposite behavior. The G01 and S04 experiments are intermediate between the other two runs and best resemble the monthly global DOD for the given configurations. The correlation between the FoO of DOD >0.2 from MODIS and from MONARCH is also very high (0.95 for the MONARCH average; between 0.84 and 0.97 for the individual runs). Whereas the results from the four experiments are very similar from November until April, the variability increases during the other months with the largest range in September. The results for the four individual runs are qualitatively similar to those for DOD, with the MB95 and K14 experiments providing, respectively, the lower and upper frames from approximately April until November, and the G01 and S04 runs being intermediate.

Figure 9 (a–c) Scatter plots of 3-hourly, daily, and monthly DOD estimated from AERONET direct-Sun V3 level 2.0 (AE <0.3) and total DOD from MONARCH (average across runs) averages for the stations listed in Appendix C; panels (d–f) are the same as (a–c) but for AERONET O'Neill V3 level 2.0 and coarse (diameters 1.2–20 µm) MONARCH DOD. The Pearson correlation coefficient (corr) and RMSE are given in the plot.

Figure 10 Annual average longwave, shortwave, and total direct radiative effect [W m⁻²] at the surface (SFC) and the top of the atmosphere (TOA) obtained from MONARCH as average across the four runs. Global average values are indicated in each panel.

Figure 11 Global averages of MONARCH-derived DRE [W m⁻²] for each run at the SFC and TOA for shortwave and longwave radiation and the total (shortwave and longwave). Shown are the relative contributions per bin, normalized with the absolute value of the DRE for all bins (left) and the DRE for all particle sizes (right). The diameter ranges and effective radii of each bin are given in, respectively, Sect. 3 and Table 6. The DRE results are also given in Appendix D.

4.3.2 Comparison of modeled DOD with AERONET

AERONET is a global network of ground-based solar photometer stations (Holben et al., 1998; O'Neill et al., 2003; Giles 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 (cloud screened), 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 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 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 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 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 (submicron) and coarse (supermicron) 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 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 time series of 3-hourly DOD from AERONET direct Sun and MONARCH for 18 selected stations in the vicinity of dust sources around the globe: four stations in northern Africa and two stations in the typical dust outflow region west of northern Africa; three stations in the Middle East; three stations in Asia; and one station each in Europe, southern Africa, Australia, North America, South America, and the Caribbean. In addition, Fig. 9 compares modeled and observed DOD and DOD_coarse for the 57 stations listed in Appendix C. The station locations are shown in Fig. C1. The time series demonstrate an overall good agreement between the average modeled and observed DOD where the temporal variability is mostly reproduced with discrepancies for individual DOD peaks. Consistent with MODIS DOD, the AERONET DOD tends to be small at the stations in southern Africa, Australia, and North and South America, and AE is often not below 0.3, i.e., at least part of the DOD is likely due to aerosols other than dust. Correlations between MONARCH and AERONET are smaller for these stations, because MONARCH DOD represents pure dust and because a mismatch between individual peaks receives more weight if the number of dust episodes is small. The Pearson correlation coefficients for all other stations range between around 0.3 (Masdar Institute, Solar Village) and 0.7 (IER Cinzana, Karachi, Ragged Point).

Taking into account the entire station list (Appendix C), the correlation is 0.61 with a RMSE of 0.31 for the total DOD and, 0.71 with an RMSE of 0.14 for DOD_coarse, based on the 3-hourly MONARCH data (Fig. 9). For total DOD, the correlation remains fairly constant when comparing daily and monthly instead of 3-hourly values. The RMSE decreases slightly with an increasing averaging period as then discrepancies for individual peaks become less relevant. A similar behavior 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 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.

4.4 Direct radiative effect

Figure 10 shows the shortwave, longwave, and total direct radiative effect (DRE) at the surface (SFC) and the top of the atmosphere (TOA). The longwave DRE at the SFC is positive and most pronounced in the dust belt ranging across northern Africa, the Middle East, and southwestern Asia. In the dustiest areas of northern Africa and the Middle East, the longwave SFC DRE reaches between 10 and 20 W m⁻². This strong sensitivity to the presence of dust is a result of the low atmospheric moisture content in this area. High near-surface atmospheric temperatures enhance the longwave downwelling radiation (Miller et al., 2014). The longwave DRE at the TOA is smaller and typically under 5 W m⁻² due to the opposing effects of scattering and absorption by dust at the TOA. The shortwave DRE is strongly negative at the SFC, with values exceeding −20 W m⁻² in the main dust regions. At the TOA, the shortwave DRE is slightly negative in most areas but slightly positive in some of the northern African dust sources related to the relatively bright underlying desert surface. This results in a negative total DRE at the SFC, with the largest (negative) values in the Sahel, the eastern Atlantic, and the Arabian Sea. At the TOA, the total DRE is positive around the main north African dust sources and slightly negative/neutral elsewhere. The globally averaged DRE at the SFC and TOA for all particle sizes as well as from each bin (relative contribution with respect to the DRE for all sizes) is summarized in Fig. 11 and listed in Table D1.

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 about 10 Wm⁻² in these regions. In contrast, Balkanski et al. (2007), for example, obtained a total DRE at the TOA of up to 20 Wm⁻² 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 was found by Gonçalves Ageitos et al. (2021a) for our refractive indices. Moreover, Miller et al. (2006) obtained a global average total DRE at the TOA of −0.39 Wm⁻² 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⁻² 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⁻² 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.

5 Conclusions and outlook

We presented the description of mineral dust in the Multiscale Online Non-hydrostatic AtmospheRe CHemistry model (MONARCH) version 2.0. MONARCH contains multiple state-of-the-art options to represent dust emissions on global and regional scales, ranging from more simplified to more complex parameterizations based on physical processes. We tested and evaluated a set of four global model configurations for the year 2012. Comparison with observations of dust optical depth from MODIS and AERONET showed a good model reproduction of key features of the observed dust cycle. Global annual dust emissions ranged between around 3500 and 6000 Tg. Differences in modeled dust emissions between the four configurations were mainly driven by the dust source description (use of a preferential source mask or not) and the particle-size distributions at emission. Dust deposition ranged between about 3450 and 5900 Tg in 2012 globally, yielding an average dust load of 29–41 Tg. The smaller range of simulated load among experiments is due to the shorter lifetime of the coarse particles included in the S04 scheme that exhibits larger emission. The total direct radiative effect obtained from the MONARCH simulations is slightly negative at the surface in dust transport regions. At the top of the atmosphere, the total direct radiative effect is positive near the main north African dust sources and slightly negative/neutral elsewhere.

The multifaceted options of MONARCH and its dust component, combined with an advanced workflow management for use in high-performance computing environments, makes it a powerful and versatile tool applicable for process studies, operational forecasting, and climate research. In the following, we outline a few ongoing activities related to the MONARCH dust component to demonstrate its capabilities.

Dust ensemble runs can be generated with MONARCH by utilizing the diverse model configurations and by perturbing model parameters related to, for example, surface winds, soil humidity, and the spatial distribution of dust emission, which are deemed to be uncertain. In Di Tomaso et al. (2017), perturbations were applied to the threshold friction velocity and the dust emission flux per size bin. In Escribano et al. (2021), different dust emission schemes – or linear combinations thereof – are used by different ensemble members. Combined meteorology and emission perturbations were shown to be necessary to produce sufficient ensemble spread in (dust) aerosol outflow regions (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 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).

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

The combination of different vegetation input data sets, drag partition approaches, and the source tagging capability allows us to represent the seasonal vegetation dynamics and provides an ideal basis to investigate the importance of dust from anthropogenic (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).

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.

Table A1 lists the relative contributions from different source regions to global annual dust emissions are summarized in Table A1. The dust source definitions are as in Kok et al. (2021b).

Table A1 Percentage contributions of different dust source regions (definitions as in Kok et al., 2021b) to global annual dust emissions with the four presented configurations.

Figure A1 Percent contribution of dust emission (left column) and deposition (right column) to their respective global and annual totals for the four experiments described in Sect. 4.1.

Appendix B: Dust optical depth and its frequency of occurrence

Figures B1 and B2 show the DOD and FoO of DOD >0.2 for MODIS and MONARCH (all-sky DOD_coarse at satellite overpass time co-located with MODIS) for the four MONARCH experiments using the MB95, G01-UST, S04, and K14 dust emission schemes.

Figure B1 Seasonally averaged MODIS Deep Blue DOD (left) and MONARCH all-sky DOD_coarse at satellite overpass times co-located with MODIS DOD for the MB95, G01-UST, S04, and K14 runs. The seasonal averages were calculated with respect to the number of valid values per grid cell in the MODIS product. Area-weighted RMSE and uncentered Pearson product-moment coefficient of linear correlation (r) are indicated in the respective panels.

Figure B2 Seasonally averaged FoO of DOD >0.2, normalized by the number of days per season for MODIS Deep Blue (left) and MONARCH all-sky DOD_coarse at satellite overpass times co-located with MODIS DOD for the MB95, G01-UST, S04, and K14 runs. The FoO was calculated with respect to the number of days in the season. Area-weighted RMSE and uncentered Pearson product-moment coefficient of linear correlation (r) are indicated in the respective panels.

Appendix C: Comparison with AERONET

The AERONET stations used for comparison with MONARCH and to obtain the global calibration factor are listed in Table 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.

Figure C2 compares 3-hourly DOD estimated from AERONET direct-Sun V3 level 2.0 (AE <0.3) and AERONET O'Neill V3 level 2.0 with, respectively, 3-hourly total DOD and DOD_coarse from MONARCH for all four experiments (MB95, G01, S04, and K14) taking into account the stations listed in Table C1.

Table C1 List of AERONET stations used for comparison with MONARCH results.

Figure C1 AERONET stations (direct-Sun V3 level 2.0) available in 2012. The station subset used for comparison with MONARCH are shown in turquoise, whereas all other stations are marked in red.

Figure C2 Scatter plots (top) of 3-hourly DOD estimated from AERONET direct-Sun V3 level 2.0 (AE <0.3) and total DOD from MONARCH by experiment (MB95, G01, S04, and K14) for the stations listed in Appendix C; the bottom is the same as the top but for AERONET O'Neill V3 level 2.0 and coarse (diameters 1.2–20 µm) MONARCH DOD. The Pearson correlation coefficient (corr) and RMSE are given in the plot.

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.

Table D1 Global averages of MONARCH-derived DRE [W m⁻²] (all particle sizes and relative contribution per bin) for each run at the SFC and TOA for shortwave and longwave radiation and the total (shortwave and longwave). Results are visualized in Fig. 11. The diameter ranges and effective radii of each bin are given in, respectively, Sect. 3 and Table 6.

Appendix E: List of symbols

a	coefficient in soil moisture correction from Fécan et al. (1999)
α	coefficient in soil moisture correction from Klose et al. (2014)
αq	vertical-to-horizontal-flux ratio
b	coefficient in soil moisture correction from Fécan et al. (1999)
β	coefficient in soil moisture correction from Klose et al. (2014)
βR	coefficient in drag partition correction from Raupach et al. (1993)
βR	ratio of roughness element to surface drag coefficients
Cα	coefficient in K14 scheme
Ce	coefficient in K14 scheme, dependent on u*st
$c_{{\mathrm{f}}_{\mathrm{1}}}$	soil moisture scaling factor
$c_{{\mathrm{f}}_{\mathrm{2}}}$	soil moisture scaling factor
cλ	coefficient
cQ	coefficient
CG01	dimensional factor in G01 scheme
cthr	scaling coefficient
cy	coefficient in S01, S04, and S11 schemes
dc	coarse particle diameter in drag partition of Marticorena and Bergametti (1995)
Dη	diameter of vegetation patches
di0, di1	lower and upper diameter limits of the particle-size bin corresponding to di
di	dust-particle diameter
ds	saltation diameter
dmax	upper diameter limit of dust in MONARCH; currently dmax=20 µm
ηci	aggregated dust fraction at diameter di
ηfi	total dust fraction at diameter di
ηmi	free dust fraction at diameter di
ηbr	bedrock cover fraction
ηsnow	snow/ice cover fraction
η	vegetation cover fraction
fclay	clay fraction
fv	drag partition correction
F	vertical dust emission flux
$f_{w_{\mathrm{B}\mathrm{wet}}}$	constant in soil moisture correction from Belly (1964)
fw	soil moisture correction factor
γ	function in S04 scheme (Eq. 5)
g	gravitational acceleration [m s−2]
hmax	maximum annual vegetation height (relevant for dust emission) for use in drag partition
h	vegetation height
hw	function in soil moisture correction from Klose et al. (2014)
κ	coefficient in S01 and S04 schemes
λ	roughness density (frontal area index)
m	coefficient in drag partition correction from Raupach et al. (1993)
mΩ	soil mass removed by a saltation impact; mΩ=ρbΩ
mp	particle mass
$m_{{\mathrm{p}}_{\mathrm{s}}}$	mass of particles with diameter ds
Msand	sand fraction
n	number of vegetation patches
Ω	soil volume removed by a saltation impact
pf	fully dispersed particle-size distribution [m−1]
pm	minimally dispersed particle-size distribution [m−1]
ps	airborne sediment particle-size distribution in the S11 scheme [m−1]
ψ*	function in K14 scheme
ψs	saturation capillary pressure head [m]

P	soil plastic pressure [Pa]
Q	horizontal saltation flux
Qs	saltation flux of a particle with diameter ds
re	effective radius [m]
ρa0	reference air density; ρa0=1.225 kg m−3
ρa	air density [kg m−3]
ρbd	bulk density of dry soil [kg m−3]
ρb	soil bulk density (approximately 1000 kg m−3)
ρl	water density [kg m−3]
ρpa	average soil particle density; ρpa=2500 kg m−3
ρp	particle density [kg m−3]
ρps	saltator density [kg m−3]
rv	equivalent volume radius [m]
S	dust source scaling function
σm	saltation bombardment efficiency; ${\mathit{\sigma }}_m=m_{\mathrm{\Omega }}/m_{{\mathrm{p}}_{\mathrm{s}}}$
${\mathit{\sigma }}_{{\mathrm{p}}_i}$	${\mathit{\eta }}_{\mathrm{mi}}/{\mathit{\eta }}_{\mathrm{fi}}$
σv	ratio of roughness-element basal to frontal area
sp	particle-size fractions in G01 scheme
sbare	bare soil fraction
${\mathit{\tau }}_{\mathrm{s}}^{\prime }$	average surface stress on exposed area
${\mathit{\tau }}_{\mathrm{s}}^{\prime \prime }$	maximum surface stress on exposed area
τ	total stress
θr	volumetric air-dry residual soil moisture content [m3 m−3]
θs	volumetric saturation soil moisture content [m3 m−3]
θ	volumetric soil moisture content [m3 m−3]
u*	friction velocity [m s−1]
u*st0	minimum standardized threshold friction velocity; $u_{*\mathrm{st}\mathrm{0}}\approx \mathrm{0.16}$ m s−1
u*st	standardized threshold friction velocity; $u_{*\mathrm{st}}=u_{*\mathrm{t}}\sqrt{{\mathit{\rho }}_{\mathrm{a}}/{\mathit{\rho }}_{\mathrm{a}\mathrm{0}}}$ [m s−1]
u*t0	threshold friction velocity for dry conditions
$u_{*{\mathrm{t}}_{\mathrm{d}\mathrm{0}}}$	minimum dry threshold friction velocity
$u_{{\mathrm{t}}_{\mathrm{d}\mathrm{0}}}$	minimum dry threshold wind velocity
u10m	10 m wind speed [m s−1]
u*NMMB	friction velocity provided by the atmospheric model
$u_{*{\mathrm{t}}_{\mathrm{dry}}}$	model threshold friction velocity for dry conditions (including optional scaling)
u*t	threshold friction velocity [m s−1]
ut	threshold wind speed [m s−1]
w	gravimetric soil moisture content [%]
wr	gravimetric air-dry residual soil moisture content [%]
X	parameter in drag partition from Marticorena and Bergametti (1995)
z0	aerodynamic roughness length
z0s	smooth aerodynamic roughness length
z0dyn	dynamic aerodynamic roughness length
z0stat	static aerodynamic roughness length