Dust Constraints from joint Observational-Modelling-experiMental analysis (DustCOMM): comparison with measurements and model simulations

Mineral dust is the most abundant aerosol species by mass in the atmosphere, and it impacts global climate, biogeochemistry, and human health. Understanding these varied impacts on the Earth system requires accurate knowledge of dust abundance, size, and optical properties, and how they vary in space and time. However, current global models show substantial biases against measurements of these dust properties. For instance, recent studies suggest that atmospheric dust is substantially coarser and more aspherical than accounted for in models, leading to persistent biases in modelled impacts of dust on the Earth system. Here, we facilitate more accurate constraints on dust impacts by developing a new dataset: Dust Constraints from joint Observational-Modelling-experiMental analysis (DustCOMM). This dataset combines an ensemble of global model simulations with observational and experimental constraints on dust size distribution and shape to obtain more accurate constraints on three-dimensional (3-D) atmospheric dust properties than is possible from global model simulations alone. Specifically, we present annual and seasonal climatologies of the 3-D dust size distribution, 3-D dust mass extinction efficiency at 550 nm, and two-dimensional (2-D) atmospheric dust loading. Comparisons with independent measurements taken over several locations, heights, and seasons show that DustCOMM estimates consistently outperform conventional global model simulations. In particular, DustCOMM achieves a substantial reduction in the bias relative to measured dust size distributions in the 0.5–20 µm diameter range. Furthermore, DustCOMM reproduces measurements of dust mass extinction efficiency to almost within the experimental uncertainties, whereas global models generally overestimate the mass extinction efficiency. DustCOMM thus provides more accurate constraints on 3-D dust properties, and as such can be used to improve global models or serve as an alternative to global model simulations in constraining dust impacts on the Earth system.

1 Introduction

Even though mineral dust accounts for a substantial fraction of the total mass of aerosol particles in the atmosphere and produces important impacts on the Earth system, global models are unable to accurately reproduce dust abundance, size, and optical properties (Kinne et al., 2006; Huneeus et al., 2011). Model difficulties in reproducing these atmospheric dust properties are largely associated with their inability to accurately simulate important dust processes, such as dust emission, transport, and deposition (e.g. Ginoux et al., 2001; Shao, 2001; Zender et al., 2003; Huneeus et al., 2011; Kok et al., 2017). Dust aerosols are emitted from source regions such as the Saharan, Middle East, and Asian deserts and are deposited after they are transported for thousands of kilometres (Duce et al., 1980; Prospero et al., 1981; Weinzierl et al., 2017). Their abundance and long-range transport allow them to play a significant role in the processes that impact global climate (Boucher et al., 2013), biogeochemistry (e.g. Mahowald et al., 2008, 2009; Ito et al., 2019), and human health (e.g. Giannadaki et al., 2014). Specifically, dust affects global climate directly by influencing the amount of radiation that can reach or leave the atmosphere and the surface (Haywood et al., 2003; Kok et al., 2017) or indirectly by changing the amount, reflectivity, and lifetime of clouds (e.g. Lohmann and Diehl, 2006; Doherty and Evan, 2014; Amiri-Farahani et al., 2017). In addition, dust also impacts global biogeochemistry through deposition of iron and phosphorous-rich micro-nutrients (Mahowald et al., 2008, 2009; Ito et al., 2019), both of which are linked to the ability of ocean and land ecosystems to absorb atmospheric carbon dioxide (e.g. Watson et al., 2000; Blain et al., 2007). Finally, dust particles are easily inhaled by humans, with smaller dust particles penetrating deeply into the lungs and leading to cardiopulmonary disease, lung cancer, and eventually death (e.g. Giannadaki et al., 2014). Therefore, obtaining accurate constraints on the many impacts of dust on the Earth system requires accurate knowledge of the sizes, abundance, and optical properties of atmospheric dust particles (Mahowald et al., 2014).

Uncertainties in dust aerosol properties directly translate into uncertainties in estimating their impact on the Earth system, such as dust radiative impacts (e.g. Huneeus et al., 2011; Zhao et al., 2013; Albani et al., 2014). Several studies have associated a large part of these uncertainties to the uncertainty in simulating the dust size distributions (e.g. Huneeus et al., 2011; Kok, 2011a; Evan et al., 2014). Specifically, global models simulate too much fine-mode dust ($\sim D\le \mathrm{2.5}\mathrm{\mu }\mathrm{m}$) and too little coarse-mode dust ($\sim D\ge \mathrm{5}\mathrm{\mu }\mathrm{m}$), both at emission and during transport in the atmosphere (e.g. Kok, 2011a; Kok et al., 2017). This bias is particularly problematic because fine dust predominantly cools the climate system by extinguishing shortwave (SW) radiation, whereas coarse dust warms it by also extinguishing longwave (LW) radiation (e.g. Tegen and Lacis, 1996; Dufresne et al., 2002). Whereas previous modelling studies affected by the size bias found that the combined (SW + LW) effect of dust is to cool the climate system (e.g. Tegen and Lacis, 1996; Tegen et al., 1996; Colarco et al., 2014), it is unclear whether the dust LW warming effect may overcome the dust SW cooling effect when the underestimation of coarse-mode particles is corrected (Kok et al., 2017). Since the dust radiative effect is sensitive to the representation of size distribution in global models, constraining the dust size distribution, and how it varies spatially, is thus important.

In addition to the sensitivity of dust size distribution, dust radiative effects are also sensitive to the shape of dust particles (e.g. Kalashnikova and Sokolik, 2004). Global models generally assume that dust particles are spherical (Ginoux et al., 2001; Miller et al., 2006; Huneeus et al., 2011), even though observations suggest that they are highly non-spherical (Okada et al., 2001; Potenza et al., 2016). This idealization in the representation of dust shape in global models is used to simplify model physics (e.g. Miller et al., 2006) and the calculation of their optical properties, but recent studies show that neglecting the asphericity of dust in models causes an underestimation of about 30 % of dust aerosol optical depth (AOD) or extinction produced per unit mass of dust (Potenza et al., 2016; Kok et al., 2017). This is largely caused by the greater surface-to-volume ratio of non-spherical particles compared to that of equal-volume spherical particles (e.g. Kalashnikova and Sokolik, 2002, 2004). The assumption of spherical dust in climate models is also problematic because the resulting underestimation of dust AOD largely masks the positive bias associated with the fine dust particles in models, which results in an overestimation of dust AOD and extinction at remote regions when the dust emissions are scaled to match the observation of AOD near the source regions (e.g. Kok et al., 2017). Hence, to properly constrain dust impacts on radiation, observational constraints must be applied to both the dust size distribution and dust shape.

Global model simulations of the global dust cycle are thus subject to numerous important biases, which have obscured a detailed understanding of the impacts of dust on the Earth system. To address the problem of size and shape biases in model simulation of dust properties, we propose a methodology to more accurately obtain three-dimensional (3-D) dust properties than is possible from global model simulations alone. Specifically, we propose a new product called the Dust Constraints from joint Observational-Modelling-experiMental analysis (DustCOMM), which combines an ensemble of global model simulations with observational and experimental constraints on dust size distribution and shape. DustCOMM builds on the results from Kok et al. (2017); however, unlike the globally averaged results obtained in Kok et al. (2017), our product constrains the climatology of 3-D global atmospheric dust properties, and it is provided on seasonal and annual timescales. Below, Sect. 2 describes the details of the methodology as well as the data used. In Sect. 3, we present the constrained spatial distribution of the dust size distribution, mass extinction efficiency, and atmospheric dust loading, which we evaluate using independent in situ measurements of dust size distributions and mass extinction efficiencies. Section 4 discusses some discrepancies between DustCOMM and measurements, the impact of dust asphericity on the DustCOMM product, and the possible use of DustCOMM to improve estimates of dust impacts in the global model simulations. Section 5 summarizes the paper. Finally, we note that all the DustCOMM dust aerosol properties (dark shaded boxes in Fig. 1) presented in this study are publicly available (Adebiyi et al., 2019a).

Figure 1 Schematic of the key steps to obtain the DustCOMM products (dark shaded boxes): constraints on the 3-D dust size distribution, 3-D mass extinction efficiency, and 2-D atmospheric loading. The variables are a function of the following: x–y–z (three-dimensional field), x–y (two-dimensional field), S (seasonally resolved), D (size-resolved), and σ (includes uncertainties). The variables in the light grey boxes are obtained from Kok et al. (2017). See Sect. 2 for details.

2 Data and methodology

We describe here all the steps we took to obtain the DustCOMM products. First, we use three sets of input datasets to create the DustCOMM products (Fig. 1): (1) the constrained globally averaged data from Kok et al. (2017); (2) six model simulations of size-resolved dust mass concentrations, from which we estimate the modelled dust size distribution; and (3) reanalysis datasets of the dust aerosol optical depth. We focus here only on describing the model simulations (Sect. 2.1) and the reanalysis products (Sect. 2.2), as details of the in situ measurements used to constrain the globally averaged datasets are described in Kok et al. (2017). Second, we describe the framework used to obtain DustCOMM dust size distribution, mass extinction efficiency, and atmospheric dust loading (Sect. 2.3). Finally, we describe the independent measurements we use to evaluate DustCOMM dust size distribution and the dust mass extinction efficiency in Sect. 2.4.

2.1 Model simulations

We use model outputs of dust aerosol properties from six leading atmospheric global models, namely the Goddard Institute for Space Studies (GISS) ModelE atmospheric general circulation model (Miller et al., 2006); the Weather Research and Forecasting model coupled with Chemistry updated by the University of Science and Technology of China (USTC) suitable for quasi-global simulation (WRF-Chem; Zhao et al., 2010, 2013; Hu et al., 2016); the Community Earth System Model (CESM; Hurrell et al., 2013); the Goddard Earth Observing System coupled with Chemistry (GEOS-Chem; see Kok et al., 2017); the ARPEGE-Climat model from the Centre National de Recherches Météorologiques Earth system model (Michou et al., 2015); and the Integrated Massively Parallel Atmospheric Chemical Transport (IMPACT; Ito and Kok, 2017, and references therein) model. We use the different simulations from global climate and chemical transport models between 2004 and 2008 (except for WRF-Chem and IMPACT, which are 2007–2016 and 2004, respectively) to capture the general model uncertainties that are associated with the dust emission, transport, and deposition processes. The GISS, CESM, and GEOS-Chem model simulations are described in Kok et al. (2017) and the references therein (see Sect. 5 of their supplementary document). Here, we supplement these simulations with three additional simulations from the WRF-Chem, ARPEGE-Climat, and IMPACT models. The WRF-Chem model simulation represents an updated USTC version of the one used in Kok et al. (2017). Further details of these three additional model simulations are thus given in the Supplement.

We obtain the spatially varying dust size distribution from each of the six model simulations, which we use to define the spatial variability of the DustCOMM dust size distribution (see Sect. 2.3.1). Specifically, the spatial variability of DustCOMM dust size distribution follows the ensemble of the six model simulations. We summarize the particle bin ranges, time periods, spatial resolutions, as well as the meteorology used for each model simulation of the dust size distribution in Table 1. All the models use discrete bins that represent the dust particles up to about 10 µm, except for the GISS, ARPEGE-Climat, and IMPACT models, which extend beyond the 10 µm diameter limit. Four of the models – WRF-Chem, CESM, ARPEGE-Climat, and IMPACT – have a lower diameter limit smaller than 0.2 µm. For consistency, we set the lower diameter limits for all the model simulations to the common diameter of 0.2 µm and correct the upper diameter limit to 20 µm, following the procedures we describe later in Sect. 2.3.1. In addition, since the time periods are different for the available model dataset (Table 1), we focus on annual and seasonal climatologies, which we obtain here from the monthly means of the model outputs.

Table 1 Details of model simulations used in this study. Shown are the particle bin ranges, time periods of simulations, the spatial resolutions, the meteorology used, and the relevant model reference. We interpolate all model simulations to 2.5^∘ by 1.9^∘ to facilitate comparison and consistency with other datasets. See Sect. 2 for details.

In order to test our hypothesis that integrating experimental and observational constraints on dust size and shape distributions can constrain 3-D dust properties more accurately than is possible from model simulations alone, we obtain a model ensemble of 3-D dust size distribution and mass extinction efficiency and 2-D dust column loading. To do so, we interpolated seasonal and annual climatologies of these dust properties to a common resolution of approximately 2.5^∘ by 2.0^∘ spatial resolution, with 35 levels from the surface to 100 hPa. In addition, we correct each modelled dust size distribution to a common particle bin spacing between 0.2 and 20 µm by assuming a power-law distribution between nearby model particle bins. After putting all the model simulations on the same footing in this manner, we thus represent the ensemble of the model dust size distribution with the mean, standard deviation, and range (minimum–maximum value) as a function of particle sizes, horizontal locations, heights, and seasons. Where necessary, the 95 % confidence interval of the model ensemble is estimated as 1.96 times the standard error (e.g. Altman and Bland, 2005). We also perform a similar aggregation and interpolation procedure on the modelled dust aerosol optical depth and column-integrated atmospheric dust loading, which are used to calculate the column-integrated dust mass extinction efficiency (MEE) for each model and thus for the model ensemble.

2.2 Reanalysis dust aerosol optical depth

We obtain the dust aerosol optical depth from four reanalysis products to constrain the atmospheric dust loading for DustCOMM (see Sect. 2.3.3). These four reanalysis products are the Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2; Gelaro et al., 2017); the Navy Aerosol Analysis and Prediction System (NAAPS; Lynch et al., 2016); the Japanese Reanalysis for Aerosol (JRAero; Yumimoto et al., 2017); and the Copernicus Atmosphere Monitoring Service (CAMS) interim Reanalysis (CAMSiRA; Flemming et al., 2017). While the description of each reanalysis product can be found in the Supplement, we give a general overview in this section.

A key advantage of these reanalysis products is that they assimilate data from several observing systems and thus provide a complete spatial and temporal coverage of atmospheric composition that captures its variabilities and trends (Buchard et al., 2017). Most of these four reanalysis products assimilate similar satellite and ground-based observations of AOD, which includes data from at least one or all of the following observing systems: the Terra and Aqua satellites of the MODerate resolution Imaging SpectroRadiometer (MODIS), the Advanced Very High Resolution Radiometer (AVHRR), the Multi-angle Imaging SpectroRadiometer (MISR), as well as ground-based observation of AOD from several Aerosol Robotic Network (AERONET) stations (Lynch et al., 2016; Flemming et al., 2017; Gelaro et al., 2017; Yumimoto et al., 2017). In addition, some reanalysis products also assimilate other aerosol constituents and reactive gases, like carbon monoxide and ozone observations from the Measurements Of Pollution In The Troposphere (MOPITT) instrument on the Terra Satellite, Solar Backscatter Ultraviolet (SBUV/2) instruments (from various National Oceanic and Atmospheric Administration (NOAA) platforms), and Microwave Limb Sounder (MLS) ozone profiles (e.g. Flemming et al., 2017). These observations are mostly bias-corrected before they are assimilated through radiatively coupled aerosol models and are used to constrain the different species that constitute the aerosol particles in the atmosphere.

