Improving dust simulations in WRF-Chem model v4.1.3 coupled with GOCART aerosol module

In this paper, we rectify inconsistencies that emerge in the WRF-Chem v3.2 code when using the Goddard Chemistry Aerosol Radiation and Transport (GOCART) aerosol module. These inconsistencies have been reported, and corrections have been implemented in WRF-Chem v4.1.3. Here, we use a WRF-Chem experimental setup configured over the Middle East (ME) to estimate the effects of these inconsistencies. Firstly, we show that the old version underestimates the PM2.5 diagnostic output by 7% and overestimates PM10 by 5% in comparison with the corrected one. Secondly, we demonstrate that sub-micron dust 5 particles’ contribution was incorrectly accounted for in the calculation of optical properties. Therefore Aerosol Optical Depth (AOD) in the old version was 25-30% less than in the corrected one. Thirdly, we show that gravitational settling procedure, in comparison with the corrected version, caused higher dust column loadings by 4-6%, PM10 surface concentrations by 2-4%, and mass of the gravitationally settled dust by 5-10%. The cumulative effect of the found inconsistencies led to the significantly higher dust content in the atmosphere in comparison with the corrected WRF-Chem version. Our results explain why in many 10 WRF-Chem simulations PM10 concentrations were exaggerated. We present the methodology for calculating diagnostics we used to estimate the impacts of introduced code modifications. We share the developed Merra2BC interpolator, which allows processing Modern-Era Retrospective Analysis for Research and Applications, Version 2 (MERRA-2) output for constructing initial and boundary conditions for chemical species, and aerosols.

near-surface concentrations of particulate matter (PM), which comprise both PM 2.5 and PM 10 (particles with diameters less than 2.5 µm and 10 µm, respectively). Dust is the major contributor to PM over the ME region (Ukhov et al., 2020a). The ME is also subjected to the inflow of dust from the nearby Sahara Desert, which is another major dust source region Kalenderski and Stenchikov, 2016). Dust deposition fertilizes ocean surface waters and the seabed (Watson et al., 2000;Zhu et al., 1997). 5 Thus, given the impact of dust on climate, technology, human health, and ecosystems, an accurate description of dust effects in numerical models is essential. In the first place, it requires careful description of the dust cycle; from emission at the earth's surface, to transport in the atmosphere, and, finally, to removal by deposition.

5
The paper is organized as follows: Section 2 describes the WRF-Chem model setup. In Section 3, a description of the inconsistencies found in the WRF-Chem code and their effects on the results are presented. Conclusions are presented in Section 4.

WRF-Chem experimental setup
To quantify the effects of introduced code modifications, we use our typical model setup which we previously adopted for sim- fer Model (RRTMG) for both short-wave (ra_sw_physics=4) and long-wave (ra_lw_physics=4) radiation is used for radiative transfer calculations. Only the aerosol direct radiative effect is accounted for. More details on the physical parameterizations used can be found at http://www2.mmm.ucar.edu/wrf/users/phys_references.html.
The detailed description of all schemes is provided in LeGrand et al. (2019).
Here, we simulate dust emissions using the original GOCART-WRF scheme (dust_opt=1) proposed in Ginoux et al. (2001). 10 Dust emission mass flux, F p (µg m −2 s −1 ) in each dustbin p=1,2,...,5 is defined by the relation:  where, C (µg s 2 m −5 ) is a spatially uniform factor which controls the magnitude of dust emission flux; S is the source function (Ginoux et al., 2001) (see Fig. 1) that characterizes the spatial distribution of dust emissions; u 10m is the horizontal wind speed at 10 m; u t is the threshold velocity, which depends on particle size and surface wetness; s p is a fraction of dust emission mass flux within dustbin p.
Sea salt size distribution in the GOCART module is approximated by four sea-salt bins (see Tab. 2). Sea salt density is 2200 5 kg/m 3 . Emission of sea salt is calculated according to Gong (2003).

