Supplement of Improved representation of the global dust cycle using observational constraints on dust properties and abundance

Abstract. Even though desert dust is the most abundant aerosol by
mass in Earth's atmosphere, atmospheric models struggle to accurately
represent its spatial and temporal distribution. These model errors are
partially caused by fundamental difficulties in simulating dust emission in
coarse-resolution models and in accurately representing dust microphysical
properties. Here we mitigate these problems by developing a new methodology
that yields an improved representation of the global dust cycle. We present
an analytical framework that uses inverse modeling to integrate an ensemble
of global model simulations with observational constraints on the dust size
distribution, extinction efficiency, and regional dust aerosol optical
depth. We then compare the inverse model results against independent
measurements of dust surface concentration and deposition flux and find that
errors are reduced by approximately a factor of 2 relative to current
model simulations of the Northern Hemisphere dust cycle. The inverse model
results show smaller improvements in the less dusty Southern Hemisphere,
most likely because both the model simulations and the observational
constraints used in the inverse model are less accurate. On a global basis,
we find that the emission flux of dust with a geometric diameter up to 20 µm (PM20) is approximately 5000 Tg yr−1, which is greater than most
models account for. This larger PM20 dust flux is needed to match
observational constraints showing a large atmospheric loading of coarse
dust. We obtain gridded datasets of dust emission, vertically integrated
loading, dust aerosol optical depth, (surface) concentration, and wet and
dry deposition fluxes that are resolved by season and particle size. As our
results indicate that this dataset is more accurate than current model
simulations and the MERRA-2 dust reanalysis product, it can be used to
improve quantifications of dust impacts on the Earth system.


This document contains a description of the atmospheric model simulations used in this study, namely CESM/CAM4, IMPACT, ModelE2, GEOS/GOCART, MONARCH, and INCA/IPSL. We also provide a description of how we account for correlation of errors in the constraint on the regional dust aerosol optical depth used in our inverse model. Finally, this document also contains a number of Supplementary  30 Figures, which are summarized below: • Figure S1. Spatial distribution of annually averaged bulk dust loading (g/m 2 ) produced by a unit (1 Tg) of total atmospheric dust loading produced by each of the nine source regions, as simulated by each of the six models in our model ensemble. • Figure S2. Inverse model results for the seasonally resolved dust emission rate. 35 • Figure S3. Inverse model results for the seasonally resolved dust AOD.
• Figure S4. Inverse model results for the seasonally resolved dust column loading.
• Figure S5. Inverse model results for the seasonally resolved dust surface concentration.
• Figure S6. Inverse model results for the seasonally resolved dust deposition flux.
• Figure S7. Bulk global dust emission rate (a) and loading (b) as a function of the maximum 40 simulated dust diameter. • Figure S8. Comparison of model simulations against measurements of NH seasonal surface concentration. • Figure S9. Comparison of model simulations against measurements of NH annual surface concentration. 45 • Figure S10. Comparison of model simulations against NH deposition flux measurements.
• Figure S11. Comparison of model simulations against measurements of SH seasonal surface concentration. • Figure S12. Comparison of model simulations against measurements of SH annual surface concentration. 50 • Figure S13. Comparison of model simulations against SH deposition flux measurements.