Although the total AOD is constrained, errors in each reanalysis model's treatment of emission, transport, and deposition of mineral dust introduce uncertainties. Dust emission and deposition in the assimilation procedure are either modelled or sometimes constrained by observations. For example, the dust emission for NAAPS is constrained by using a regional source tuning that is, in turn, constrained by space-based and ground-based AOD observations (Lynch et al., 2016). Other reanalysis products use dust emissions that are parameterized and model-dependent (e.g. Yumimoto et al., 2017). In general, wet deposition is partially constrained by the assimilated global satellite-based precipitation information, such as from the NOAA Climate Prediction Center MORPHing technique data (CMORPH). Dry deposition is still mostly model dependent, but may also be adjusted based on assimilated AOD. For all the reanalysis products, aerosol transport in the atmosphere is constrained by the assimilation of several meteorological observations of winds and temperature. Hence, in order to constrain the dust AOD, the assimilation procedure takes advantage of the best features in both the observations and model simulations.

Similarly to our treatment of the model simulations described in Sect. 2.1 above, we use annual and seasonal climatologies of dust AOD obtained from monthly averages of the reanalysis products. We use the reanalysis dust AOD from 2004 to 2008 for each reanalysis product except for JRAero, for which we use 2011–2015. In order to combine the different reanalysis dust AOD products, we interpolate each product to approximately 2.5^∘ by 2.0^∘ spatial resolution and estimate the ensemble mean and standard error over each location (see Sect. 2.3.3).

2.3 Constraining DustCOMM products

Our aim is to create a new product – DustCOMM – that constrains the spatial variability of three major properties of atmospheric dust which determine many of its impacts on the Earth system, namely (1) the atmospheric dust size distribution, (2) the dust mass extinction efficiency, and (3) the column-integrated atmospheric dust loading. We do so by combining observational, experimental and theoretical constraints on dust properties and abundance with global model simulations of the size-resolved spatially varying dust concentration (Fig. 1). After we present a general overview of the methodology here, we describe the details of the methodology and the calculation of the associated uncertainty estimates in the following sub-sections.

We obtain the first constrained product in our dust climatology, the dust size distribution, by bias-correcting the six global model simulations (see Sect. 2.3.1; left panel of Fig. 1). Specifically, we bias-correct these model simulations using the constraint on the globally averaged dust size distribution from Kok et al. (2017), which was obtained from measurements of the emitted dust size distribution and model simulations of the globally averaged dust lifetime. Model simulations of the size-resolved dust lifetimes were used because this cannot be readily constrained with observations or measurements. Similarly, we use the constraints on the globally averaged size distribution from Kok et al. (2017) to correct modelled size distributions because dust size distribution measurements are insufficient to constrain the dust size distribution for every location. After correcting the model simulations of the dust size distribution, we combine them into a single multi-model constraint on the 3-D dust size distribution. To do this, we estimate the sub-bin distributions by fitting the dust size distribution after the bias correction with a generalized analytical function based on brittle fragmentation theory (Kok, 2011a). We then use the resulting distributions from the multiple models to obtain a constraint on the atmospheric dust size distribution, for each horizontal location and height level.

We use these constrained size distributions to obtain our second product, namely the size-integrated 3-D dust mass extinction efficiency (Sect. 2.3.2; middle panel in Fig. 1). Specifically, we combine the constrained 3-D dust size distribution with the constraint on the size-resolved globally averaged single-particle dust extinction efficiency at 550 nm obtained from Kok et al. (2017). This size-resolved single-particle dust extinction efficiency leverages measurements of the dust index of refraction and also accounts for the non-spherical shape of dust particles. As we did for the size distribution, we use the globally averaged dust extinction efficiency here because measurements of dust shapes and index of refraction are currently insufficient to constrain this for every location. As with the size distribution, we also constrain the mass extinction efficiency over each horizontal location and height level.

We obtain our third product – the column-integrated atmospheric dust loading – by combining the constraint on dust mass extinction efficiency with dust aerosol optical depth from multiple reanalysis products (Sect. 2.3.3; right panel in Fig. 1). Using four state-of-the-art reanalysis products (see Sect. 2.2), we calculate the ensemble average of dust aerosol optical depth, accounting for systematic and random errors. We propagate the errors in the dust mass extinction efficiency and dust aerosol optical depth to obtain the mean and the uncertainty of the column-integrated atmospheric dust loading over each horizontal location.

We estimate all DustCOMM products at a horizontal resolution of 2.5^∘ × 1.9^∘ with 35 levels, that is, up to 100 hPa.

2.3.1 Constraining the 3-D atmospheric dust size distribution

We constrain the spatially varying atmospheric dust size distributions by combining constraints on the globally averaged dust size distribution with an ensemble of simulations of the 3-D spatial variability of the dust size distribution (Fig. 1). We obtain the globally averaged atmospheric size distribution, ${\left[\frac{\mathrm{d}\widehat{V}\left(D\right)}{\mathrm{d}D}\right]}_{\mathrm{g}}$, from Kok et al. (2017; see their Fig. 2a). This globally averaged size distribution was obtained by combining constraints on the size distribution of emitted dust particles with simulations of the size-resolved dust lifetime. While details can be found in Kok et al. (2017), a summary of their globally averaged size distribution is given here as

$$\begin{array}{}\text{(1)}& {\left[{\displaystyle \frac{\mathrm{d}\widehat{V}\left(D\right)}{\mathrm{d}D}}\right]}_{\mathrm{g}}={\left[{\displaystyle \frac{\mathrm{d}{\widehat{V}}_{\mathrm{emit}}\left(D\right)}{\mathrm{d}D}}\right]}_{\mathrm{g}}\cdot {\left[{\displaystyle \frac{\stackrel{\mathrm{̃}}{T}\left(D\right)}{\stackrel{\mathrm{‾}}{\stackrel{\mathrm{̃}}{T}}}}\right]}_{\mathrm{g}},\end{array}$$

where the long-square parentheses []_g indicate quantities that are globally averaged, quantities with $\widehat{}$ accents are partially constrained by observations, and quantities with $\stackrel{\mathrm{̃}}{}$ accents are obtained from model simulations. As reported in Kok et al. (2017; hereafter referred to as K17), the constrained globally averaged size distribution of emitted dust particles, ${\left[\frac{\mathrm{d}{\widehat{V}}_{\mathrm{emit}}\left(D\right)}{\mathrm{d}D}\right]}_{\mathrm{g}}$, is based on an analysis of different measurements of the emitted dust size distribution, while the size-resolved globally averaged dust lifetime, ${\left[\stackrel{\mathrm{̃}}{T}\left(D\right)\right]}_{\mathrm{g}}$, is based on an ensemble of global model simulations; ${\left[\stackrel{\mathrm{‾}}{\stackrel{\mathrm{̃}}{T}}\right]}_{\mathrm{g}}$ is the mass-weighted mean of ${\left[\stackrel{\mathrm{̃}}{T}\left(D\right)\right]}_{\mathrm{g}}$. The constrained K17 globally averaged size distribution is normalized such that ${\int }_{\mathrm{0}}^{D_{\max }}{\left[\frac{\mathrm{d}\widehat{V}\left(D\right)}{\mathrm{d}D}\right]}_{\mathrm{g}}\mathrm{d}D=\mathrm{1}$, where D_max represents the maximum geometric diameter above which the contribution to extinction is negligible (D_max=20 µm; see Sect. 2.3.1).

We use this constrained K17 globally averaged atmospheric dust size distribution (Eq. 1) to bias-correct our spatially varying model simulations of the annually averaged dust size distribution. This is necessary because models generally underestimate coarse dust particles, largely because they assume too much fine dust in the emitted dust size distribution (Kok, 2011a). For each model, we force the simulated globally averaged dust size distribution to match the K17 constraint on the globally averaged size distribution (see Supplement Fig. S1), such that

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\widehat{f}}_{k,i}\left(x,y,z,D_{k,i}\right)& ={\stackrel{\mathrm{̃}}{f}}_{k,i}\left(x,y,z,D_{k,i}\right)\cdot {\mathit{\alpha }}_{k,i},\\ \text{where}{\mathit{\alpha }}_{k,i}& ={\displaystyle \frac{{\int }_{D_{k,i-}}^{D_{k,i+}}{\left[\frac{\mathrm{d}\widehat{V}\left(D\right)}{\mathrm{d}D}\right]}_{\mathrm{g}}\mathrm{d}D}{{\left[{\stackrel{\mathrm{‾}}{\stackrel{\mathrm{̃}}{f}}}_{k,i}\left(D_{k,i}\right)\right]}_{\mathrm{g}}}}.\end{array}\end{array}$$

The annually averaged 3-D distribution of the dust size distribution for each particle bin i simulated by model k is ${\stackrel{\mathrm{̃}}{f}}_{k,i}\left(x,y,z,D_{k,i}\right)$, and the corresponding simulated globally averaged dust mass fraction is ${\left[{\stackrel{\mathrm{‾}}{\stackrel{\mathrm{̃}}{f}}}_{k,i}\left(D_{k,i}\right)\right]}_{\mathrm{g}}$; x is the dimension for longitude, y is for latitude, and z is for height. Further, the numerator, ${\int }_{D_{k,i-}}^{D_{k,i+}}{\left[\frac{\mathrm{d}\widehat{V}\left(D\right)}{\mathrm{d}D}\right]}_{\mathrm{g}}\mathrm{d}D$, is the constraint obtained from Kok et al. (2017), while $D_{k,i-}$ and $D_{k,i+}$, respectively, denote the lower and upper geometric diameter limits of particle bin i of model k, and $i=\mathrm{1},\mathrm{2},\mathrm{\dots },N_k$, with N_k as the total number of dust particle bins for a given model simulation k. In Eq. (2), we multiply each simulated dust size distribution ${\stackrel{\mathrm{̃}}{f}}_{k,i}\left(x,y,z,D_{k,i}\right)$ by a correction factor α_k,i that is estimated as the ratio of the fractional contributions of the K17 globally averaged constraint to that obtained from the same model. This correction is done for each bin i of each model k, defined between $D_{k,i-}$ and $D_{k,i+}$. The resulting corrected spatially varying dust size distribution, ${\widehat{f}}_{k,i}\left(x,y,z,D_{k,i}\right)$, is normalized such that the discrete sum over each location and height equals unity, that is, $\sum _{i={\mathrm{1}}_k}^{N_k}{\widehat{f}}_{k,i}\left(x,y,z,D_{k,i}\right)=\mathrm{1}$.

Each model simulation in the ensemble has a particle size range and spacing that differ from other models (see Table 1 and Sect. 2.1 for details). In order to combine the corrected size distributions from the different models into a single estimate and to quantify the uncertainty across the different models, each corrected size distribution must be in a consistent size range and spacing with other models. We therefore process the corrected size distributions over a given location as follows: (1) we correct and scale each model's lower and upper diameter limits to the common diameter range of 0.2–20 µm (see Sect. 2.3.1); and (2) we estimate the sub-bin distribution for each model's bias-corrected size distribution by fitting a generalized analytical function, extending the Kok et al. (2017) theoretical expression of dust size distribution to the 3-D dataset (see Sect. 2.3.1).

Correcting model simulations to a common diameter range

For all simulations in the model ensemble, we set the lower and upper diameter limits to common limits defined by D_min=0.2 µm and D_max=20 µm, respectively. The lower diameter limit (D_min) is based on the lowest common diameter included in all the model simulations used in our analysis (Table 1). In addition, possible contaminations by other aerosol species are significantly more likely below 0.2 µm in measurements of dust aerosol particles (e.g. Dubovik et al., 2000). For these reasons, we set the lower diameter limit to D_min=0.2 µm, consistent with previous studies (e.g. Mahowald et al., 2014; Kok et al., 2017). Further, we set the upper diameter limit to D_max=20 µm, because most global models generally do not incorporate dust particles beyond 20 µm and also because the observational constraints on the size distribution from Kok et al. (2017) are limited to this maximum diameter. Although advances in airborne observations in recent years have led to measurements of larger dust particles with D>D_max in the atmosphere, which has shown that the contributions of D>20 µm to shortwave and longwave extinctions are non-negligible (e.g. Ryder et al., 2013b, 2019; Weinzierl et al., 2009, 2017), there is still a scarcity of these measurements, such that an observational constraint on dust particles with D>D_max would be very uncertain (e.g. Mahowald et al., 2014).

To correct each model simulation to the common diameter range of [D_min, D_max], we first create a new particle bin for the lower and/or upper diameter limit, and then we use the K17 constraints on the globally averaged size distribution (Eq. 1) to estimate the equivalent fraction of dust mass in that bin. This dust mass fraction is estimated in a way that is consistent with the size distribution obtained earlier from Eq. (2). Specifically, for simulations with a lower diameter limit ($D_{k,{\mathrm{1}}_k-}$) less than D_min, we estimate the equivalent dust mass fraction for the bin between D_min and $D_{k,{\mathrm{1}}_k+}$ (where $D_{k,{\mathrm{1}}_k+}$ is the upper diameter limit of bin 1, such that $D_{k,{\mathrm{1}}_k+}>D_{\min }$) by scaling the mass in the nearest bin with a factor that depends on the globally averaged size distribution. For instance, the first particle bin of the CESM (Table 1) has a range of $\left[D_{k,{\mathrm{1}}_k-},D_{k,{\mathrm{1}}_k+}\right]=\mathrm{0.1}-\mathrm{1.0}\mathrm{\mu }\mathrm{m},$ such that we create a new particle bin defined by $\left[D_{\min },D_{k,{\mathrm{1}}_k+}\right]=\mathrm{0.2}-\mathrm{1.0}\mathrm{\mu }\mathrm{m}$ and estimate the equivalent dust mass fraction in that new bin. For all model simulations, we can denote this procedure mathematically as

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{\widehat{f}}_k& \left(x,y,z,\left[D_{\min },D_{k,{\mathrm{1}}_k+}\right]\right)=\\ & {\widehat{f}}_k\left(x,y,z,\left[D_{k,{\mathrm{1}}_k-},D_{k,{\mathrm{1}}_k+}\right]\right)\cdot {\mathit{\delta }}_{D_{\min }},\\ & \text{where}{\mathit{\delta }}_{D_{\min }}={\displaystyle \frac{{\int }_{D_{\min }}^{D_{k,{\mathrm{1}}_k+}}{\left[\frac{\mathrm{d}\widehat{V}\left(D\right)}{\mathrm{d}D}\right]}_{\mathrm{g}}\mathrm{d}D}{{\int }_{D_{k,{\mathrm{1}}_k-}}^{D_{k,{\mathrm{1}}_k+}}{\left[\frac{\mathrm{d}\widehat{V}\left(D\right)}{\mathrm{d}D}\right]}_{\mathrm{g}}\mathrm{d}D}},\end{array}\end{array}$$