Dust emission tuning
To adjust to regional conditions, dust emission in the model is calibrated to fit observed AOD and aerosol volume size distributions (AVSD) obtained from the AErosol RObotic NETwork (AERONET; Holben et al. (1998)). AERONET AOD observations represent the total AOD with contributions from all types of aerosols. But because in the ME dust is more prevalent than all 10 other aerosols, we focus on dust emission only. More detailed information on dust emission tuning is provided in Ukhov et al. (2020a). For this study, we choose three AERONET sites: KAUST Campus, Mezaira, and Sede Boker located within the domain ( Fig. 1). We utilize level 2.0 (cloud screened and quality assured) AERONET AOD data. Note that from here onwards, we assume that AOD is given or calculated at 550 nm; see Appendix C.
2.1.1 Tuning the C parameter 15 To adjust dust emissions we first tune the C factor from Eq. 1, as practiced in our own studies (Kalenderski et al., 2013;Jish Prakash et al., 2015;Khan et al., 2015;Kalenderski and Stenchikov, 2016;Anisimov et al., 2017;Parajuli et al., 2019Parajuli et al., , 2020Ukhov et al., 2020a) and in the studies of other authors Kumar et al., 2014;Flaounas et al., 2017;Rizza et al., 2017). Our test runs indicate that for the ME, C = 0.5 achieves a good agreement between simulated and observed AOD at the KAUST Campus, Mezaira, and Sede Boker AERONET sites. Therefore this sub-optimal value C = 0.5 is retained 20 in all subsequent runs.