Description of CAM4 simulations
We used the Community Atmosphere Model version 4 (CAM4) within the Community Earth System Model version 1.0.5 (CESM 1.0.5) to simulate the emission, transport and deposition of dust mineral 55 aerosols (Hurrell et al., 2013) in order to calculate their radiative effects. Within CAM4, dust aerosols are treated as externally mixed and separated by particle size into four size bins (0.1-1.0, 1.0-2.5, 2.5-5.0, 5.0-10.0 µm) using a bulk aerosol model (BAM) scheme. We extended the CAM4 simulation results by adding a 10-20 µm bin following the procedure in Adebiyi et al. (2020; where ̃C AM4 and ̃G C denote normalized (i.e., per 1 Tg dust loading) dust fields (loading, concentration, and deposition) in respectively CAM4 and GOCART. As such, Eq. S1 uses the difference between the spatial distributions of the 12-20 and 6-12 µm particle bin dust fields in GOCART to approximate the difference between the spatial distributions of the 10-20 and 5-10 µm particle bin dust fields in CAM4. 65 Dust aerosol was subdivided into eight mineral species: calcite, feldspar, gypsum, hematite, illite, kaolinite, quartz, and smectite (Scanza et al., 2015;Scanza et al., 2018), with each mineral being advected as its own tracer in the model. The dust emission mineralogy was prescribed from offline soil mineralogy maps (Claquin et al., 1999). For this study, we also configured the Dust Entrainment and Deposition Model (DEAD) (Zender et al., 2003) to include the dust mobilization parameterization of Kok et al. 70 (2014), which removes the offline soil erodibility map requirement and improves the model dust concentration and deposition comparisons with observations (Hamilton et al., 2019). Dust aerosol optical depth (AOD) was tuned to attain a global annual average of 0.03, as suggested by Ridley et al. (2016). All simulations used a model resolution of 1.9° x 2.5° (latitude by longitude) with 56 vertical levels from the surface to 2hPa and were run from 2003-2008. The last five years (2004)(2005)(2006)(2007)(2008) were used for analysis. 75 The model dynamics were forced offline by the European Centre for Medium-Range Weather Forecast's (ECMWF) global atmospheric reanalysis ERA-Interim meteorology (Dee et al., 2011).

Description of IMPACT simulations
This study used the Integrated Massively Parallel Atmospheric Chemical Transport (IMPACT) model to calculate the concentration of mineral dust aerosols in 4 size bins (diameters: 0.1-1.26, 1.26-2.5, 2.5-5, 80 and 5-20 µm) (Ito et al., 2020 and references therein using a physically-based emission scheme in conjunction with satellite products of vegetation cover and soil moisture in the model (Kok et al., 2014;Ito and Kok, 2017). The satellite-based estimates of fractional vegetation area in conjunction with land cover type were used to account for suppressing effects of vegetation on dust emission in barren and open shrublands (Webb et al., 2014;Ito and Kok, 2017). The chemical composition of mineral dust aerosols can change dynamically from that in the 95 originally emitted aerosols due to reactions with gaseous species. The aging of mineral dust aerosols from hydrophobic to hydrophilic can enhance their gravitational settling, dry and wet deposition (Wang and Penner, 2009;Ito et al., 2020). Atmospheric concentrations and depositions of aerosols have been evaluated extensively on global and regional scales (Myriokefalitakis et al., 2018;Ito et al., 2019). To isolate the contribution of each of the specified source regions to the global dust loading, we tagged the 100 dust originating from each source region.

Description of ModelE2.1 simulations
The NASA Goddard Institute for Space Studies (GISS) ModelE2.1 calculates the emission, transport and removal of dust aerosols in five size bins, whose diameters span 0.1-2 µm, 2-4 µm, 4-8 µm, 8-16 µm and 16-32 µm, respectively. The smallest bin corresponds to clay sized particles, while the remainder 105 represent silt sizes. We did not use the largest bin (16-32 µm) because it exceeds the 20 µm maximum diameter used in the inverse model and instead generated a 16-20 µm bin based on the 8-16 µm bin and the GOCART 12-20 µm bin, as follows (also see description of the CAM4 simulations and Adebiyi et al., 2020): ModelE2.1 has horizontal resolution of 2° latitude by 2.5° longitude and 40 vertical layers that extend to 0.1 hPa just above the stratopause (Bauer et al., 2020;Kelley et al., 2020).
Dust emission occurs where loosely bound soil particles are abundant and vulnerable to wind erosion in the absence of vegetation. Regions of dust emission are prescribed using the source map of Ginoux et al. (2001) that combines a satellite vegetation mask with a geomorphological criterion identifying sources as 115 basins where soils are replenished by erosion from surrounding highlands. This map results in a distribution of dust sources that corresponds to well-known dust "hot spots": localized regions where aerosol optical depth is frequently high (Prospero et al., 2002). Soil particles are emitted when the wind speed exceeds a prescribed threshold that increases with soil moisture. The model also represents emission by wind gusts that are ephemeral compared to the duration over which surface wind speed is 120 updated by the model (Cakmur et al., 2004). The cubic dependence of emission upon wind speed is based upon wind tunnel measurements (Shao et al., 1993). The mass of emitted dust for a given wind speed is calibrated using a worldwide collection of measurements . In particular, the emitted ratio of clay and silt particles is prescribed to match retrievals of the aerosol size distribution at several AERONET stations in regions of high dust concentration. This results in a ratio consistent with 125 measurements of the emitted size distribution compiled by Kok (2011).
Dust aerosols are transported by the model winds, which in this study are relaxed to four-times daily NCEP values (Kalnay et al., 1996). The relaxation time at all model levels is 1000 s (just under 20 minutes), a duration chosen to match the observed circulation including vertical motion. The NCEP winds are linearly interpolated in time between their successive 6-hour values. 130 Dust is removed from the atmosphere through turbulent removal near the surface and gravitational settling (Wesely and Hicks, 1977;Koch et al., 1999). In addition, removal occurs within clouds when particles nucleate cloud droplets and ice crystals. Dust is also removed both within and below clouds when particles are captured through collisions with precipitation. Falling droplets that completely reevaporate during descent release the dust particle (Bauer and Koch, 2005;Schmidt et al., 2006). 135 More information about the calculation of the dust life cycle is given by Miller et al. (2006) and Perlwitz et al. (2015). The dust distribution has been compared to an extensive collection of measurements, for example by Cakmur et al. (2006), Miller et al. (2006), and Huneeus et al. (2011).

Description of GEOS/GOCART simulations
This study used the NASA Goddard Earth Observing System (GEOS) global Earth systems model, which 140 incorporates coupled components for atmospheric and oceanic circulation and atmospheric chemistry and composition, and state-of-the-art data assimilation capabilities for both atmospheric dynamics and aerosols (Rienecker et al., 2008). Aerosol processes are represented using a version of the Goddard Chemistry Aerosol Radiation and Transport (GOCART) (Chin et al., 2002;Colarco et al., 2010). GOCART is run online within GEOS in a configuration that replays meteorology from a previous 145 atmospheric analysis, in this case from the Modern-Era Retrospective analysis for Research and Applications, version 2 (MERRA-2) (Gelaro et al., 2017;Randles et al., 2017). In this study, simulations were run on a cubed-sphere horizontal grid at ~100 km resolution with 72 vertical hybrid-sigma levels extending from the surface to 1Pa (around 80 km). GOCART models the emissions, mixing, chemistry, and deposition of major aerosol components in the troposphere, including sulfates, nitrates, dust, black 150 carbon, organic carbon, and sea-salt aerosols. Dust emissions were calculated using an updated version of the scheme described in Ginoux et al. (2001). Dust is described with 8 size bins spaced between 0.1 -10 µm radius. The four submicron particles size bins (radius < 1µm) are grouped in one unique bin for emission and transport purposes. The initial dust particle size distribution was tuned to follow the brittle fragmentation theory approach of Kok (2011). Dust optical properties are modeled using spheroidal shape 155 assumptions, after Colarco et al. (2014).

Description of MONARCH simulations
The Multiscale Online Non-hydrostatic AtmospheRe CHemistry model (MONARCH; previously known as NMMB/BSC-CTM) is a model developed at the Barcelona Supercomputing Center (e.g., Pérez et al., 2011;Badía et al., 2017). MONARCH contains advanced chemistry and aerosol packages, and it is coupled 160 online with the Non-hydrostatic Multiscale Model (NMMB), which allows for running either global or high-resolution (convection-permitting) regional simulations (Janjic et al., 2001;Janjic and Gall, 2012). with friction velocity u* and threshold friction velocity u*t. The fraction of particles in size bin d is given by 180 fk and follows Kok (2011). S is the topographic source function described by Ginoux et al. (2001). The threshold friction velocity for a smooth and dry surface, u*t0, is calculated as the minimum value obtained using the formulation from Shao and Lu (2000). Threshold friction velocity, u*t, is then obtained by applying corrections for the effects of soil moisture (fw) and vegetation cover (fv) as * = * 0 • • 185 The correction factor for soil moisture, fw, follows Ginoux et al. (2001). fv is based on Raupach et al. (1993) and requires frontal area index λ as an input, which we estimate based on monthly satellite-based retrievals of photosynthetic and non-photosynthetic vegetation cover (Guerschman et al., 2015) using the logarithmic relationship (Shao et al., 1996) = −c λ ln(1 − η) 190 with fractional vegetation cover η and cλ = 0.35. The photosynthetic vegetation fraction is also used in the land-surface modules, which warrants consistency within MONARCH. Potential dust source areas are specified using the frequency of occurrence of MODIS Deep Blue dust optical depth larger than 0.2 (Ginoux et al., 2012).
The model includes eight dust size transport bins ranging up to 20 µm in diameter. The effective and 195 volume radii of the bins used in the radiative and sedimentation schemes respectively are time-invariant and based on a lognormal distribution with mass median diameter of 2.524 μm and geometric standard deviation of 2.
The model runs performed with MONARCH are conducted at a horizontal resolution of 1.0° latitude by 1.4° longitude with 48 vertical layers and a computational time step of 3 min. Turbulence, surface layer, 200 dust emission, sedimentation and dry deposition routines are called every 4 computational times steps (12 min), moist convection, microphysics and wet scavenging routines every 8 time steps (24 min), and shortand longwave radiation routines are called every 20 time steps (60 min). MONARCH runs are initialized using ERA Interim reanalysis data (Dee et al., 2011). The meteorological fields are re-initialized daily, while dust fields and soil moisture are transferred between the daily runs. We used one year of spinup for 205 soil moisture and one month of spinup for the dust fields before the 5-year simulation time period 2004 -2008. To isolate the contribution of each of the specified source regions to the global dust loading, we tagged the dust originating from each source region.  Hourdin et al., 2013) with nudged monthly wind-fields from ERA-Interim (Dee et al., 2011), forced by sea-surface 215

Description of INCA simulations
temperatures and sea-ice concentrations and coupled to a land surface model (ORCHIDEE) (Krinner et al., 2005) and the aerosol and chemistry model INCA (INteraction with Chemisty and Aerosols) (Hauglustaine et al., 2004). The model covers both the troposphere and stratosphere with 79 hybrid-sigma layers, and with a spatial resolution of 2.5x1.25 degrees. A cross-comparisons between nudged and nonnudged simulations shows 15% larger emissions in the latter case, for aerosols that are emitted by 220 interactive surface-winds with a single coarse dust mode.
The aerosol size-distribution is represented by a multi-modal scheme which describes each aerosol species with one or more log-normal modes (Balkanski et al., 2004;Schulz et al., 2009). The INCA version used here implements black-carbon, nitrate, sulfate, ammonium, and organic aerosols with a combination of accumulation and coarse log-normal modes (both soluble and insoluble). In contrast, sea-225 salt is described by three soluble modes (accumulation, coarse and super-coarse). We use prescribed dimethyl sulfide (DMS), and secondary organic aerosols are not accounted for. Mineral dust is described by a new scheme (Di Biagio et al., 2020) based on 4 insoluble modes with following mass median diameters (MMD) at emission of 1.0, 2.5, 7.0 and 22.0 µm, and sigmas of 1.8, 2.0, 1.9, and 2.0 respectively. The prescribed size distribution at emission is partitioned among the four modes (0.57%, 230 4.2%, 30.8%, 62.4%), which ensures consistency with Kok et al. (2017) and measurements from the Fennec field campaign Experiment Ryder et al. (2013).
The different dust diagnostics X ( , , , ; ) estimated by INCA are given for each mode M. For the inter-model comparison based on size intervals, a post-processing of these diagnostics remap the results into a bins/sectional scheme, X ( , , , ; ) with upper bin limits at 2, 3.6, 6, 12, 20, 30, 40, 100 µm, 235 using a monthly time resolution. Of these bins, only the bins with upper limit up to 20 µm are used in this study. The size distribution at emission has prescribed invariant mass median diameters MMD, therefore the binning process, X(M)->X(B), is accurate for emitted fluxes. The concentrations are distributed amongst the binned intervals based on MMD ( , , , ; ) fields for the 4 modes, and the loadings have been estimated using these binned concentrations. Wet and dry deposition fluxes were binned starting 240 from the modal diagnostics using the MMD associated to the surface layer. This approximates the exact deposition fluxes which depend on the vertical size distribution which varies temporally and spatially. Therefore, in order to force mass conservation for each dust size bin, globally and seasonally-averaged correction factors c(B) were applied to both wet and dry deposition for each bin, such that X'(B) = c(B)X(B), and were X' is the effective wet or dry deposition flux. We used as constraint for these 245 correction factors that deposited mass per bin, integrated over the globe and the season, equals the total emitted mass per bin, also integrated over the globe and the season. These correction factors ranged between 0.5 and 1.5, and as such are relatively small compared to average error between model results and measurements of deposition fluxes (Figs. 7c, 10c). The methodology and code used to transform the modal description into the bin description is described in Checa-Garcia (2020). 250

Treatment of error correlation in the DAOD constraints
As described in the main text, we assume that the magnitude of errors that are completely random between seasons and regions ( � ,rand ), systematic errors that are correlated between different seasons for the same region ( � seas ), and systematic errors that are correlated across regions for a given season ( � ,reg ) are all equal. That is, we assume that 255 (S.1) Specifically, we obtain a realization of the seasonally-averaged DAOD for a given source region as follows ̅ = ̅ ,mean + �0, � ,rand ; = , � + �0, � seas ; = � + �0, � ,reg ; = �,

(S.2)
where ̅ ,mean denotes the mean DAOD for region p and season s (see Table 2). Further, (0, ; ) denotes a normal distribution with mean of 0 and standard deviation , which is sampled at its cumulative distribution at the point x, which denotes a randomly-drawn number between 0 and 1. The random 260 number , quantifies the point on the normal distribution of the random error and thus has different values for different seasons s and different DAOD-constrained regions p; the random number denotes the systematic error that stays constant for different seasons at a given region, and is thus equal for the different seasons but varies for different DAOD-constrained regions p; finally, the random number denotes the systematic error that stays constant for different regions for a given season, and is thus equal 265 for the different DAOD-constrained regions p but varies for different seasons s. Figure S1. Spatial distribution of annually-averaged bulk dust loading (g/m 2 ) produced by a unit (1 Tg) of total 270 atmospheric dust loading produced by each of the nine source regions, as simulated by each of the six models in our model ensemble. Figure S2. Inverse model results for the seasonally resolved PM20 dust emission flux.