The modelled dust size distribution is relatively invariant for fine particles because of the consistent emitted dust size distribution (Kok, 2011a, b) and because removal processes for fine dust (wet deposition) do not strongly depend on particle size (e.g. Zender et al., 2003). Therefore, we simply estimate ${\mathit{\delta }}_{D_{\min }}$ in Eq. (3) as the ratio between the fractional values of the K17 globally averaged size distribution in the desired new bin [$D_{\min },D_{k,{\mathrm{1}}_k+}$] and in the model's original bin [$D_{k,{\mathrm{1}}_k-}{D}_{k,{\mathrm{1}}_k+}$].

We also create a new bin with the upper diameter equal to D_max for model simulations with an upper diameter limit ($D_{k,N_k+}$) that differs from D_max. We do so by scaling the nearest bin by a factor (${\mathit{\delta }}_{D_{\max }}$) that also depends, in part, on the constrained K17 globally averaged size distribution. Because the main removal process for large dust particles (D>10 µm) is dry deposition, which depends strongly on particle size, the relative contribution to the size distribution of different particle bins of large particles has substantial spatial variability. To account for this, we use simulations of bins with $D>D_{k,N_k+}$ from other model simulations in order to estimate what model k would have predicted for a hypothetical $\left[D_{k,N_k+},D_{\max }\right]$ particle bin. That is,

$$\begin{array}{}\text{(4a)}& \begin{array}{rl}{\widehat{f}}_k& \left(x,y,z,\left[D_{k,N_k+},D_{\max }\right]\right)=\\ & {\widehat{f}}_k\left(x,y,z,\left[D_{k,N_k-},D_{k,N_k+}\right]\right)\cdot {\mathit{\delta }}_{D_{\max }}\left(x,y,z\right),\\ & \text{where}{\mathit{\delta }}_{D_{\max }}\left(xyz\right)={\displaystyle \frac{{\int }_{D_{k,N_k+}}^{D_{\max }}{\left[\frac{\mathrm{d}\widehat{V}\left(D\right)}{\mathrm{d}D}\right]}_{\mathrm{g}}\mathrm{d}D}{{\int }_{D_{k,N_k-}}^{D_{k,N_k+}}{\left[\frac{\mathrm{d}\widehat{V}\left(D\right)}{\mathrm{d}D}\right]}_{\mathrm{g}}\mathrm{d}D}}\cdot {\mathit{\beta }}_{\mathrm{r}}\left(xyz\right).\end{array}\end{array}$$

The factor β_r thus quantifies the ratio of the mass fractions between the model's largest particle bin $\left(\left[D_{k,N_k-},D_{k,N_k+}\right]\right)$ and the newly created particle bin to extend the simulation to D_max=20 µm $\left(\left[D_{k,N_k+},D_{\max }\right]\right)$, as estimated from the GISS and ARPEGE-Climat simulations, which have particle bins extending to D_max (Table 1). We denote these latter model simulations with a subscript r for the purpose of clarity and to separate them from the model simulation that is being adjusted to the [D_min, D_max] size range, which is denoted by a subscript k in Eq. (4a) above. We thus estimate β_r as

$$\begin{array}{}\text{(4b)}& \begin{array}{rl}{\mathit{\beta }}_{\mathrm{r}}\left(x,y,z\right)=& {\displaystyle \frac{{\widehat{f}}_r\left(x,y,z,\left[D_{r,N_r-},D_{r,N_r+}\right]\right)}{{\int }_{D_{r,N_r-}}^{D_{r,N_r+}}{\left[\frac{\mathrm{d}\widehat{V}\left(D\right)}{\mathrm{d}D}\right]}_{\mathrm{g}}\mathrm{d}D}}/\\ & {\displaystyle \frac{{\widehat{f}}_r\left(x,y,z,\left[D_{r,j_r-},D_{r,j_r+}\right]\right)}{{\int }_{D_{r,j_r-}}^{D_{r,j_r+}}{\left[\frac{\mathrm{d}\widehat{V}\left(D\right)}{\mathrm{d}D}\right]}_{\mathrm{g}}\mathrm{d}D}},\end{array}\end{array}$$

where $\left[D_{r,N_r-},D_{r,N_r+}\right]$ is the bin in model r with dust mass that overlaps in size with the new bin $\left[D_{k,N_k+},D_{\max }\right]$ we want to estimate for model k; and $\left[D_{r,j_r-},D_{r,j_r+}\right]$ is the bin that similarly overlaps with $\left[D_{k,N_k-},D_{k,N_k+}\right]$. To account for the bin–range mismatch between the model simulation that resolved dust up to D_max (with subscript r) and the model simulation being adjusted to the dust size range up to D_max (with subscript k), we normalize each bin mass fraction by its contribution to the constrained K17 globally averaged size distribution. For cases where model r is the same as model k (i.e. for GISS and ARPEGE-Climat), β_r reduces to 1 everywhere. It should be noted that the correction of Eq. (4) takes into account the potential difference in the dust deposition between models k and r by considering the differences in the spatial variability of dust loading between similar bins of $\left[D_{k,N_k-},D_{k,N_k+}\right]$ and $\left[D_{r,j_r-},D_{r,j_r+}\right]$. After the dust mass fractions are corrected, they are re-normalized such that the discrete sum between D_min and D_max equals unity over each location and height.

This procedure described above (Eqs. 3 and 4) can be used to correct either the original modelled dust size distribution (Sect. 2.1) or the bias-corrected modelled dust size distribution of Eq. (2).

Estimating the sub-bin distribution of the dust size distribution

After setting the corrected dust size distribution from each model to a common diameter range, [D_min, D_max], we next estimate the sub-bin distribution in order to combine estimates from different models into one dust size distribution product. To do this, we fit a generalized theoretical function of the dust size distribution to the estimated bias-corrected dust size distribution from each model over each location and height level (Eq. 2). Although fitting log-normal modes are appropriate for several other aerosol species, Mahowald et al. (2014) highlighted that dust size distributions are usually not log-normal and are thus better characterized by a generalized function based on mechanistic understanding of dust emission and deposition processes. Therefore, we describe in this section the generalized function and the fitting procedure used to constrain the dust size distributions.

We define the generalized function for the atmospheric size distribution by considering the theoretical expressions that characterize the processes affecting the dust size distribution. The degree of the impact of any of these processes on the dust size distribution will depend on the location. For example, the impact of emission processes on the shape of the dust size distribution is expected to be large close to major dust source regions but less farther from source regions. Furthermore, farther from dust source regions, deposition processes are expected to have more impact on the size distribution. We therefore assume that the atmospheric size distribution over any location is proportional to the dust size distribution at emission, $\frac{\mathrm{d}{V}_{\mathrm{emit}}\left(D\right)}{\mathrm{d}D}$, the size-resolved dust lifetime in the atmosphere, T(D), and any other changes to the dust particle size distribution during transport, A(D) (e.g. Weinzierl et al., 2009; Schladitz et al., 2011; Kok et al., 2017). That is,

$$\begin{array}{}\text{(5a)}& {\displaystyle \frac{\mathrm{d}{V}_{\mathrm{atm}}\left(D\right)}{\mathrm{d}D}}\propto {\displaystyle \frac{\mathrm{d}{V}_{\mathrm{emit}}\left(D\right)}{\mathrm{d}D}}\cdot T\left(D\right)\cdot A\left(D\right).\end{array}$$

For the dust size distribution at emission, Kok (2011) suggested that $\frac{\mathrm{d}{V}_{\mathrm{emit}}\left(D\right)}{\mathrm{d}D}$ can be represented by a simple theoretical expression based on brittle fragmentation theory, which shows good agreement with measurements (e.g. Mahowald et al., 2014; Rosenberg et al., 2014). To better represent the variability in dust emission affecting the emitted size distribution in the different simulations, here we generalize this expression such that

$$\begin{array}{}\text{(5b)}& {\displaystyle \frac{\mathrm{d}{V}_{\mathrm{emit}}\left(D\right)}{\mathrm{d}D}}={\displaystyle \frac{\mathrm{1}}{C_{\mathrm{v}}}}\cdot \left[\mathrm{1}+\mathrm{erf}\left({\displaystyle \frac{\ln \left(\frac{D}{D_{\mathrm{v}}}\right)}{\sqrt{\mathrm{2}}\ln \left({\mathit{\sigma }}_{\mathrm{v}}\right)}}\right)\right]\cdot e^{\left[-{\left(\frac{D}{\mathit{\omega }}\right)}^{\mathit{\alpha }}\right]},\end{array}$$

where D_v and σ_v are, respectively, the geometric median diameter by volume and the geometric standard deviation of a typical desert soil, ω denotes the propagation distance of main cracks in dust aggregates during fragmentation, α is a tunable parameter primarily affecting the large dust particles, and C_v is a normalization constant.

The second term in our generalized dust size distribution describes the size-resolved dust lifetime, which global model results compiled in Kok et al. (2017) suggest can analytically be approximated as an exponential function of particle diameter, such that

$$\begin{array}{}\text{(5c)}& T\left(D\right)\cong T_{\mathrm{0}}\cdot e^{-\left(\frac{D}{\mathit{\kappa }}\right)},\end{array}$$

where T₀ is a constant associated with the lifetime for vanishingly small dust particles, which is determined by depositional processes, and κ is a constant that scales the exponential decay of the dust lifetime with particle size. This exponential decay of dust lifetime with size is caused by the increase in the gravitational settling speed with particle size (e.g. van der Does et al., 2016, 2018).

Finally, we account for other changes to the dust size distribution during transport by assuming that such changes are likely described by power-law distribution (e.g. Seinfeld and Pandis, 2016). Maring et al. (2003) highlighted that between emission and deposition, changes in dust size distribution cannot be accounted for by simple preferential removal of dust particles by gravitational settling. Since such changes in the dust size distribution are difficult to account for, we represent them with a parameter that can affect the entire size range. In addition, T(D) and $\frac{\mathrm{d}{V}_{\mathrm{emit}}\left(D\right)}{\mathrm{d}D}$ represent expressions that describe the globally averaged size distributions, and applying them to a specific location requires an additional parameter that captures the loss rate as a function of location. To represent all other changes to the dust size distribution between emission and deposition, we thus define

$$\begin{array}{}\text{(5d)}& A\left(D\right)\propto D^b.\end{array}$$

Combining Eqs. (5b) and (5d), we obtain

$$\begin{array}{}\text{(5e)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}{V}_{\mathrm{atm}}\left(D\right)}{\mathrm{d}D}}& ={\displaystyle \frac{\mathrm{1}}{C_{\mathrm{v}}}}\left[\mathrm{1}+\mathrm{erf}\left({\displaystyle \frac{\ln \left(\frac{D}{D_s}\right)}{\sqrt{\mathrm{2}}\ln \left({\mathit{\sigma }}_s\right)}}\right)\right]\\ & e^{\left[-{\left(\frac{D}{\mathit{\omega }}\right)}^{\mathit{\alpha }}\right]}\cdot T_{\mathrm{0}}{e}^{-\left(\frac{D}{\mathit{\kappa }}\right)}\cdot D^b.\end{array}\end{array}$$

We combine the two exponential terms in Eq. (5e) in order to reduce the number of fitting parameters. It is worth noting that both parameters κ and ω are sensitive to the larger particles, as they remain highly uncertain and poorly constrained by observation (e.g. Mahowald et al., 2014). The parameter ω depends on the soil moisture, mineralogy, and other processes (e.g. Mahowald et al., 2014; Rosenberg et al., 2014; Kok et al., 2017), while the parameter κ depends on the dust wet and dry deposition rates, as the dust particles are transported away from the source (e.g. Han and Zender, 2010; van der Does et al., 2016). To combine them, we define Λ to account for the uncertainty in the atmospheric large-size dust particles over every location. The generalized theoretical function for atmospheric size distribution therefore becomes

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}{V}_{\mathrm{atm}}\left(D\right)}{\mathrm{d}D}}& ={\displaystyle \frac{\mathrm{1}}{C_{\mathrm{v}}^{\ast }}}\cdot \left[\mathrm{1}+\mathrm{erf}\left({\displaystyle \frac{\ln \left(\frac{D}{D_s}\right)}{\sqrt{\mathrm{2}}\ln \left({\mathit{\sigma }}_s\right)}}\right)\right]\\ & \cdot e^{\left[-{\left(\frac{D}{\mathrm{\Lambda }}\right)}^{\mathit{\alpha }}\right]}\cdot D^b,\end{array}\end{array}$$

where $C_{\mathrm{v}}^{\ast }$ is a new normalization constant that is obtained from requiring that the integral over Eq. (6) from D_min to D_max yields unity.

To determine the parameters in Eq. (6) for each height, horizontal location, season, and model simulation, we fit the generalized size distribution of Eq. (6) to the corresponding corrected dust size distribution from Eq. (2) above. To do this, we minimize the chi-squared (${\mathit{\chi }}_k^{\mathrm{2}}$) value for each height, for each location, and for each model k, such that

$$\begin{array}{}\text{(7)}& {\mathit{\chi }}_k^{\mathrm{2}}=\sum _i^{N_k}{\left[\log \left(\underset{D_{k,i-}}{\overset{D_{k,i+}}{\int }}{\displaystyle \frac{\mathrm{d}{V}_{\mathrm{atm}}}{\mathrm{d}D}}\mathrm{d}D\right)-\log \left({\widehat{f}}_{k,i}\right)\right]}^{\mathrm{2}}.\end{array}$$

In each case, we estimate the constrained dust size distribution, $\frac{\mathrm{d}{\widehat{V}}_{\mathrm{atm}}}{\mathrm{d}D}(\mathrm{x},\mathrm{y},\mathrm{z})$, based on the parameters we determine from Eq. (7). In order to restrict the fitted function to physically realistic dust size distributions, we set the following bounds for the five parameters of Eq. (6): D_s=0.25 to 6.0 µm; σ_s=1.6 to 4.0; $\mathrm{\Lambda }=\mathrm{1}\mathrm{to}\mathrm{30}\mathrm{\mu }\mathrm{m};\mathit{\alpha }=\mathrm{1}\mathrm{to}\mathrm{6}$; and $b=-\mathrm{10}\mathrm{to}\mathrm{4}$, consistent with previous studies (e.g. Kok, 2011a; Kok et al., 2017; Rosenberg et al., 2014). The probability distribution of these parameters for all heights, horizontal locations, and model simulations of the annually averaged dust size distribution is shown in Supplement Fig. S2. Finally, we note here that although our generalized theoretical function of Eq. (6) builds on the brittle fragmentation theory of Kok (2011), it adds analytical expressions of dust deposition and dust changes during transport that allow us to better fit different shapes of dust size distribution over different locations.

2.3.2 Constraining the 3-D dust mass extinction efficiency