Tuning the s p fractions
We also tune s p fractions from Eq. 1 to better reproduce the AVSDs provided by AERONET retrievals using the spectral deconvolution algorithm (SDA) (O'Neill et al., 2003). The AERONET provides column integrated AVSD dV/dlnr (µm 3 /µm 2 ) on 22 logarithmically equidistant discrete points in the range of radii between 0.05 and 15 µm. For AVSDs we use the AERONET V3, level 2.0 product (Dubovik and King, 2000 As is the case with any partial differential equation solver, WRF-Chem requires the IC&BC for meteorological parameters and These profiles are obtained from the NALROM model (Liu et al., 1996)  to those used previously in WRF-Chem. To calculate the chemical IC&BC using MERRA-2 output, we develop an interpolator Merra2BC , which uses gaseous and aerosol fields from MERRA-2 reanalysis to construct the IC&BC required by the WRF-Chem simulation. For more details regarding the Merra2BC interpolator, see Appendix A.
To quantify the effect of each inconsistency we perform a WRF-Chem run where all the other corrections we discuss here are implemented, with the exception that we focus on a given time. The relative difference (%) of a specific set of variables in this run with respect to the ALL_OK run is presented as a measure of sensitivity to the chosen correction. All WRF-Chem runs are performed for 1-12 August, 2016. At the end of this section we estimate the cumulative effect of all inconsistencies. For this purpose we performed WRF-Chem simulation over the period from June 1 to December 31 of 2016. The subroutine sum_pm_gocart in module_gocart_aerosols.F calculates PM 2.5 and PM 10 surface concentrations using the following formulas:

10
where ρ is the dry air density (kg/m 3 ), DU ST 1,2,3,4 and SEAS 1,2,3 are the mixing ratios (µg/kg) of the dust in the first four bins and sea-salt in the first three bins, respectively. The contribution of the dust and sea salt bins to PM 2.5 and PM 10 is defined by the mapping coefficients d_25, d_10 for dust and s_25 for sea salt, see eq. 2. Black and organic carbon and sulfate also contribute to PM, but over the ME region their contributions are small in comparison to dust and sea salt, and we omit them for the sake of brevity. 15 We suspect that the default mapping coefficients are calculated incorrectly. Therefore we recalculated them assuming that dust and sea salt volume size distributions are functions of natural logarithm of particle radius. E.g., interpolation in the logarithm space is more accurate than in the radius space, as aerosol size distributions are smoother functions of logarithm than radius. The updated values of mapping coefficients s_25, d_25, d_10 along with their default values are presented in Tab. 3.
Effectively, the contributions in PM 2.5 of sea salt SEAS 2 decreases, while that of dust DU ST 2 increases. The contribution of 20 DU ST 2 in PM 10 decreases. The effects of using the updated mapping coefficients in place of default ones in PM calculation are shown in Fig. 2. We calculate the PM 2.5 and PM 10 concentrations in the lowest model layer using Eq. 2. Surface concentrations of dust and sea salt are computed using the procedure presented in Appendix E. With the default mapping coefficients, the model on average yields 7% lower PM 2.5 and 5% higher PM 10 concentrations over the ME.

Calculation of Aerosol optical Properties
For modeling in the ME, the treatment of optically active dust within the model is vitally important. AOD is calculated based on aerosol number-density and aerosol optical properties, which depend on the aerosol size and refractive index. In WRF-Chem, a parameterized Mie theory (Ghan and Zaveri, 2007) (Zaveri et al., 2008). Therefore, we further refer to them as MOSAIC bins (M OS 1,2,3,4,5,6,7,8 ).
To correctly calculate the dust optical properties we implement two corrections in the subroutine optical_prep_gocart() in 5 module_optical_averaging.F that computes the volume-averaged refractive index needed for Mie calculations.

Effect of Small Particles
In WRF-Chem's GOCART aerosol module, dust particle sizes span two orders of magnitude, from 0.1 to 10 µm; see Tab. 2.
However, we find that dust particles with radii between 0.1 and 0.46 µm are incorrectly accounted for in the Mie calculations of aerosol optical properties. Their mixing ratio is mapped on coarser MOSAIC bins than is required. Since finer particles have 10 a stronger effect on AOD per unit mass in comparison to coarser particles, the model AOD decreases. As a result, when tuning dust emission we push the model to emit more dust into the atmosphere, in order to fit the observed AOD. We rectify this error by correcting mapping fractions of DU ST 1 into MOSAIC bins, see Tab. 4. range and is therefore not accounted for in the mass redistribution. Also, the mass from DU ST 4 is only partially accounted for. After the changes, the dust mass from DU ST 1 bin is redistributed between finer M OS 3,4,5,6 bins compared to the original WRF-Chem where all DU ST 1 mixing ratio was mapped on the coarser M OS 5,6 bins.

Bin Concentration Interpolation
Originally, the subroutine optical_prep_gocart() redistributes dust and sea salt mass from GOCART into MOSAIC bins, using the assumption that dust size distribution is a function of particle radius. Consistent with Sect. 3.1, here we conduct interpolation assuming that dust distribution is a function of natural logarithm of radius. This modification causes changes in the mass redistribution between the GOCART and MOSAIC bins (see Tab. 5) and increases the contribution of small dust particles 5 into the AOD. Because the dust size distribution is a smoother function of the logarithm of a radius than the radius itself, interpolation is more accurate in logarithms than in radii. Table 5. Dust mass redistribution between GOCART and MOSAIC bins based a) on the assumption that bin concentration is a function of radius, and b) on the assumption that bin concentration is a function of natural logarithm radius.  To estimate the effect of these two corrections, we develop the WRF-Chem simulation NON_LOG_046, where only these two inconsistencies are not fixed, and compare the resulting AOD with that from the ALL_OK run. The AOD values are computed as described in Appendix C. As expected, the AOD increases after the corrections. Fig. 3 compares the AOD obtained from struggles to capture the observations. We find the worst correlation (R=0.42) and highest mean bias (-0.19) with AERONET AOD at the Mezaira station, which is located in a major dust source region (see Fig. 1). We obtain higher correlations with AERONET AOD of 0.66 and 0.75 for KAUST Campus and Sede Boker stations, respectively. Both of these stations are located outside the main dust source regions. Figure 4 shows the averaged AOD fields obtained from the ALL_OK and NON_LOG_046 runs, as well as their relative difference (%). We conclude that due to these two inconsistencies, averaged AOD obtained from the NON_LOG_046 run is lower by 25-30% on average over the ME in comparison with the ALL_OK run. Over Libya, Egypt, Oman, Iran, Azerbaijan, Turkmenistan, and Pakistan, the difference is even higher, reaching 30-35%. 15 We find that in the original WRF-Chem code the gravitational settling of dust and sea salt is calculated incorrectly. The default finite-difference scheme (implemented in the subroutine settling() file module_gocart_settling.F) does not account for change in air density when it calculates deposition mass flux. Thus, in the course of the gravitational settling the total mass of dust and The change of aerosol mixing ratio due to gravitational settling at downward directed velocity w is given by the following differential equation:

Gravitational Settling
where q is the aerosol mass mixing ratio (µg/kg) and ρ is the dry air density (kg/m 3 ). Using the first-order upwind scheme, this equation can be discretized into the following form: where ∆z k is the depth of the k model level, ∆t is the model time step. Subscript k denotes the model levels and superscript n 10 is the time-level. Taking into account that the calculation of gravitational settling is split from the calculation of the continuity equation, we assume ρ n+1 k ≈ ρ n k and get the following solution: Equation 5 is solved for each model column from the top to the bottom.
To validate the modified finite-difference scheme, we zero dust emissions across the whole domain, except for the 200 this area are included. We prohibit the inflow of dust from the domain boundaries by zeroing the corresponding boundary conditions, and we zero the initial dust concentrations to simplify calculation of the dust mass balance, which we compute using the following balance relation: Dust in the atmosphere = Emitted dust − (Grav. settled dust + Dry deposited dust) The amount of dust in the atmosphere is controlled by dust emission and dust deposition. The latter comprises gravitational 5 settling and dry deposition. For the sake of clarity, we refrain from introducing other dust removal processes, such as subgrid wet deposition (conv_tr_wetscav=0). The procedure of calculation of these diagnostics using the WRF-Chem output is provided in Appendix F.  calculating gravitational settling, as discussed above. This is in contrast with Fig. 5b, where we see perfect agreement between 15 the amounts of deposited dust plus dust in the atmosphere and emitted dust until the dust plume reaches the boundaries of the domain. Thus, this inconsistency in the gravitational settling subroutine is significant, as the error of 2.16% of total emitted mass accumulates within ≈20 hours. For a larger domain this imbalance will be more significant. This effect is especially important in the low-latitude desert regions. Zhang et al. (2015); Dipu et al. (2013); Huang et al. (2010) reported that in dry subtropics the boundary layer height can reach 6-7 km, which promotes the transport of dust particles to this altitude. When dust particles are settling from higher altitudes, a larger mass imbalance is accumulated.
We estimate the effect of the gravitational settling error by comparing averaged total dust column loadings (see Fig. 6a), 5 accumulated gravitationally settled dust (see Fig. 6b), and averaged dust and sea salt PM 10 surface concentrations (see Fig. 6c) obtained in ALL_OK and NOT_FIXED_GRAV_SETTLING runs, where the latter corresponds to the run with error in gravitational settling. We perform a comparison in terms of relative differences (%) in the runs, with and without corrections. Dust column loadings, gravitationally settled dust, and PM 10 surface concentrations are calculated according to the methodology described in Appendix D, F3, E, respectively. According to Fig. 6a,b,c, we observe higher negative values of relative difference 10 over non-dust source regions (see Fig. 1), i.e., over Sudan, Turkey, Yemen, Eritrea, Djibouti, and Ethiopia. By contrast, the relative differences over dust source regions, which include Egypt and the eastern part of Arabian Peninsula, are close to zero.
Coarse dust particles have shorter lifetimes in the atmosphere because of their higher deposition velocities. Thus, coarse dust particles are mostly deposited in the dust source regions, which explains close to zero values of relative difference in this region. Fine dust particles have longer atmospheric lifetime and thus can be transported over longer distances. The discrepancies 15 in the descriptions of the life cycle of fine dust explain larger relative errors in non-dust regions, as mentioned above.
Thus, we can conclude that in the NOT_FIXED_GRAV_SETTLING run, the total dust column loading is higher by 4-6% over the ME in comparison with the ALL_OK run. The computed total amount of dust in the atmosphere (see Appendix F4) was 6.41 and 6.72 Tg for ALL_OK and NOT_FIXED_GRAV_SETTLING runs, respectively. Hence, the amount of dust in the atmosphere is around 4.8% higher. The total amount of gravitationally settled dust is 5-10% higher on aver-20 age in NOT_FIXED_GRAV_SETTLING run. The biggest difference (15-25%) is observed in Sudan, Yemen, Eritrea, Djibouti, Ethiopia, and Turkey. The computed total amount of gravitationally settled dust (see Appendix F3) was 11 and 11.55 Tg for ALL_OK and NOT_FIXED_GRAV_SETTLING runs, respectively. Hence, the amount of gravitationally settled dust is around 5% higher in the NOT_FIXED_GRAV_SETTLING run. Dust and sea salt PM 10 surface concentrations (see eq. 2 and Appendix E) are higher by 2-4% on average over the ME in comparison with the ALL_OK run. We observe even bigger differences 25 (6-10%) over Eritrea, Djibouti, Ethiopia, and Turkey.