After obtaining the constrained atmospheric dust size distributions (Sect. 2.3.1 above), we combine it with constraints on size-resolved single-particle extinction efficiency at 550 nm to obtain constraints on the 3-D dust mass extinction efficiency (${\widehat{\mathit{ϵ}}}_{\mathit{\tau }}-{\mathrm{m}}^{\mathrm{2}}{\mathrm{g}}^{-\mathrm{1}}$). That is (see also Kok et al., 2017),

$$\begin{array}{}\text{(8)}& {\widehat{\mathit{ϵ}}}_{\mathit{\tau }}\left(x,y,z\right)=\underset{D_{\min }}{\overset{D_{\max }}{\int }}{\displaystyle \frac{\mathrm{d}{\widehat{V}}_{\mathrm{atm}}\left(x,y,z,D\right)}{\mathrm{d}D}}{\displaystyle \frac{\mathrm{3}}{\mathrm{2}{\mathit{\rho }}_{\mathrm{d}}D}}{\widehat{Q}}_{\mathrm{ext}}\left(D\right)\mathrm{d}D,\end{array}$$

where $\frac{\mathrm{d}{\widehat{V}}_{\mathrm{atm}}\left(x,y,z,D\right)}{\mathrm{d}D}$ is the constrained atmospheric dust size distribution at a given location and height with sub-bin distribution (Eq. 6); ${\mathit{\rho }}_{\mathrm{d}}=\mathrm{2.5}\pm \mathrm{0.2}$ g cm⁻³ is the globally averaged density of dust aerosols (Fratini et al., 2007; Reid et al., 2008; Kaaden et al., 2009; Sow et al., 2009; Kok et al., 2017), and its error range is expected to account for the spatial and temporal variability of dust density (e.g. Tegen and Fung, 1994; Li et al., 2008); and ${\widehat{Q}}_{\mathrm{ext}}\left(D\right)$ is the globally averaged size-resolved single-particle extinction efficiency at 550 nm wavelength, with the extinction cross section normalized by πD²∕4 – the projected area of a sphere with diameter D.

We obtain ${\widehat{Q}}_{\mathrm{ext}}$ from Kok et al. (2017), which constrained the dust extinction efficiency by combining measurements of the dust index of refraction and probability distribution of dust particle shape with the single-scattering database of Meng et al. (2010). Specifically, Kok et al. (2017) estimated the globally averaged values of the real and imaginary dust indexes of refraction as $n=\mathrm{1.53}\pm \mathrm{0.03}$ and $\log \left(-k\right)=-\mathrm{2.5}\pm \mathrm{0.3}$ (Sokolik et al., 1993; Patterson et al., 1977; Dubovik et al., 2002; Kandler et al., 2009; Kim et al., 2011; Denjean et al., 2016), and both are assumed to be normally distributed. Dust particle shapes were represented by the dust aspect ratio – the ratio of the major and minor axes of an ellipsoid best fit to the irregularly shaped 2-D image of a dust particle – and the height-to-width ratio. Kandler et al. (2007) showed that the deviation of measured dust aspect ratios from a sphere can be approximated by a log-normal distribution, with typical values ranging from 1 – a perfect sphere – to about 3, and a median between ∼1.5 and 1.9. Based on aggregates of measurements (Okada et al., 2001; Reid et al., 2003a; Kandler et al., 2007; Chou et al., 2008; Kandler et al., 2009, 2011; Scheuvens et al., 2011; Scheuvens and Kandler, 2014), Kok et al. (2017) estimated the median and geometric standard deviation for the distribution of the dust aspect ratio as 1.7±0.3 and 0.6±0.2, respectively. Similarly, based on limited available measurements of the dust height-to-width ratio (Okada et al., 2001; Chou et al., 2008; Veghte and Freedman, 2014), Kok et al. (2017) used a mean value of 0.333 (see details in the supplementary document of Kok et al., 2017). By combining these constraints on the optical properties and shape of the ensemble of dust particles in Earth's atmosphere with the single-scattering database of Meng et al. (2010), Kok et al. (2017) obtained a constraint on the globally averaged size-resolved extinction efficiency ${\widehat{Q}}_{\mathrm{ext}}\left(D\right)$, which explicitly accounts for the enhancement of extinction by the asphericity of dust. Specifically, they found that accounting for dust asphericity enhances the extinction produced by a unit mass loading of dust by 29±5 % over the extinction calculated from Mie theory for spherical dust particles, which is used in most climate models.

We use this constrained K17 globally averaged ${\widehat{Q}}_{\mathrm{ext}}$ to constrain ${\widehat{\mathit{ϵ}}}_{\mathit{\tau }}$ (Eq. 8) for every location. We thus neglect any regional variation in ${\widehat{Q}}_{\mathrm{ext}}$ because measurements of dust shapes and index of refraction are currently insufficient to constrain ${\widehat{\mathit{ϵ}}}_{\mathit{\tau }}$ on a regional basis. In addition, since measurements of dust refractive index needed to constrain ${\widehat{\mathit{ϵ}}}_{\mathit{\tau }}$ at other wavelengths are also scarce, we limit our estimate here only to the 550 nm wavelength. We use 550 nm as the wavelength of choice because measurements to validate our estimate of ${\widehat{\mathit{ϵ}}}_{\mathit{\tau }}$ and the observational constraints to estimate the dust atmospheric loading are mostly available at mid-visible wavelength.

2.3.3 Constraining the 2-D atmospheric dust loading

We now combine the above-estimated mass extinction efficiency at 550 nm (Sect. 2.3.2) with dust aerosol optical depth at the same wavelength to constrain the atmospheric column dust loading ($\widehat{L}-\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}$) (Kaufman et al., 2005; Kok et al., 2017). Because our constraints on dust size distributions are normalized to unity, and also to ensure that our estimates of dust loading produce the same extinction as those from reanalysis dataset or satellite measurements, we use this approach to estimate the atmospheric dust loading, such that

$$\begin{array}{}\text{(9)}& \widehat{L}\left(x,y\right)={\displaystyle \frac{{\widehat{\mathit{\tau }}}_{\mathrm{d}}\left(x,y\right)}{{\widehat{\mathit{ϵ}}}_m\left(x,y\right)}},\end{array}$$

where ${\widehat{\mathit{ϵ}}}_m$ (unit: m² g⁻¹) is the vertically integrated 2-D mass extinction efficiency calculated from ${\widehat{\mathit{ϵ}}}_{\mathit{\tau }}$ and ${\widehat{\mathit{\tau }}}_{\mathrm{d}}\left(x,y\right)$ is obtained from an ensemble of reanalysis dust aerosol optical depth products.

The ensemble dust AOD climatology is obtained from the average of four different reanalysis products (MERRA-2, JRAero, NAAPS, and CAMSiRA; see Sect. 2.2 for details). This individual reanalysis dataset assimilates several satellite and ground-based measurements from multiple platforms, including MODIS (Terra and Aqua), AVHRR, and MISR satellites, as well as from the ground-based AERONET stations (Lynch et al., 2016; Mccarty et al., 2016; Flemming et al., 2017; Yumimoto et al., 2017). As such the assimilation procedure takes advantage of the best features in both the observations and model simulations, thus producing column-integrated dust AOD that is largely representative of what is observed, based on validation studies (e.g. Buchard et al., 2017).

Despite the advantage of assimilating observational datasets, estimating a realistic overall error in the dust AOD across the reanalysis datasets is difficult yet important. Here, we estimate the total error (σ_d) by considering both the systematic error (σ_sys) and the random error (σ_rand) inherent in the reanalysis-derived dust AOD. As such, we estimate the uncertainty in dust AOD as ${\mathit{\sigma }}_{\mathrm{d}}\left(x,y\right)=\sqrt{{\mathit{\sigma }}_{\mathrm{sys}}^{\mathrm{2}}\left(x,y\right)+{\mathit{\sigma }}_{\mathrm{rand}}^{\mathrm{2}}\left(x,y\right)}$. We define the σ_rand as the standard error between the four datasets, which represents that part of the total uncertainty that does not correlate across the four reanalysis dust AOD datasets. For instance, σ_rand may be associated with differences in the assimilating systems for the different reanalysis products. In contrast, the σ_sys is expected to correlate between the four datasets since most of the reanalysis datasets use similar observational datasets. Hence, we assume that the σ_sys will be proportional to the mean dust AOD, such that $\frac{{\mathit{\sigma }}_{\mathrm{sys}}\left(x,y\right)}{{\mathit{\tau }}_{\mathrm{d}}\left(x,y\right)}=C_{\mathrm{d}}$. We estimate the proportionality constant, C_d, by requiring that the relative error is the same as the relative error obtained from annually averaged climatology of dust AOD from Ridley et al. (2016), which leveraged observational datasets similar to those used for the reanalysis dataset but propagated many of the relevant uncertainties. From that, we deduce that $C_{\mathrm{d}}=\sqrt{\frac{{\mathit{\sigma }}_{\mathrm{o}}^{\mathrm{2}}}{{\mathit{\tau }}_{\mathrm{o}}^{\mathrm{2}}}-\frac{{\mathit{\sigma }}_{\mathrm{rand}}^{\mathrm{2}}}{{\mathit{\tau }}_{\mathrm{d}}^{\mathrm{2}}}}$, where τ_o and σ_o are the mean and error estimates of the observationally constrained dust AOD from Ridley et al. (2016), respectively. We estimate C_d=0.26 for annual climatology (between 2004 and 2008), averaged over regions that are constrained by Ridley et al. (2016) and account for about 95 % of the global dust AOD. Similarly, we estimate 0.31, 0.22, 0.24, and 0.28 for December–February, March–May, June–August, and September–November seasonal climatologies, respectively.

2.3.4 Quantifying the uncertainties in DustCOMM products

For each DustCOMM product above – the constrained dust size distribution, dust mass extinction efficiency, and dust atmospheric loading – we describe here how we estimate the most likely value and quantify the uncertainty over each location. Specifically, we use a non-parametric procedure based on the bootstrap method (Efron and Gong, 1983; Chernick, 2007). We use this method because the complexity of the equations (Eqs. 1–9) prevents a parametric quantification of error, and the bootstrap approach allows us to nonetheless propagate the uncertainty in the various physical variables used to estimate each product. Using this method, we further assume that the set of input variables in relevant equations above are independent and are represented by defined probability distributions. Thus, we estimate the probability distribution of the resulting products by randomly sampling (with replacement) the probability distribution of each of the input variables for a large number of times (n≈1500).

In practice, the procedure uses the following steps to determine the dust size distribution, mass extinction efficiency, and atmospheric loading, and their uncertainties.

• We randomly select a realization of the globally averaged size distribution from Kok et al. (2017), which in turn was obtained in that study by randomly selecting a realization of the emitted dust size distribution and the dust lifetime (Eq. 1).

• We use this randomly selected constrained K17 globally averaged size distribution to correct a randomly selected model simulation (Eq. 2).

• After this model simulation is corrected, we then scale the resulting 3-D dust size distribution between D_min and D_max following Eqs. (3) and (4).

• We thereafter estimate the constrained dust size distribution, $\frac{\mathrm{d}{\widehat{V}}_{\mathrm{atm}}}{\mathrm{d}D}\left(xyzD\right)$, and obtain the sub-bin distribution by fitting the generalized theoretical expression (Eq. 6) and minimizing the chi square over each location (Eq. 7).

• We randomly select a realization of the globally averaged size-resolved single-particle extinction efficiency, ${\widehat{Q}}_{\mathrm{ext}}\left(D\right)$, from Kok et al. (2017). This realization is also similarly estimated by randomly selecting from the distribution of the dust index of refraction and dust shape distribution parameters, as explained in Sect. 2.3.2.

• We then use the randomly selected ${\widehat{Q}}_{\mathrm{ext}}\left(D\right)$ and $\frac{\mathrm{d}{\widehat{V}}_{\mathrm{atm}}}{\mathrm{d}D}\left(xyzD\right)$ to estimate the dust mass extinction efficiency over each location, ${\widehat{\mathit{ϵ}}}_{\mathit{\tau }}\left(xyzD\right)$, following Eq. (8). This uses a randomly selected dust density value (ρ_d) from its assumed normal distribution.

• Similarly, assuming a normal distribution for the dust AOD, we randomly estimate the ${\widehat{\mathit{\tau }}}_{\mathrm{d}}\left(x,y\right)$ value within the range of its uncertainty, σ_d(x, y).

• We use this ${\widehat{\mathit{\tau }}}_{\mathrm{d}}\left(x,y\right)$ and the vertically integrated value of dust mass extinction efficiency, ${\widehat{\mathit{ϵ}}}_m\left(x,y\right)$, to estimate the atmospheric dust loading, $\widehat{L}\left(x,y\right)$, following Eq. (9).

• We repeated steps 1–8 for n=1500 times, thereby producing a probability distribution for $\frac{\mathrm{d}{\widehat{V}}_{\mathrm{atm}}}{\mathrm{d}D}\left(xyzD\right)$, ${\widehat{\mathit{ϵ}}}_{\mathit{\tau }}\left(xyzD\right)$, and $\widehat{L}\left(x,y\right)$ for each location and height. We report the mean, median, 1-sigma uncertainty range (68 % of the distribution) and the 95 % confidence interval (95 % CI) of those distributions (Adebiyi et al., 2019a).

The above procedure propagates various uncertainties in the estimation of each product. These include the measurement uncertainties and the uncertainties in model simulations. First, the measurement uncertainties are associated with the K17 globally averaged size distribution and the globally averaged extinction efficiency (Fig. 1), and these are propagated equally to every location. In addition, we estimated the correlated systematic error in the dust AOD (Sect. 2.3.3), associated with the assimilated observational dataset, and this is also propagated. Second, the uncertainty in model simulations is associated with the spread of the model dust size distribution, which is different for every location. This model uncertainty is, in turn, a result of many processes, such as dust emission, deposition, and transport processes in the models (Ginoux et al., 2001; Huneeus et al., 2011; Zhao et al., 2013). Our procedure constrains these model uncertainties (see Supplement Fig. S1) while retaining the spatial distribution of the model ensemble.

To quantify the size-resolved discrepancies in the DustCOMM size distribution, we quantify the bias with respect to independent measurements as follows (e.g. Lee et al., 2009):

$$\begin{array}{}\text{(10)}& {\mathit{\psi }}_i^{\mathrm{bias}}={\displaystyle \frac{\mathrm{1}}{N_m}}\sum _{m=\mathrm{1}}^{N_m}{\log }_{\mathrm{10}}\left({\displaystyle \frac{M_{i,m}^f}{O_{i,m}^f}}\right),\end{array}$$

where m sums over the N_m in situ measurements of the dust size distribution available in the literature (see Table 2), $O_{i,m}^f$ is the mth measurement of the mass fraction contained in measurement bin i, and $M_{i,m}^f$ is the corresponding mth DustCOMM dust mass fraction for the same diameter range as measured and collocated with the measurement – i.e. $M_{i,m}^f={\int }_{D_{i-}}^{D_{i+}}\frac{\mathrm{d}{\widehat{V}}_{\mathrm{atm}}}{\mathrm{d}D}\mathrm{d}D$. ${\mathit{\psi }}_i^{\mathrm{bias}}$ is the log-mean normalized bias and it represents the average number of orders of magnitude bias for each bin i.