Case study
In the previous sections, we separately quantified the effect of each inconsistency in the WRF-Chem code and explained the associated physical links using short-term runs. In this section, we conduct a seven-month case-study to demonstrate the cumulative effect of all inconsistencies. We ran two WRF-Chem simulations from June 1 to December 31, 2016, using the  To adjust dust emissions in ALL_OLD run we tuned the C factor from Eq. 1. Our test runs indicated that C = 0.8 provides the best agreement between simulated and observed AOD. For the ALL_OK run we used C = 0.5 as before. Comparison of daily averaged AOD time-series obtained from the ALL_OK and ALL_OLD runs with the AERONET AODs at KAUST Campus, Mezaira, Sede Boker is presented in Fig. 7. AODs from both experiments are in good agreement with the AERONET AOD.
The Pearson correlation coefficients and mean biases (see Appendix B) with respect to AERONET AOD are in the ranges 5 of 0.62-0.75 and -0.03-0.07, correspondingly, for all AERONET sites. Thus, in ALL_OLD run the incorrect mapping of dust particles with radii between 0.1 and 0.46 µm causes stronger dust emissions in comparison with ALL_OK run. The stronger dust emissions lead to increased dust surface concentrations and to increased dust content in the atmosphere.  Figure   15 9 demonstrates the averaged over the summer (June, July, August) of 2016 total dust column loadings (g/m 2 ) and their relative differences (%) obtained from the ALL_OK and ALL_OLD runs. In some locations dust content in the atmosphere from the ALL_OLD run is higher by 80% in comparison with ALL_OK run. The total mass of dust in the atmosphere in the ALL_OK run yields 6.68 Tg in comparison with 10.92 Tg in the ALL_OLD run, so the difference exceeds 60%.

Effect of Initial and Boundary Conditions
We specifically conduct a sensitivity simulation to examine the impact of boundary conditions on PM 10 surface concentration over the ME. In this simulation boundary conditions are constructed using the developed Merra2BC interpolator (Ukhov 5 and Stenchikov, 2020) (see Appendix A) and we zero the initial concentrations of dust and sea salt. The emissions of dust and sea salt within the domain are turned off (dust_opt=0, seas_opt=0). In this instance, PM 10 concentrations are entirely determined by the inflow from the lateral boundaries. The averaged PM 10 surface concentrations are presented in Fig. 10.
PM 10 concentrations are calculated using Eq. 2. Figure 10 shows the inflow of PM 10 from Africa, Central Asia and from the Indian ocean. Dust is the major contributor to the PM 10 transported from Africa and Central Asia, whereas sea salt contributes 10 to PM 10 transported over the Indian ocean.