Table 2 Overview of in situ measurements used to evaluate DustCOMM and model ensemble estimates of the dust size distribution and dust mass extinction efficiency (see Sect. 2.4 for details). The label PSD in the last column indicates that we take dust size distribution values from the study. MSE and MEE similarly indicate that we take dust mass scattering or extinction efficiency values.

We also estimate the performance of DustCOMM mass extinction efficiency by quantifying the reduced chi square (${\mathit{\chi }}_{\mathit{\epsilon }}^{\mathrm{2}}$) defined as the chi square per degree of freedom (e.g. Bevington et al., 1993):

$$\begin{array}{}\text{(11)}& {\mathit{\chi }}_{\mathit{\epsilon }}^{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{{\mathit{\upsilon }}_{\mathit{\epsilon }}}}\sum _{m=\mathrm{1}}^{N_m}\left({\displaystyle \frac{O_m^{\mathit{\epsilon }}-{\widehat{\mathit{ϵ}}}_{\mathit{\tau },m}}{{\mathit{\sigma }}_m^{\mathit{\epsilon }}}}\right),\end{array}$$

where $O_m^{\mathit{\epsilon }}$ is the mth measurement of the dust mass extinction efficiency with error defined as ${\mathit{\sigma }}_m^{\mathit{\epsilon }}$, ${\widehat{\mathit{ϵ}}}_{\mathit{\tau },m}$ is the corresponding mth DustCOMM dust extinction efficiency (${\widehat{\mathit{ϵ}}}_{\mathit{\tau }}$) collocated with the measurement, and υ_ε is the number of degrees of freedom given as N_m−1. A value of ${\mathit{\chi }}_{\mathit{\epsilon }}^{\mathrm{2}}\approx \mathrm{1}$ in Eq. (11) indicates there is agreement between DustCOMM and observations that is in accordance with the measurement errors, while ${\mathit{\chi }}_{\mathit{\epsilon }}^{\mathrm{2}}\gg \mathrm{1}$ indicates that DustCOMM estimates do not fully capture the observations (e.g. Andrae et al., 2010).

To facilitate comparison between DustCOMM and model evaluations, Eqs. (10) and (11) are also used to evaluate the performance and calculate the discrepancies between the measurements and the model ensemble.

2.3.5 DustCOMM at other timescales

While we describe above the procedure that constrains the annually averaged dust size distribution, dust mass extinction efficiency, and atmospheric dust loading, a similar procedure to that highlighted above can also be used to constrain the three products at any other timescale, such as at daily, monthly, or seasonal timescales. For this study, we only consider the seasonally averaged and annually averaged products.

First, to constrain the dust size distribution at any specific timescales, we correct an ensemble of model size distributions at that timescale in a way similar to Eq. (2) above. However, unlike Eq. (2) that uses the constrained globally averaged size distribution, here we use the constrained annually averaged dust size distribution over every location. That is,

$$\begin{array}{}\text{(12)}& \begin{array}{rl}{\widehat{f}}_{k,i}^{\mathrm{t}}\left(x,y,z,D_{k,i}\right)& ={\stackrel{\mathrm{̃}}{f}}_{k,i}^{\mathrm{t}}\left(x,y,z,D_{k,i}\right)\\ & \cdot {\displaystyle \frac{{\int }_{D_{k,i-}}^{D_{k,i+}}\frac{\mathrm{d}{\widehat{V}}_{\mathrm{atm}}}{\mathrm{d}D}\left(x,y,z,D\right)\mathrm{d}D}{{\stackrel{\mathrm{̃}}{f}}_{k,i}\left(x,y,z,D_{k,i}\right)}},\end{array}\end{array}$$

where $\frac{\mathrm{d}{\widehat{V}}_{\mathrm{atm}}}{\mathrm{d}D}$ is the DustCOMM annually averaged dust size distribution at a given 3-D location, obtained from the procedure described in Sect. 2.3.1, while $\stackrel{\mathrm{̃}}{f}$ and ${\stackrel{\mathrm{̃}}{f}}^{\mathrm{t}}$ are the annually averaged and specific time-averaged model simulations of the dust size distribution, respectively. Using an ensemble of model simulations, as we do above, the resulting corrected time-averaged dust size distributions, ${\widehat{f}}^{\mathrm{t}}$, are also taken through the steps highlighted in Sect. 2.3.4 to calculate the mean and the uncertainty of the constrained dust size distribution ($\frac{\mathrm{d}{V}_{\mathrm{atm}}^{\mathrm{t}}}{\mathrm{d}D})$ at that particular timescale.

Second, to constrain the dust mass extinction at any specific timescale (${\widehat{\mathit{ϵ}}}_{\mathit{\tau }}^{\mathrm{t}}$), we combine the constrained dust size distribution at that timescale, $\frac{\mathrm{d}{V}_{\mathrm{atm}}^{\mathrm{t}}}{\mathrm{d}D}$, with the globally averaged extinction efficiency, ${\widehat{Q}}_{\mathrm{ext}}$. This similarly follows Eq. (8) above. We note here that the uncertainty range of ${\widehat{Q}}_{\mathrm{ext}}$ also accommodates the location-dependent and time-dependent variability in the dust index of refraction and dust particle shape, consistent with previous studies (e.g. Dubovik et al., 2002). Hence, using ${\widehat{Q}}_{\mathrm{ext}}$ propagates the uncertainty in the measurements that determine the dust mass extinction efficiency estimate at that timescale. Finally, we constrain the dust loading at any specific timescale (${\widehat{L}}^{\mathrm{t}}$) using the constrained ${\widehat{\mathit{ϵ}}}_{\mathit{\tau }}^{\mathrm{t}}$ and dust AOD at that same timescale, similarly following Eq. (9).

2.4 Description of measurements used for evaluation

We use several types of published measurements to evaluate the dust size distribution and dust mass extinction efficiency from both DustCOMM and the model ensemble. We select 21 studies that measured dust properties – 14 of these reported dust size distributions and 11 of these reported dust mass extinction or scattering efficiencies (Table 2). These measurements were taken both near and far from dust-dominated regions (Table 2 and Supplement Fig. S3). While some measurements were taken close to (or over) some of the Northern Hemisphere deserts – such as the Saharan, Middle East, and Asian deserts – no measurements were taken close to the Southern Hemisphere deserts. Of the 21 studies, 12 obtained measurements near the Sahara, while one measurement each was taken near the Middle East (Sde Boker, Israel) and Asian (Qinghai Province, China) deserts. Other measurements represent dust properties at different distances of transport away from the dust sources.

Except for four measurements, most of the data are taken during airborne field campaigns that often occur over a wide geographical area, several altitude levels, and several days (Table 2). As such, studies often report measurements that represent the averages of the dust properties taken during the campaign. Details of the flight path, showing the locations where dust particles are encountered, are not always reported. To use these measurements, we therefore define a representative location and altitude for each measurement based on the area where the majority of dust was encountered. In addition, since the measurements often represent the average of several days and sometimes multiple months, we also compare them against seasonal averages of the DustCOMM and model ensemble estimates.

Below, we give a broad overview of the measurements of the dust size distribution and mass extinction efficiency, and further information on each study, including the instruments used, can be found in the Supplement.

2.4.1 Dust size distribution measurements