Conclusions
In this paper, we discuss the inconsistencies found in the WRF-Chem v3.2 model coupled with the GOCART aerosol module.
All of these inconsistencies are rectified in the WRF-Chem v4.1.3 code release. Here, we demonstrate the effect of the code rectification on WRF-Chem model performance. We also demonstrate the methodology we employ to calculate diagnostics, which we then use to estimate the effects of the changes made. To make these assessments, we configure the WRF-Chem domain over the ME and run it with 10 km grid resolution. The runs under discussion in this paper were performed over the period of 1-12 August, 2016. The effect of each inconsistency was estimated using specifically designed WRF-Chem runs where only one model inconsistency was activated.  The contribution of dust and sea salt into PM 10 surface concentration was also higher by 2-4% on average over the ME.
The cumulative effect of all inconsistencies was estimated in the seven-month case-study conducted for June 1 -December 31, 2016, when both AERONET AODs and PM 10 surface observations are available. The comparison of runs with and without proposed changes shows that the run without corrections yields higher dust loadings and total dust mass in the atmosphere by 80% and 60%, respectively. This seven-month case-study shows that the cumulative response to all code modifications applied simultaneously is stronger than the sum of their partial contributions. For instance, AOD underestimation causes higher dust emissions, which causes higher dust surface concentrations and increased production of dust in the atmosphere due to the error in gravitational settling. As a consequence, PM 10 surface concentrations further increases. Finally, an already high PM 10 surface concentration becomes even higher due to the incorrect calculation of PM 10 . Thus, the proposed improvements help 5 to explain the considerable bias towards higher PM 10 concentrations found in Flaounas et al., 2017;Su and Fung, 2015;Nabavi et al., 2017;Rizza et al., 2017;Eltahan et al., 2018).
In the course of improving the simulation of natural and anthropogenic aerosols and chemicals, we developed the capability to use MERRA-2 reanalysis for constructing WRF-Chem initial and boundary conditions for chemical species and aerosols.
The interpolation utility Merra2BC was coded for this purpose. Boundary conditions constructed using MERRA-2 reanalysis 10 more realistically account for the trans-boundary transport of aerosols. Merra2BC is made available to the community.
We believe the detailed quantification of the effects of the recent WRF-Chem code improvements are in line with open-source principles. The results of this work aim at better understanding of the model sensitivities to physical parameterizations. This work will add a greater understanding of model performance, and will be especially helpful for those who use the WRF-Chem model coupled with the GOCART aerosol module to carry out dust simulations over regions where dust plays an important role.
Code and data availability. The standard version of WRF-Chem is publicly available online at https : //github.com/wrf −model/W RF .
The full MERRA-2 reanalysis data set including aerosol and gaseous collections is publicly available online (https://disc. gsfc.nasa.gov/daac-bin/FTPSubset2.pl). Depending on the requirements, all or one of the following aerosol and gaseous collec-  Indexing for the vertical coordinate is from top to bottom, i.e., the first layer is the top layer of the atmosphere (P T OP ), while the 72nd layer is adjacent to the earth's surface.
In WRF-Chem, the pressure field is not given in wrf input_d01 and wrf bdy_d01 files. Hence, the pressure field must be restored using surface pressure P SF C taken from met_em_... * files created by metgrid.exe during the preprocessing stage.
Pressure at the top of the model wrf _p_top and η-values on half levels (znu) are taken from the wrf input_d01 file. The procedure of reconstructing the pressure from met_em_... * files using the python code is demonstrated in Fig. A1: Figure A1. A python script, which reconstructs the pressure using the met_em_... * files. nx, ny, nz -number of grid nodes in WRF-Chem domain.

A2 Mapping chemical species between MERRA-2 and WRF-Chem
Merra2BC file conf ig.py contains multiplication factors to convert MERRA-2 mass mixing ratios of gases given in kg/kg into ppmv. Aerosols are converted from kg/kg to ug/kg. When using the GOCART aerosol module in WRF-Chem simulation, 5 all MERRA-2 aerosols and gases are matched with those from WRF-Chem. We simply multiply by a factor of 10 9 to convert MERRA-2 aerosols mixing ratios given in kg/kg into ug/kg. In the case of gases, we need to multiply MERRA-2 mass mixing ratios by a ratio of molar masses M air /M gas multiplied by 10 6 to convert kg/kg into ppmv, where M gas and M air are molar masses (g/mol) of the required gas and air (28.97 g/mol), respectively. If another aerosol module is chosen in WRF-Chem, then different multiplication factors should be used.

A3 Interpolation procedure
A brief description of the interpolation procedure applied to the initial conditions is presented in Fig. A2.
For boundary conditions the procedure is similar, except that additional updates of domain boundary tendencies are required and interpolation is performed for each step, where boundary conditions are applied.

A4 Typical workflow 15
Here are the steps describing how to work with Merra2BC interpolator: Figure A2. Interpolation procedure applied to initial conditions. 4. Edit conf ig.py file which contains: (a) mapping of chemical species and aerosols between MERRA-2 and WRF-Chem; (b) paths to wrf input_d01, wrf bdy_d01, met_em . .. * files; (c) path to the downloaded MERRA-2 files; 5. real.exe sets default boundary and initial conditions for some chemical species. Merra2BC adds interpolated values to 5 the existing values, which may cause incorrect concentration values. To avoid this, run "python zero_f ileds.py", which will zero the required fields; 6. Run "python main.py", which will do the interpolation. As a result, files wrf input_d01, wrf bdy_d01 will be updated by the interpolated from MERRA-2 values; 7. Modify the WRF-Chem namelist.input file at section &chem: set have_bcs_chem = .true. to activate updated boundary conditions and, if needed, chem_in_opt = 1 to activate updated initial conditions; 8. Run wrf.exe.

Appendix B: Statistics
The following statistical parameters were used to quantify the level of agreement between estimations and observations:

5
Pearson correlation coefficient (R): Mean bias (bias): O i their averages and N is the 10 number of data.

Appendix C: AOD calculations
WRF-Chem does not calculate AOD at 550 nm (only at 300, 400, 600, 1000 nm variables TAUAER1, TAUAER2, TAUAER3, TAUAER4, respectively), but, instead, it outputs the extinction coefficient at 550 nm (variable EXTCOF55). The AOD at 550 nm (AOD 550 ) for (i, j) vertical column can be calculated by summing throughout the vertical column of product of 15 multiplication of the EXTCOF55 by the ∆z: where ∆z i,j,k is the depth (m) of the (i, j, k) cell, which can be computed using the formula: where P H is the geopotential and P HB is the perturbed geopotential and g=9.81 m/s 2 is the gravitational acceleration. 20 Variables P H and P HB are taken from the WRF-Chem output.
To facilitate comparison with the model output the 550 nm AOD is calculated using the following relation: where α is the Ångström exponent for the 440-675 nm wavelength range provided by AERONET, τ λ is the optical thickness at wavelength λ, and τ λ0 is the optical thickness at the reference wavelength λ 0 .

25
Appendix D: Column loadings WRF-Chem stores dust column loadings (µg/m 2 ) using variables DU ST LOAD_1, 2, 3, 4, 5. Column loadings for (i, j) vertical column of other aerosols or chemical species can be computed by vertically summing throughout the vertical column of product of multiplication of the mass mixing ratio q (µg/kg) by the cell depth ∆z (m) (see eq. C2) and dry air density (kg/m 3 ).
WRF outputs variable ALT , which is inverse dry air density (m 3 /kg): WRF-Chem outputs gases concentrations expressed in ppmv. Conversion from ppmv into the mass mixing ratio can be calculated using the following formula: M ass mixing ratio = ppmv · 10 −6 · M gas /M air , where M gas and M air are molar masses (g/mol) of the required gas and air (28.97 g/mol), respectively.

Appendix E: Surface concentrations
Surface concentration (µg/m 3 ) of an aerosol at (i, j) vertical column can be computed by multiplication of the mass mixing ratio (µg/kg) at the first model level (q 1 ) by the corresponding dry air density (kg/m 3 ) at the first model level (1/ALT 1 ): Surf ace concentration i,j = q i,j,1 · 1/ALT i,j,1 To obtain gas surface concentration (µg/m 3 ), (ppmv) needs to be converted to the mass mixing ratio; see eq. D2.

Appendix F: Dust mass balance
In the WRF-Chem's GOCART aerosol module, dust emissions along with three types of removal processes (dry deposition, gravitational settling, and wet scavenging) are implemented. Here, for the sake of clarity, we refrain from consideration of wet scavenging. To calculate the dust mass balance, assuming there is no flow of dust through the domain boundaries, we need to calculate the amount of dust emitted from the domain area, the amount of dust that was deposited by gravitational settling and 20 dry deposition, and the amount of dust in the atmosphere. By default, WRF-Chem stores instantaneous values of dust emission and deposition fluxes. We modified the WRF-Chem code to accumulate the dust emission and deposition fluxes.

F1 Grid column area
In WRF, one of the following four projections can be used: the Lambert conformal, polar stereographic, Mercator, and latitudelongitude projections. These projections are implemented using map factors. In the computational space, the grid lengths ∆x between grid points vary with position on the grid. Map factors mx i,j and my i,j for both the x and y components are used for the transformation from computational to physical space, and computed by geogrid.exe during the preprocessing stage. mx i,j and my i,j are defined as the ratio of the distance in computational space to the corresponding distance on the earth's surface (Skamarock et al., 2008): Map factors mx i,j and my i,j for each (i, j) vertical column are stored in wrf input_d01 file in variables M AP F AC_M X and M AP F AC_M Y , respectively. Thus, the area of (i, j) column S i,j (m 2 ) in physical space is calculated using formula: where S i,j is the area of the (i, j) column (m 2 ); see eq. F2. Here we divide S i,j by ∆x · ∆y to account for the fact that in the subroutine gocart_dust_driver() dust emission are calculated in the computational space where grid cells have dimensions ∆x and ∆y.
F3 Gravitational settling and dry deposition 20 The subroutines settling() implemented in module_gocart_settling.F and gocart_drydep_driver() implemented in module_gocart_drydep.F are used to calculate gravitational settling and dry deposition of dust. By default, instantaneous gravitational and dry deposition fluxes (µg/m 2 s) are stored in variables GRASET _1, 2, 3, 4, 5 and DRY DEP _1, 2, 3, 4, 5, respectively. Thus, multiplying these fluxes on each timestep by the timestep ∆t and the scaling coefficient 10 −9 , and by adding the resulting value to the previous value, we obtain accumulated gravitational and dry deposition mass per unit area expressed in (kg/m 2 ). 25 Hence, deposition of the dust from the first dust bin due to gravitational settling (Grav. settled dust 1 , kg) and dry deposition (Dry. deposited dust 1 , kg) is calculated using the following formulas: Grav. settled dust 1 = i,j where S i,j is the area of the (i, j) column (m 2 ); see eq. F2.

F4 Dust in the atmosphere
There are two approaches to calculate the amount of dust in the atmosphere (Dust in the atmosphere, kg). In the first approach we use dust column loadings (variables DUSTLOAD_1,2,3,4,5, µg/m 2 ). Thus, the mass of dust in the first dust bin 5 is given: Dust in the atmosphere 1 = 10 −9 · i,j where S i,j is the area of the (i, j) column (m 2 ); see eq. F2.
In the second approach we calculate the mass of air in each grid cell, multiply it by the dust mass mixing ratio (for example DU ST 1 , µg/kg), and sum over all grid cells in the domain: 10 Dust in the atmosphere 1 = 10 −9 · i,j S i,j · k DU ST 1 i,j,k · ∆z i,j,k · 1/ALT i,j,k , where ∆z i,j,k is the depth (m) (see eq. C2) and ALT i,j,k is the inverse dry air density (m 3 /kg) in the grid cell (i, j, k).
Gaseous concentrations expressed in ppmv need to be converted into mass mixing ratios (µg/kg); see eq. D2.
Author contributions. A. Ukhov planned and performed the calculations, wrote the manuscript, and led the discussion. R. Ahmadov, G.
Grell, and G. Stenchikov participated in the discussion and reviewed the manuscript.