Dust size distribution measurements are taken using a variety of instruments with different sizing methodologies (e.g. Reid et al., 2003b). These instruments generally fall within the categories of sample collectors (e.g. D'Almeida, 1987; McConnell et al., 2008), cascade impactors (e.g. Chou et al., 2008; Kandler et al., 2009) and aerodynamic particle sizers (e.g. Otto et al., 2007), and optical particle counters or spectrometers (e.g. Chou et al., 2008; Clarke, 2004; Otto et al., 2007). The first category of instruments, sample collectors, are usually installed behind filters or thermal denuders to remove non-dust particles. The aerosol samples are then analysed using electron or light microscopy techniques, where they are counted and sized either manually or using an automated software. This type of measurement yields dust size distribution with respect to geometric diameters. For the second category of instruments, cascade impactors and particle sizers, aerosol particles are usually accelerated through a jet outlet and sometimes collected on a substrate. Using these instruments, the aerosols are sized based on the mass-to-drag characteristics of the particles. Dust particle sizes measured using these types of instruments are associated with the aerodynamic diameter. Finally, the optical particle counters generally determine particle sizes in optical diameters based on the amount of light they scatter. Another category is the imaging probe whereby the particle image is detected by a linear photodiode array providing a 2-D projection of the particle (Baumgardner et al., 2017; Ryder et al., 2018). For many of the studies we use here, these instruments are sometimes combined to verify the accuracy of the measurements (e.g. Ryder et al., 2013a). For all dust size distribution measurements, the studies that used aerodynamic or optical sizing instruments eventually report the measured size distribution in geometric diameters.

An important consideration is the elevation at which the dust size distributions are measured. With the exception of two studies (D'Almeida, 1987; Kandler et al., 2011) that took measurements at ground stations, most measurements were performed solely aboard aircrafts with in-cabin or wing-mounted instruments. Ground stations were equipped with stationary instruments to collect aerosol samples or stationary optical particle counters to measure size distributions directly. For aircraft measurements, size distributions are often measured during flight segments at constant altitude – also called horizontal legs. For the dust size distributions, our criteria for selection of studies are as follows: (1) the measured size range of the data should extend into the coarse dust (D>5 µm) size range; (2) the study should report the original in situ measurements, instead of (log-normal) fits to the actual measurements; and (3) each study's measurements should be taken with commonly used instrumentation in order to ensure some consistency with measurements taken by other studies.

Regardless of the instrument used, most dust size distribution measurements are subject to uncertainties associated with measurement type or presence of other aerosol species, such as biomass burning aerosols. The contamination by other aerosol species is common for fine-mode dust particles, especially dust particles less than ∼0.5 µm (e.g. Dubovik et al., 2000; Clarke, 2004), since to the instruments these aerosols are indistinguishable from dust particles of the same size. This causes a high bias in the fine mode of measured dust size distributions (e.g. Clarke, 2004). Another important measurement error arises from assumptions made about the non-sphericity of dust particles. For example, during the microscopy analysis, particle diameters are usually determined as the volume-equivalent geometric diameters based on 2-D images (Chou et al., 2008). Because of the asphericity of dust aerosols, this could introduce biases (e.g. Huang et al., 2020; Okada et al., 2001). Since dust particles have a small height-to-width ratio (Okada et al., 2001), the resulting size distribution may overestimate dust particle diameters. In the case of cascade impactors and particle sizers, unusual dust particle shapes and the possibility of particle bouncing off the substrate may lead to significant bias, especially for coarse-mode particles. For in-cabin measurements, studies have shown that the loss rate of coarse dust particles can be substantial due to the aircraft's instrument inlet, therefore leading to lower sampling rate and size bias (e.g. von der Weiden et al., 2009). For dust measurements that used optical particle counters, irregularly shaped dust particles are often assumed to be spherical in order to convert them to volume-equivalent geometric diameters, but light scattering between spherical and non-spherical particles is different. In addition, optical particle counters also make assumptions about the refractive index to derive the dust size distribution and are affected by the non-monotonic increase in the intensity of scattered light with particle size (e.g. Weinzierl et al., 2011; Ryder et al., 2018). Unlike the optical particle counters that require assumption regarding dust refractive index and shape to convert scattered light intensity to particle size, the imaging probes are not subject to these uncertainties (Baumgardner et al., 2017; Ryder et al., 2018). Nevertheless, these assumptions often lead to biases that many studies try to account for to various degrees (e.g. Ryder et al., 2013a, b, 2018).

2.4.2 Mass extinction efficiency measurements

In the literature, the term dust mass extinction efficiency is sometimes used interchangeably with the mass scattering efficiency (MSE; e.g. Hand and Malm, 2007). This is because, for a typical solar wavelength at 550 nm, dust particles scatter more radiation than they absorb for D≤10 µm. Despite the strong scattering by these particles, larger particles (D≥10 µm) often exhibit substantial absorption relative to scattering in the visible wavelength (e.g. Ryder et al., 2018). In order to put all the measurements on the same equal footing, we convert the reported dust MSE in some of these studies to dust MEE by using a measured scattering albedo value of 0.95±0.03 (e.g. Haywood, 2003; Clarke, 2004; Ryder et al., 2018).

MEEs that are reported in the literature are generally derived using two methods: regression and theoretical methods (e.g. Hand and Malm, 2007). The regression method calculates the dust MEE as the slope between the dust extinction coefficient (m⁻¹) and the dust mass concentration (g m⁻³). In this case, the dust samples are typically collected using filters, while aerosol extinction is measured using nephelometers. The difficulty, however, is that measured total aerosol extinction from the nephelometer may be influenced by several aerosol species other than dust particles. Some studies ignore the impact of other aerosol species and derive the dust MEE using the total aerosol extinction and the collected dust mass concentration (e.g. Li et al., 1996). Others take advantage of the linear relationship between the aerosol extinction and mass concentrations in order to separate the column MEE into constituents that correspond to each aerosol species, using a multivariate linear regression method (e.g. Andreae et al., 2002; Maring et al., 2003). Such calculations therefore require that all the aerosol species contributing to the extinction are included. With this in mind, the regression-derived MEE is therefore subject to several systematic and random errors, including instrument uncertainties (Hand and Malm, 2007).

The theoretical method calculates the dust MEE using the measured size distributions of dust mass or number concentration (Seinfeld and Pandis, 2016). This may take the form of calculating the dust MEE directly using the dust size distribution and the estimate of single-particle extinction efficiency or indirectly by first calculating the size-resolved dust extinction coefficient, using dust size distribution, and then combining the result with dust mass concentration. In either case, the dust density, shape, and index of refraction are needed. While assumptions of dust density and index of refraction are typically based on previously reported measurements, dust shapes are generally assumed to be spherical, which is contrary to observations (e.g. Okada et al., 2001; Kandler et al., 2007). This is a major disadvantage that may result in an underestimation of the derived dust extinction efficiency (e.g. Kok et al., 2017). Another source of error is associated with the instrument used to measure the aerosol size distribution, which may assume certain mixing properties of the observed aerosols. For mobility measurements (differential mobility analyser, DMA), optical measurements (optical particle counter, OPC), or aerodynamic measurements (aerodynamic particle sizer, APS), aerosols are often assumed to be internally mixed (e.g. Quinn et al., 2002; Clarke, 2004). In contrast, for an impactor, aerosols are often assumed to be externally mixed (e.g. Chiapello et al., 1999; Osborne et al., 2008).

Despite the differences between both methods used to derive dust MEE from observed quantities, previous studies have highlighted that they both produce similar values within measurement uncertainties (e.g. Maring et al., 2000; Quinn et al., 2004). In addition, for measurements where only the mean dust MEE/MSEs are reported, but not the uncertainty estimates, we estimate here in this study what the measured uncertainty estimate could be by assuming that its relative uncertainty (that is, the ratio of the presumed uncertainty to the reported mean) is proportional to the mean relative uncertainty that is calculated from other measurements. While this estimated uncertainty may likely not be representative of the specific field campaign from which the measurement was taken, they are likely representative of the seasonal values over the region.

3 Results

In this section, we present the DustCOMM products obtained using the methodology and data described above. We first present the dust particle size distribution (PSD; Sect. 3.1) and then the dust MEE (Sect. 3.2). In each case, we evaluate the DustCOMM and model ensemble products against available in situ measurements. We show that DustCOMM products generally reproduce observations better than model ensemble estimates. We then compare the spatial variability of the DustCOMM products against the model ensemble. In Sect. 3.3, we compare the atmospheric dust loading obtained from both DustCOMM and the model ensemble, and we examine the spatial distribution of the uncertainty in all DustCOMM products in Sect. 3.4.

3.1 Dust size distribution

3.1.1 Evaluation of DustCOMM against measurements

We evaluate DustCOMM and the model ensemble PSD against available in situ measurements taken during field campaigns (Figs. 2 and 3). We compare these location-based measurements against season-averaged DustCOMM and model ensemble estimates. The reason for using the seasonal averages is justified in Sect. 2.4 above. An additional justification for the comparison between the individual measurements and the season-averaged DustCOMM and model ensemble estimates is that the variability of the normalized dust PSD within each season is relatively small, especially for dust with D≤10 µm (e.g. McConnell et al., 2008; Mahowald et al., 2014). Furthermore, most of these measurements are campaign averages often over a variety of cases that could be representative of the season-averaged size distribution.

Figure 2 Comparison of normalized dust size distributions between published in situ measurements (blue and purple dots; see Table 2) and season-averaged DustCOMM (black lines) and model ensemble (red lines) estimates. The grey shading shows the 95 % confidence interval for the DustCOMM dust size distributions, whereas the pink shading shows the range of the model ensemble size distributions. The size distributions are normalized between 2.5 and 10 µm. The comparisons are made at the nearest model grid points to the representative location and height level of the measurements.

Figure 3 Same as Fig. 2 above but as a function of height, which increases from bottom to top. The measurements plotted in the left panels are from Otto et al. (2007) taken during the ACE-2 campaign (June/July 1997) in the vicinity of the Canary Islands; the measurements plotted in the middle panels are from Ryder et al. (2013a) taken during the Fennec project (June 2011) near the Canary Islands; and the measurements plotted in the right panels are from Ryder et al. (2018) taken during the AER-D campaign in August 2015 near Cabo Verde.

Model simulations of dust PSD generally show substantial errors when compared against measurements. In each of the 12 studies used in Fig. 2, the model ensemble overestimates the observed fine-mode particles (defined here as D≤2.5 µm) and underestimates the coarse-mode particles (defined here as D≥5 µm). In some of the cases, the overestimation extends above D=2.5 µm and the underestimation below D=5 µm. Nevertheless, these differences are apparent in all the comparisons and are consistent with previous studies indicating more coarse-mode dust particles are in the atmosphere than models account for (e.g. van der Does et al., 2016, 2018; Kok et al., 2017; Ryder et al., 2018).

In contrast, the DustCOMM dust PSD shows overall better agreement against measurement than the model ensemble (Fig. 2). This improved agreement includes a substantial reduction of the underestimation of coarse-mode dust as well as a reduction of the overestimation of some fine-mode particle sizes. Although DustCOMM better reproduces the measurements for D≥0.5 µm, it shows poorer agreement for D≤0.5 µm (e.g. Fig. 2e, h, i, j), underestimating the measurements by about 1 to 2 orders of magnitude. For example, during DARPO (Fig. 2e; Wagner et al., 2009) and BACEX (Fig. 2h; Jung et al., 2013), the differences between DustCOMM PSD and the measurements are about 2 orders of magnitude. The D≤0.5 µm size range is also the size range in which measurements of dust PSD are potentially contaminated by the presence of other aerosol species (see Sects. 2.3.1 and 4.1). In addition to the disagreement for D≤0.5 µm, there is also some disagreement for D≥10 µm (e.g. Fig. 2d, e, h), although for fewer cases. Overall, the DustCOMM dust PSDs significantly better represent the measurements in the $\mathrm{0.5}\le D\le \mathrm{20}\mathrm{\mu }\mathrm{m}$ size range than the model ensemble.

DustCOMM also shows better agreement than the model ensemble against measurements of the dust PSD as a function of altitude (Fig. 3). We highlight here measurements taken from three campaigns: (1) the ACE-2 campaign (June–July 1997) in the vicinity of the Canary Islands (Otto et al., 2007); (2) the Fennec project (June 2011) between the Canary Islands and Mauritania/Mali (Ryder et al., 2013a); and (3) the AER-D campaign in August 2015 near Cabo Verde (Ryder et al., 2018). All three cases show that a significant fraction of coarse-mode dust particles, including with D≥10 µm, is transported off the coast of northern Africa. We compare these measurements at selected altitudes of 2500 m (2700 m for ACE-2), 4000, 5500, and 6000 m (7000 m for ACE-2). Similar to Fig. 2 above, the DustCOMM dust PSD agrees better with the measurements than the model ensemble for these measurements at similar 2-D locations but at different altitudes. For dust particles with D≤0.5 µm, the DustCOMM size distributions also differ from the measurements by about an order of magnitude (similar to Fig. 2) for altitude at 2500 m. However, this difference increases to more than 2 orders of magnitude above ∼4000 m altitude.

In summary, the overall differences between the in situ measurements and DustCOMM are significantly smaller than the differences between the measurements and the model ensemble, especially for D≥0.5 µm. To quantify this, we report the log-mean bias in each bin following Eq. (10) and using all the measurements shown in Figs. 2 and 3. DustCOMM shows an overall reduction in the bias relative to the model ensemble, except for dust particles with D≤0.5 µm (Fig. 4). For D≤0.5 µm, the model shows an average (95 % CI) positive log-mean bias of 0.26 (−0.08 to +0.6), while DustCOMM shows an average negative log-mean bias of −0.92 (−1.18 to −0.73). In contrast, DustCOMM shows a remarkable reduction in the average log-mean bias in the $\mathrm{0.5}\le D\le \mathrm{10}\mathrm{\mu }\mathrm{m}$ size range; for instance, the bias for the 5–10 µm bin is ∼90 % less than it is for the model ensemble. DustCOMM also shows a substantially reduced bias in the $\mathrm{10}\le D\le \mathrm{20}\mathrm{\mu }\mathrm{m}$ size range, although the bias here remains substantially negative, indicating a persistent underestimation of these coarse particles. On average, DustCOMM reduces the log-mean bias for dust particles with D≥0.5 µm by about 46 % relative to the model ensemble.

Figure 4 Average log-mean bias between measurements and DustCOMM (grey) or model ensemble (pink) estimates of dust size distributions (shown in Figs. 2 and 3) for different particle bins. The vertical bars represent the 95 % confidence interval.

3.1.2 Global comparison between DustCOMM and the model ensemble

Considering that the DustCOMM dust PSD agrees better with in situ measurements than the model ensemble, we now compare the differences between DustCOMM and model ensemble PSDs. Specifically, we first compare the differences in the shape of the globally averaged dust size distribution between DustCOMM and the model ensemble (Sect. 3.1.2). Second, we examine the changes in the spatial variability of the DustCOMM and model dust mass fraction as a function of particle size range (Sect. 3.1.2).

Differences in dust size distribution

As we already concluded based on in situ measurements, climate models globally overestimate fine-mode dust particles (D≤2.5 µm) and underestimate coarse-mode dust particles (D≥5 µm), relative to globally averaged DustCOMM dust PSD (compare black and coloured lines in Fig. 5a). On average, simulations in our model ensemble overestimate the dust mass fraction of the fine mode by ∼14 % and underestimate that of the coarse mode by ∼15 %. The degree of this deviation from DustCOMM depends on the model and can be as much as 50 % in the fine mode or 37 % in the coarse mode.

Figure 5 (a) Comparison between DustCOMM (black line) and model simulations (coloured lines) of the globally averaged dust particle size distribution (PSD). The grey shading denotes the 95 % confidence interval for the DustCOMM product. (b) The probability distribution of the volume median diameter (µm) of the PSD for DustCOMM (black line) and the individual model simulations (coloured lines) over all locations and height levels.

While the globally averaged dust PSDs clearly show marked differences, it is also important to quantitatively examine the variability of the dust PSD for all locations. The variability of dust PSDs in the atmosphere is influenced by dust emission, transport, and deposition processes, and it can be assessed by considering metrics such as the volume median diameter (e.g. Maring et al., 2003; Formenti et al., 2011; Mahowald et al., 2014). Thus, the probability distributions of the volume median diameters (VMDs) for the model simulations are generally biased towards smaller VMD values, with different peak diameters for each model. WRFChem and IMPACT show the lowest VMD at ∼1.9 µm, and ARPEGE-Climat shows the highest VMD at ∼5.5 µm (Fig. 5a). In contrast, the DustCOMM VMD peaks around 5 µm. The probability distribution also shows that the DustCOMM VMD lies between approximately 2.5 and 6.5 µm at most heights and locations (Fig. 5b). This range is consistent with the range of measured VMDs (3–6 µm) for coarse-mode dust particles generally reported in the literature and compiled by Reid et al. (2003; see their Table 1). It also falls within the range of values measured at near-source regions and farther downwind. For instance, the VMD calculated from dust particle size distributions measured at Cabo Verde, off the coast of northern Africa (Ryder et al., 2018), is about 5.5 µm. Farther downstream where dust particles are likely to deposit after long-range transport, the VMD value near Puerto Rico is approximately 4 µm (Maring et al., 2003). It is noteworthy however that some studies (e.g. Carlson and Caverly, 1977; Weinzierl et al., 2009) have reported measured VMD values that exceed 13 µm, but these studies often include giant-mode dust particles with D≥20 µm, whereas we limited our analysis to dust with D≤20 µm (see Sect. 2.3.1). Overall, DustCOMM shows better consistency with observations of VMD than model simulations.

Changes in spatial variability of dust mass fraction

Although coarse-mode particles dominate the dust mass fraction near source regions and fine-mode particles dominate the dust mass fraction in the far remote regions, there are considerable changes in the spatial variability of the dust mass fraction between DustCOMM and the model ensemble (left and middle panels of Fig. 6). As highlighted above in Sect. 3.1.2, there is a general decrease in DustCOMM dust mass fraction for particles between 0.2–2.5 and 2.5–5 µm, relative to the model ensemble (right panel of Fig. 6). In contrast, there is an overall increase in DustCOMM dust mass fraction for particles between 5–10 and 10–20 µm. These changes cause DustCOMM to produce generally better agreement against in situ measurements than the model ensemble, as shown in Sect. 3.1.1 above. Overall, the most significant changes in DustCOMM dust mass fraction, relative to the model ensembles, are near dust-dominated regions, resulting in a decrease of up to 26 % and an increase of up to 29 % for dust particles between 2.5–5 and 10–20 µm, respectively.

Figure 6 Differences in the spatial variability of the dust mass fraction between DustCOMM and the model ensemble. Shown are the spatial distributions of the vertically integrated dust mass fractions for different particle bins for DustCOMM (left panel), the model ensemble (middle panel), and the difference between the two (right panel; DustCOMM−model ensemble).

These changes in the dust mass fraction gradually decrease away from the dust-dominated regions. This is evident, for example, over the North Atlantic basin, where dust from the Sahara is transported to the Caribbean and South America. Models generally simulate fewer large dust particles (D≥5 µm) and thus transport only a small fraction to the Caribbean. But observational evidence shown earlier in Fig. 2h and l indicates that dust in Barbados includes a significant fraction of coarse dust. Thus, the east–west gradient and the overall increase in the DustCOMM dust mass fraction over the North Atlantic help resolve the underestimation of long-range transported coarse particles, such as near Barbados (Fig. 6; e.g. Weinzierl et al., 2017).

Figure 7 Vertical distributions of dust mass fractions as a function of particle size for the individual model simulation (coloured lines) and DustCOMM (black lines) estimates. The grey shading shows the 95 % uncertainty confidence interval for DustCOMM.

The vertical distribution of the DustCOMM dust mass fraction shows differences with the model ensemble that are consistent with the globally averaged differences (Fig. 7) – that is, DustCOMM dust mass fractions are lower than for the model ensemble for particles between 0.2–2.5 and 2.5–5 µm and higher for particles between 5–10 and 10–20 µm. It is noteworthy here that vertical changes in the dust PSD in DustCOMM are based on model simulations, causing a similarity in the shape of the vertical profile of the dust mass fraction between DustCOMM and the model ensemble. Finally, similar changes in the spatial variability of the annually averaged dust mass fraction are apparent in the seasonally averaged values.

3.2 Dust mass extinction efficiency

3.2.1 Evaluation of DustCOMM against measurements

We evaluate the dust MEE (m² g⁻¹) of DustCOMM and the model ensemble against measurement (Fig. 8). These measurements span from those taken near dust source regions such as the Saharan, Middle East, and Asian deserts to those taken farther downwind from source regions (Table 2). Higher values of dust MEE are expected where fine-mode dust particles dominate, because smaller dust particles scatter light more efficiently per unit mass at visible wavelengths. In contrast, dust MEE decreases as the coarse-mode fraction increases. Thus, observed dust MEE values generally range between ∼0.3 and 0.8 m² g⁻¹ at approximately 550 nm.

Figure 8 Comparison of measurements (blue dots) of dust MEEs (m² g⁻¹) against column-integrated DustCOMM (black bars) and model ensemble (red bars) estimates. Vertical bars on the measurements represent reported uncertainty. For the DustCOMM and model ensemble estimates, the black and red boxes show 1 standard error, whereas the vertical dotted lines show the 95 % confidence interval; the middle horizontal bar and star show the median and mean values, respectively. The DustCOMM and model ensemble values are season-averaged values corresponding to the observation time period (see Table 2 for details). These seasons are labelled DJF – Dec–Feb, MAM – Mar–May, JJA – Jun–Jul, and SON – Sep–Nov; ANN represents an annually averaged value. The model ensemble MEE is calculated from the ratio between individual model dust aerosol optical depth and the dust mass loading, while the DustCOMM MEE is calculated using the constrained dust size distributions and single-particle extinction efficiency that takes into account the asphericity of dust aerosols. ${\mathit{\chi }}_{\mathit{ϵ}}^{\mathrm{2}}$ is the reduced chi-squared value (Eq. 10b) and quantifies the performance of a model in representing observations (e.g. Andrae et al., 2010).

DustCOMM shows better agreement with measurements of dust MEE than the model ensemble (Fig. 8). DustCOMM dust MEE estimates are within the measurement uncertainty range for most of the 11 studies used here. Notable exceptions are the comparison at Sde Boker, Israel (Andreae et al., 2002), and Qinghai Province, China (Li et al., 2000), where both DustCOMM and the model ensemble underestimate the measured MEE. Nevertheless, the DustCOMM estimates better reproduce the lower values of dust MEE near dust sources, and the higher values farther downstream. For example, lower dust MEE values near the Sahara, between Niamey and the Canary Islands (generally below 0.6 m² g⁻¹), and higher values farther downstream, such as over Barbados, are better reproduced by DustCOMM. DustCOMM dust MEE also compares well against measurements at the same location but for different seasons. An example is the measurements over Cabo Verde, off the coast of northern Africa (Haywood et al., 2003; Ryder et al., 2018), taken in September 2000 and August 2015. For both cases, DustCOMM estimates compare better with the observed dust MEE, while the model ensemble overestimates the values in both cases.

DustCOMM also reproduces the observed dust MEE values with a strong spatial gradient, measured during the same campaign (INDOEX) over the Arabian Sea and Indian Ocean (Quinn et al., 2002). Dust particles emitted from Middle East deserts can get transported over the Arabian Sea and are deposited over the Indian Ocean, where strong precipitation occurs year-round (e.g. Kulshrestha et al., 1996). Since dust MEE increases with distance from source regions due to deposition of larger dust particles, the measured dust MEE values increase from 0.5 m² g⁻¹ measured in the Arabian Sea to 0.75 m² g⁻¹ in the Indian Ocean, south of the Equator. DustCOMM captures much of this gradient and is in better quantitative agreement than the model ensemble estimate (Fig. 8).

DustCOMM also shows better agreement than the model ensemble against the observed dust MEE averaged over all measurements (see the last column of Fig. 8). DustCOMM shows a very small difference with the mean of the measurement estimates (0.007 m² g⁻¹; 95 % CI is −0.04 to −0.08), whereas the model ensemble mean (95 % CI) overestimates the measurements by 0.12 (−0.17–0.4) m² g⁻¹ – that is, about a 94 % reduction in the mean bias. We further assess DustCOMM performance by calculating the reduced chi square (${\mathit{\chi }}_{\mathit{ϵ}}^{\mathrm{2}}$; Eq. 11); a value of ${\mathit{\chi }}_{\mathit{ϵ}}^{\mathrm{2}}>\mathrm{1}$ highlights the degree that a model does not fit the observations within the uncertainty range (e.g. Andrae et al., 2010). DustCOMM shows a ${\mathit{\chi }}_{\mathit{ϵ}}^{\mathrm{2}}$ value of 1.19, in comparison to the model ensemble with ${\mathit{\chi }}_{\mathit{ϵ}}^{\mathrm{2}}$ value of 8.70 (Fig. 8), thereby showing a substantial improvement.

3.2.2 Global comparison between DustCOMM and model ensemble

After showing that DustCOMM better reproduces measurements of dust MEE than the model ensemble, we now compare the spatial variability of the DustCOMM and model ensemble dust MEE. To do so, we estimate the column-integrated dust MEE for DustCOMM and model ensembles over each location (Fig. 9a and b). Both DustCOMM and model estimates show smaller values of dust MEE over dust-dominated regions and higher values farther downwind – like over the Inter-Tropical Convergence Zone (ITCZ), the eastern Pacific Ocean, and the polar regions. Although DustCOMM and model ensemble estimates are thus spatially similar, important differences exist. Near dust-dominated regions, DustCOMM dust MEE values are lower than model ensembles, but farther downstream, DustCOMM values are higher than model ensembles. This regional difference in dust MEE values corresponds to a similar difference in dust mass fraction, with a fractional increase in coarse-mode dust over dust-dominated regions than farther downstream (compare Figs. 9 and 6). In addition, there is also a gradual east-to-west change in the dust MEE values as coarser dust particles are deposited away from dust sources, consistent with similar changes in dust mass fraction shown earlier in Fig. 6. The globally averaged DustCOMM dust MEE values are lower than predicted by the model ensemble. The global mean values of dust MEE for DustCOMM and model ensembles are 0.68 (min–max: 0.22–1.1) m² g⁻¹ and 0.95 (min–max: 0.30–1.98) m² g⁻¹, respectively.

Figure 9 Spatial distributions of column-integrated dust MEE (m² g⁻¹), weighted by the dust vertical distribution, for (a) DustCOMM, (b) the model ensemble, and (c) the difference between the two (DustCOMM−model ensemble).

3.3 Global comparison of atmospheric dust load between DustCOMM and models

After obtaining the DustCOMM dust MEE as described in the previous section, we combine this with the reanalysis-derived dust AOD (Eq. 9; see also Sect. 2.3.3) to obtain the atmospheric dust loading. We find that the DustCOMM dust column loading is generally larger than the model ensemble estimate (Fig. 10a and b). DustCOMM shows substantially larger dust column loading than the model ensemble over desert regions, such as the Middle East and Asian deserts. The relative increase in dust load in DustCOMM over the Asian desert is more than twice the increases over the Middle East desert. DustCOMM also shows larger dust column loading over most parts of the northern African desert, except some parts that include the north-western section and the coastal regions, which show smaller dust column loading than the model ensemble. Although reanalysis-derived mean dust AOD over northern Africa is substantially lower than the model ensemble, it is within the uncertainty estimates, which is higher over this region (see Supplement Fig. S4; see also Ridley et al., 2016). In addition, DustCOMM estimates over the Australian deserts show a lower dust column loading than the model ensemble, similarly corresponding to lower reanalysis-derived dust AOD (Figs. 10c and S4). Overall, globally averaged DustCOMM dust column loading is about 46 % higher than the model ensemble.

Figure 10 Spatial distributions of the atmospheric dust column loading (g m⁻²) for the (a) DustCOMM and (b) model ensemble estimates and (c) the difference between the two (DustCOMM−model ensemble).

3.4 Spatial distribution of DustCOMM relative uncertainty

We examine here the spatial distribution of the DustCOMM relative uncertainty – that is, the uncertainty characterizing 68 % of the distribution of each variable over each location divided by the mean value of that variable at that location. We do this for the dust mass fraction for the particle bins shown in Figs. 6 and 7, the dust MEE, and the dust load (Fig. 11).

Figure 11 Spatial distributions of DustCOMM relative uncertainties for (a–d) the dust mass fraction in the diameter ranges of 2–25, 25–5, 5–10, and 10–20 µm; (e) the dust MEE; and (f) dust load. The relative uncertainties are calculated as the ratio of the uncertainty characterizing 68 % of the distribution of each variable, divided by the mean value.

The relative uncertainties in the DustCOMM fine-mode fraction (D=0.2–2.5 µm) are higher mostly near emission regions (Fig. 11a), while the relative uncertainties in the coarse-mode fractions are higher over remote regions, especially for D=10–20 µm (Fig. 11d). These uncertainties are, in part, directly associated with the uncertainties in the measurement constraints. The globally averaged constrained dust size distribution (Eq. 1) has a higher relative uncertainty for the D≤1 µm and D≥10 µm diameter range than for the $\mathrm{1}\le D\le \mathrm{10}\mathrm{\mu }\mathrm{m}$ diameter range (see Fig. 2 in Kok et al., 2017), and we propagate these uncertainties over every location. In addition, the spatial distribution for the relative uncertainties in the dust mass fraction is similar to that of the model ensembles (Supplement Fig. S5), which is also propagated into the DustCOMM product.

The relative uncertainties in DustCOMM dust MEE are mostly higher over dust-dominated regions (Fig. 11e). The dust MEE is influenced by the uncertainty in the constrained globally averaged extinction efficiency, which in turn is partially due to uncertainties in the in situ emission measurements of the index of refraction and dust particle shapes (see Fig. 1b in Kok et al., 2017), all of which are propagated into the DustCOMM dust MEE. In addition, the relative uncertainties in the dust MEE are also affected by the uncertainty in the dust size distribution. Thus, the spatial distribution of dust MEE relative uncertainty is particularly informed by the uncertainties in the fine-mode and coarse-mode dust particles (compare Fig. 11a and d with e). For the most part, uncertainties in the fine-mode dust fraction appear to dominate the uncertainties in dust MEE, more than the uncertainties in the coarse-mode dust fractions.

The relative uncertainties in the DustCOMM dust column loading are mostly higher over remote regions, where the mean dust load is small (Fig. 11f). Though the dust column loading is influenced by the uncertainties in dust MEE, the spatial distribution of the relative uncertainties in dust load is largely informed by the uncertainties in the reanalysis dust AOD (see Supplement Fig. S4).

4 Discussion

We presented the DustCOMM products in the previous section, where we showed that both the dust PSD and the dust MEE are reproduced more accurately than by an ensemble of model simulations. Despite the overall agreement with observations, there are some disagreements highlighting potential limitations of our methodology. In this section, we discuss these disagreements between DustCOMM and measurements and provide possible insights into these discrepancies (Sect. 4.1). We also discuss the impact of dust sizes and asphericity on DustCOMM dust mass extinction efficiency (Sect. 4.2), and we highlight the limitations in using modelling constraints as part of DustCOMM estimates (Sect. 4.3). We end by highlighting how our constrained DustCOMM products can be used by the research community to potentially improve estimates of dust impacts on the Earth system (Sect. 4.4).

4.1 Cause of discrepancy between DustCOMM and size distribution measurements

The evaluation of the DustCOMM PSD shows an underestimation of dust with D≤0.5 µm and D≥10 µm (Figs. 2, 3, and 4). This is in contrast to the ensemble of model simulations overestimating the dust mass fractions for D≤0.5 µm and underestimating the dust mass fraction substantially more than DustCOMM for D≥10 µm. Although the comparison between date-specific individual measurements and season-averaged DustCOMM dust PSD is expected to induce errors, this difference cannot explain the apparently systematic difference between measurements and the DustCOMM dust PSD for both D≤0.5 µm and D≥10 µm (Fig. 4). We provide here possible reasons for this disagreement between DustCOMM and observations.

First, DustCOMM's underestimation of dust with D≤0.5 µm may be caused by contamination of the measured size distributions by other aerosol species for D≤0.5 µm. Studies have shown that a substantial fraction of aerosols with D≤0.5 µm are not mineral dust, even in dust-dominated regions (Chou et al., 2008; Kandler et al., 2009; Weinzierl et al., 2009). For example, during the Saharan Mineral Dust Experiment (SAMUM) over southern Morocco, Kandler et al. (2009) showed that more than 50 % of the measured particles with D≤0.5 µm are ammonium sulfates or a mixture of sulfate and dust. Even when strict measurement techniques are used to separate other non-mixing aerosol components, the aerosol mixing state for D≤0.5 µm often leads to outer coating of available dust particles, thus leading to a higher particle volume that overestimates the true dust size (Weinzierl et al., 2009). In addition, campaign logistics often require that some measurements of dust properties are taken close to major cities, where contaminations by other aerosol species, such as biomass-burning aerosols or urban pollutions, are possible (e.g. McConnell et al., 2008; Wagner et al., 2009). For example, Clarke (2004) highlighted that the presence of biomass-burning aerosols (e.g. soot) led to a variability of about 2 orders of magnitude for measured size distributions with diameters less than D≤0.6 µm during the ACE-Asia campaign. This variability is consistent with the average difference between our estimates and the observations for D≤0.5 µm. After separating out the contamination of the soot mode from the dust size distribution, their resulting dust PSD generally agrees with our estimate within the uncertainty range (Fig. 2f). Thus, the large variability of the measured size distribution is indicative of the potential problems with the representation of dust particles with D≤0.5 µm.

Second, the constraint on the globally averaged dust size distribution could also underestimate the contribution from dust with D≤0.5 µm. A key input to this constraint is the emitted dust size distribution, but there is a dearth of measurements of the mass fraction of emitted dust with D≤0.5 µm, leading to uncertainty in constraining the globally averaged emitted dust size distribution with D≤0.5 µm (Kok et al., 2017). Moreover, the measurements of emitted dust size distribution with D≤0.5 µm that do exist (e.g. Fratini et al., 2007; Sow et al., 2009; see Fig. 1c in Kok et al., 2017) indeed show a larger dust mass fraction than represented in the constraint on the globally averaged emitted dust size distribution. Therefore, more measurements of the size distribution of emitted dust particles extending to very fine sizes are needed.

Third, the underestimation of dust with D≥10 µm by both DustCOMM and the model ensemble might be caused by biases in both global model simulations and the constraints on the global dust size distribution used by DustCOMM. Similar to D≤0.5 µm, the experimental constraint on the emitted dust size distribution with D≥10 µm also has a large uncertainty because of limited available measurements (Kok, 2011a). In addition, since spatial and temporal variabilities of large dust particles (D≥10 µm) strongly depend on the model simulation of dust emission and deposition processes, uncertainties in these processes will influence the constraints on DustCOMM dust size distribution. For example, if the giant mineral dust particles are transported far away from the source regions as suggested by observations (e.g. van der Does et al., 2018), the lack of this mechanism would result in a negative bias of the simulated dust atmospheric lifetime (e.g. Huneeus et al., 2011). And since modelling constraints of globally averaged dust lifetime are used to constrain the globally averaged size distribution (Eq. 1), such systematic negative bias may have contributed to the underestimation of dust particles with D≥10 µm. Although our methodology partly constrains dust deposition globally, it does not constrain regional variability in dust deposition, and we expect that such uncertainties may increase as a function of distance away from dust-dominated regions. We note here that regional observational constraints on dust lifetime are currently not available, and stronger modelling constraints that may account for the underestimation of coarse dust particles in the atmosphere are a subject for future work.

4.2 Impacts of dust sizes and asphericity on DustCOMM dust mass extinction efficiency

The dust MEE is partially determined by the dust size distribution (Eq. 8). Despite the good agreement between DustCOMM and the measurements of dust MEE (Fig. 8), the size discrepancies in the dust size distribution for particles with D≤0.5 µm and D≥10 µm (Figs. 2, 3, and 4) affect the estimation of dust MEE. Dust with D≤0.5 µm has a large single-particle MEE, whereas dust with D≥10 µm has a small single particle MEE (see Supplement Fig. S6). Consequently, errors due to the possible overestimation of both size fractions at least partially cancel each other.

In addition to the impact of dust sizes, dust asphericity also has a substantial impact on the dust MEE. The DustCOMM constraint on dust MEE leverages measurements of dust shape to represent dust particles as an ensemble of tri-axial ellipsoids (Meng et al., 2010; Kok et al., 2017). In contrast, most models use Mie theory, which approximates dust as spherical particles. Thus, the difference between single-particle dust MEE used in DustCOMM and calculated using Mie theory shows the impact of dust asphericity is substantial for both small and lager dust particles, increasing extinction for particles with D≥1 µm (Supplement Fig. S6). This implies that in typical global model simulations, which contain too many fine-mode dust particles and approximate dust as spherical, the overestimation of the dust extinction due to the fine size bias could (partially) cancel out the underestimation of the dust extinction due to the treatment as dust spherical shapes, leading to nonetheless reasonable agreement with measured dust MEE. However, for DustCOMM, both the size bias and dust asphericity are accounted for, thereby producing better agreement with measurements (Fig. 8). In addition, accounting for dust asphericity could allow dust particles to stay longer in the atmosphere because asphericity reduces dust settling speed (Ginoux, 2003), which may in turn lead to a more accurate estimation of dust deposition mass fluxes onto land and ocean ecosystems (e.g. van der Does et al., 2016, 2018).

4.3 Limitations in using modelling constraints

We used modelling constraints in DustCOMM where observational constraints were either not available or insufficient. For example, modelling constraints are used for the regional differences in dust size distribution and extinction efficiency because the measurements to constrain these parameters on a regional basis across the different dust source regions are currently insufficient. To further reduce the uncertainty associated with using modelling constraints, we used an ensemble of six model simulations. In addition to the uncertainties associated with model simulations of dust emission and deposition processes that may influence the constraints on dust size distributions as highlighted in Sect. 4.1, there are other limitations in the modelling constraints that can influence DustCOMM estimates.

First, one such limitation is the uncertainty in the dust mass spatial distribution of the model ensemble, which directly determines the spatial distribution of dust mass for DustCOMM estimates. Variability in dust emission rates influences the distribution of simulated size-resolved atmospheric dust loading and consequently the 3-D dust mass fractions. In addition, ensemble model simulations of dust emission and transport are driven by different meteorological datasets (Table 1), which represent the actual historical meteorology with various degrees of accuracy (e.g. Evan, 2018). Dust transport is also influenced by model resolution and sub-grid parameterizations of wind and turbulence, which differ between models (e.g. Zender et al., 2003; Cakmur et al., 2004). Although averaging over multiple models and over long time periods reduces random errors, systematic errors that affect different models similarly would affect the model ensemble (e.g. Ridley et al., 2012) and would impact the spatial distribution of dust mass (e.g. Johnson et al., 2012; Ridley et al., 2016). In addition, uncertainties in the vertical distribution of size-resolved dust mass fractions directly affect DustCOMM dust size distributions. Since we use the globally averaged size-resolved extinction efficiency to constrain the dust MEE over every location (Eq. 8), the spatial distribution of dust MEE is thus partially determined by the dust size distribution, effectively propagating any uncertainty in model simulations of the spatial distribution.

Second, some errors may have been introduced while scaling and fitting the different model dust size distributions to a common diameter range (Sect. 2.3.1). For the scaling procedure (Sect. 2.3.1), the variances of the dust mass fraction in all the bins, including the newly created ones, are of similar orders of magnitude, and thus errors introduced through this process are small relative to the magnitude of errors in the dust mass fraction. In addition, the resulting dust size distributions are dependent on the specific function and set of parameters used in the fitting procedure (Sect. 2.3.1), which may also introduce some errors.

Third, our constraints on the dust atmospheric loading use ensemble estimates of reanalysis-derived dust AOD, which depends in part on the assimilated aerosol observations, in part on the numerical simulation of dust sources and sinks, and in part on the numerical simulation of other aerosol species. Although some of the reanalysis products try to constrain these dust processes using space-based observations (e.g. Lynch et al., 2016; see Supplement), the impact of the uncertainties associated with each process on the DustCOMM estimates of the atmospheric dust loading is beyond the scope of the study.

Finally, this study primarily uses climatologies of modelled dust size distribution between 2004 and 2008, except for WRF-Chem and IMPACT (see Table 1), and it also scales dust mass loading using the 2004–2008 reanalysis products (see Sects. 2.3.3 and 2.2). Thus, any application of our methodology to a different time period is expected to have some errors. While these errors are expected to be small for the dust size distribution and dust mass extinction efficiency, they may have a substantial impact on the dust mass loading, depending on the inter-annual variability in the reanalysis-derived products and also on the assimilated observations.

4.4 Possible use of DustCOMM to improve estimates of dust impacts on the Earth system

Given that DustCOMM estimates of dust aerosol properties are in better agreement with measurements than the model ensemble, DustCOMM could be used to obtain improved constraints on dust impacts on the Earth system than is possible from current global models. Specifically, DustCOMM dust properties could be used as an alternative to global model simulations in constraining dust impacts, such as the dust direct radiative impact or dust impacts on biogeochemistry and human health. For instance, dust radiative heating rates in the atmosphere strongly depend on the ability of dust particles to absorb shortwave and longwave radiation (e.g. Perlwitz and Miller, 2010). In turn, such absorption depends on the dust size distribution, which strongly influences the optical parameters like the dust absorption optical thickness (e.g. Tegen and Lacis, 1996). With improved constraints on the dust size distribution and therefore the dust optical properties, DustCOMM could be used to determine the dust (shortwave and longwave) heating rates in the atmosphere more accurately than is possible with current global model simulations. As a result, our constraints on dust size distribution could be used to better quantify radiative effects of dust, especially in the longwave spectrum, which have remained very uncertain (Dufresne et al., 2002; Di Biagio et al., 2017a; Kok et al., 2017; Song et al., 2018). Furthermore, since recent studies associate much of the biases in dust properties, such as the dust aerosol optical depth, deposition fluxes, and surface dust concentration, with model biases in dust size distribution (Evan et al., 2014; Huneeus et al., 2011), DustCOMM estimates can therefore serve as a better alternative. For example, DustCOMM's improved constraints on atmospheric dust loading and dust size distribution could contribute to better estimates of size-resolved dust concentration near the surface (e.g. Whicker et al., 2018). Over the ocean, such constraints on size-resolved dust concentration could potentially be used for constraints on dust deposition fluxes that are more accurate than is possible from global model simulations.

In addition to being used as an alternative to global model simulations, DustCOMM could also be used to improve the simulation of dust aerosol properties in global models. Incorporating DustCOMM products into the simulation process can potentially be achieved when the aerosol module is coupled with the global model in either the so-called online or offline modes (e.g. Tegen, 2003; Pérez et al., 2011; Han et al., 2012). In the online mode, the simulated dust size distributions could be adjusted (“nudged”) to match the DustCOMM constraints on dust size distribution, similar to what is often done with meteorological fields (e.g. Kooperman et al., 2012). Alternatively, the 3-D dust size distribution could also be corrected offline after the simulated size distribution is obtained but before dust impacts such as on radiation are estimated (e.g. Weaver et al., 2002). Specifically, the modelled dust size distribution can be corrected by minimizing the differences between the DustCOMM and modelled size distributions for a specific timescale (see Sect. 2.3.5). Whether simulated dust properties are corrected in the online or offline modes, using DustCOMM to bias-correct global model simulations could produce better estimation of dust impacts, such as dust impacts on radiation, clouds and precipitation, biogeochemistry, and human health.

An example of dust impacts that can be substantially improved by the DustCOMM product are dust radiative effects. These radiative effects are sensitive to dust particle sizes and shapes, which are both constrained substantially more accurately in DustCOMM than in models (Figs. 2–6). Smaller dust particles (D≤2.5 µm) scatter more shortwave radiation and cool the climate, while larger dust particles (D≥5 µm) absorb more longwave radiation and warm the climate. Thus, correcting both biases of too much fine dust and not enough coarse dust in models (Figs. 4 and 5), as we do here in DustCOMM, decreases the shortwave cooling and increases the longwave warming (e.g. Otto et al., 2011; Kok et al., 2017). Using the 3-D DustCOMM size distribution to correct modelled dust properties could yield more accurate estimates of dust radiative effects.

In addition, simulated dust impacts on clouds and precipitation can also be improved using DustCOMM dust aerosol properties. For the interactions of dust particles with clouds, it is important to know the number of particles that are activated above a given particle size as cloud condensation nuclei or ice nuclei (e.g. Andreae and Rosenfeld, 2008; DeMott et al., 2015). Therefore, in regions with significant dust loading, accurate estimates of dust size distribution can be key to accurate simulations of precipitation initiation and aerosol–cloud interactions, including dust aerosol indirect and semi-direct effects (e.g. Sassen, 2003; Doherty and Evan, 2014). Since DustCOMM represents the dust size distribution more accurately than model simulations, it could be used to improve the simulated dust impacts on clouds, precipitation, and aerosol–cloud interactions.

Another key advantage of DustCOMM over global model simulations is that it propagates many observational, experimental, and modelling uncertainties of dust properties, which can be propagated into the calculation of dust impacts on the Earth system. For instance, experimental uncertainties associated with the emitted dust size distributions are propagated into the DustCOMM 3-D dust size distribution, and experimental uncertainties in the dust index of refraction and dust particle shapes are propagated into the DustCOMM mass extinction efficiency at 550 nm wavelength (e.g. Kok et al., 2017). In addition, our methodology propagates the uncertainty due to the spread in model predictions of the dust spatial distribution, although substantial biases in the model ensemble might exist (see Sect. 4.3 for example).

Finally, it is worth noting that DustCOMM can be readily updated as more accurate constraints on dust properties and abundance become available. Current constraints in DustCOMM can also be expanded to include more information about dust properties. For instance, a next step could be to include constraints on the dust vertical concentration profile over every location, in order to more accurately estimate dust deposition, and dust concentration at the surface and in 3-D. For this, lidar-based retrieval of vertical dust extinction profiles from Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observations (CALIPSO) can be combined with the corresponding constraints on dust mass extinction efficiency from this study to obtain constraints on the dust vertical concentration profile. Another addition could be constraining the relative contribution of each dust source region to the 3-D dust load, which can be combined with constraints on optical properties of dust emitted from each region (Di Biagio et al., 2017b, 2019; Green et al., 2018) to obtain more accurate quantifications of dust radiative impacts. Given that dust particles with D≥20 µm can contribute substantially to dust extinction in both the shortwave and longwave spectra (Ryder et al., 2019), future versions of DustCOMM could be extended to a diameter range beyond 20 µm as more measurements of dust size distribution with D≥20 µm become available.

5 Summary and conclusions

In this study, we presented a new dataset of atmospheric dust aerosol properties called Dust Constraints from joint Observational-Modelling-experiMental Analysis – DustCOMM. DustCOMM combines observational and experimental constraints on dust properties and abundance with an ensemble of global model simulations of dust spatial distribution to obtain more accurate 3-D annual and seasonal climatologies of dust properties and abundance than is possible with global model simulations alone. Here, we presented three DustCOMM products: the three-dimensional (3-D) dust size distribution, 3-D dust mass extinction efficiency, and two-dimensional (2-D) dust loading. First, we obtained constraints on the 3-D dust size distribution by combining constraints on the globally averaged dust size distribution with an ensemble of model simulations of the spatial variability of the dust size distribution. Second, we combined the resulting 3-D dust size distribution with constraints on the size-resolved globally averaged dust extinction efficiency, which accounts for the substantial asphericity of dust aerosols, to constrain the 3-D dust mass extinction efficiency. Finally, we used the resulting column-integrated dust mass extinction efficiency with an ensemble of reanalysis-derived dust aerosol optical depth to constrain the atmospheric dust column loading.

By comparing DustCOMM estimates of dust size distribution and dust mass extinction efficiency against independent in situ measurements, we showed that DustCOMM reproduces observations substantially better than an ensemble of model simulations (Figs. 2–4 and 8). Models generally overestimate the contribution of fine-mode dust (D≤2.5 µm) and underestimate the contribution of coarse-mode dust (D≥5 µm), consistent with previous studies (e.g. Mahowald et al., 2014; Kok et al., 2017). In contrast, the DustCOMM size distribution is in substantially better agreement with measurements for different locations, heights, and seasons over the $\mathrm{0.5}\le D\le \mathrm{20}\mathrm{\mu }\mathrm{m}$ size range. However, there remain some discrepancies between DustCOMM and measurements, notably an underestimation of dust with D≤0.5 µm. Potential reasons for these discrepancies include contamination of measured dust size distribution by other aerosol species for D≤0.5 µm and biases in observational and modelling constraints for D≤0.5 µm (Sect. 4.1). Because DustCOMM underestimates the measurements for D≤0.5 µm, it shows a more negative bias (∼50 % more) over the full size range (between D=0.2 and 20 µm), although the error is markedly lower (∼15 %) when compared to the ensemble of model simulations. Overall for D≥0.5 µm, DustCOMM shows a bias against measured size distributions that is significantly less (about 46 % less) than for an ensemble of global model simulations.

DustCOMM similarly shows better agreement against measurements of the dust mass extinction efficiency (MEE) than an ensemble of model estimates. Because DustCOMM predicts a coarser dust size distribution, as supported by the comparison against in situ size distribution measurements, it yields a global-mean dust MEE that is about 28 % lower than that from the model ensemble, driven by large reductions in MEE over dust-dominated regions, where coarse particles dominate. For specific locations and seasons, DustCOMM estimates consistently show smaller errors relative to dust MEE measurements than an ensemble of model results, including in regions with strong spatial gradients in dust loading. On average, there is a negligible difference (∼1 %) between DustCOMM and measurements of MEE, while the model ensemble overestimates MEE by about 23 % relative to measurements.

DustCOMM estimates of spatially varying dust properties and abundance can be used to constrain various dust impacts on the Earth system in a manner that is more robust than possible with current global models. This is because DustCOMM reproduces dust properties more accurately than global model simulations and also because DustCOMM explicitly propagates uncertainties in experimental, observational, and modelling constraints used in obtaining the DustCOMM products, and these uncertainties can be propagated in calculations of dust impacts on global climate, biogeochemistry, and human health.

Appendix A: List of some acronyms

GISS	Goddard Institute for Space Studies (GISS) ModelE atmospheric general circulation model
WRF-Chem	Weather Research and Forecasting model coupled with Chemistry
CESM	Community Earth System Model
GEOS-Chem	Goddard Earth Observing System coupled with Chemistry
IMPACT	Integrated Massively Parallel Atmospheric Chemical Transport
INDOEX	Indian Ocean Experiment Intensive Field Phase
SHADE	Saharan Dust Experiment
ACE-Asia	Asian Pacific Regional Aerosol Characterization Experiment
TRACE-P	Transport and Chemical Evolution over the Pacific
ACE-2	Aerosol Characterisation Experiment
DABEX	Dust and Biomass-burning Experiment
AMMA	African Monsoon Multidisciplinary Analysis
DODO	Dust Outflow and Deposition to the Ocean project
SAMUM	Saharan Mineral Dust Experiment
DARPO	Desert Aerosols over Portugal
BACEX	Barbados Aerosol Cloud Experiment
SALTRACE	Saharan Aerosol Long-Range Transport And Aerosol–Cloud-Interaction Experiment
AER-D	AERosol Properties – Dust
MERRA-2	Modern-Era Retrospective analysis for Research and Applications, Version 2
NAAPS	Navy Aerosol Analysis and Prediction System
JRAero	Japanese Reanalysis for Aerosol
CAMSiRA	Copernicus Atmosphere Monitoring Service interim Reanalysis