Simulating Atmospheric Dust With a Global Variable‐Resolution Model: Model Description and Impacts of Mesh Refinement

In this study, a global variable‐resolution modeling framework of atmospheric dust and its radiative feedback is established and evaluated. In this model, atmospheric dust is simulated simultaneously with meteorological fields, and dust‐radiation interactions are included. Five configurations of global mesh with refinement at different resolutions and over different regions are used to explore the impacts of regional refinement on modeling dust lifecycle at regional and global scales. The model reasonably produces the overall magnitudes and spatial variabilities of global dust metrics such as surface mass concentration, deposition, aerosol optical depth, and radiative forcing compared to observations and previous modeling results. Two global variable‐resolution simulations with mesh refinement over major deserts of North Africa (V16 km‐NA) and East Asia (V16 km‐EA) simulate less dust emissions and smaller dry deposition rates inside the refined regions due to the weakened near‐surface wind speed caused by better resolved topographic complexity at higher resolution. The dust mass loadings over North Africa are close to each other between V16 km‐NA and the quasi‐uniform resolution (∼120 km) (U120 km), while over East Asia, V16 km‐EA simulates higher dust mass loading. Over the non‐refined areas with the same resolution, the difference between global variable‐resolution and uniform‐resolution experiments also exists, which is partly related to their difference in dynamic time‐step and the coefficient for horizontal diffusion. Refinement at convection‐permitting resolution around the Tibetan Plateau (TP) simulates less dust due to its more efficient wet scavenging from resolved convective precipitation around the TP against coarse resolution.

Although limited-area simulation at relatively high resolution can generally well characterize dust and its impacts on a regional scale, lateral boundary conditions may introduce some numerical issues due to their inconsistency with simulation and constrain regional feedback to large-scale circulation, such as interactions of dust-cloud and dust-radiation (e.g., Giorgi & Mearns, 1999;Laprise et al., 2008;Leung et al., 2013;Prein et al., 2015; Y. Q. Wang et al., 2004).In addition, the limited-area model cannot be used for forecasting alone because it needs a global model to provide lateral boundary conditions.In recent years, global variable-resolution atmospheric models with high-resolution regional refinement have been used in various research topics.These models exhibit the capacity to avoid the challenges posed by lateral boundary conditions and regional feedback issues, while concurrently achieving a reduction in computational expenses when compared to globally uniform high-resolution models (e.g., Zhao et al., 2016Zhao et al., , 2019)).Most of global variable-resolution models have been applied in studying regional weather and climate characteristics (e.g., Burakowski et al., 2019;Gettelman et al., 2018;X. Y. Huang et al., 2016;Michaelis et al., 2019;Rauscher et al., 2013;Rhoades et al., 2016Rhoades et al., , 2018;;Xu et al., 2021;Zarzycki et al., 2014Zarzycki et al., , 2015;;Zhao et al., 2016Zhao et al., , 2019)).
Several studies have also attempted to simulate regional aerosols and their impacts with global variable-resolution atmospheric models (e.g., Rahimi et al., 2019;Wu et al., 2018).Wu et al. (2018) used the variable-resolution Community Earth System Model (VR-CESM) with a regionally refined high-resolution of 0.125° or ∼14 km to examined the impacts of black carbon and dust on snow melting, the hydrologic cycle, and surface air temperature in the Rocky Mountains.Rahimi et al. (2019) used the VR-CESM with refined resolution of 0.125° (∼14 km) and found that black carbon and dust led to an enhancement of the South Asian monsoon through radiation-dynamic feedback that enhanced precipitation in May and June.These global variable-resolution models with aerosol effects in previous studies are based on climate models and are hydrostatic models.The hydrostatic approximation could be considered a reliable and accurate assumption for meso-scale and global-scale weather systems.However, for relatively small scale phenomena where the vertical and horizontal scales of motion are comparable, departures from the hydrostatic balance could be significant (Yang et al., 2017).Therefore, hydrostatic models are not suitable for regional refinement at a horizontal resolution of a few kilometers.However, regional refinement at a few kilometers (so-called convection-permitting scale) may be needed to study aerosol impacts, which can properly resolve aerosol heterogeneous spatial distributions and aerosol-cloud-radiation-precipitation interactions, particularly over regions with complex topography.
To date, very few studies have attempted to simulate and investigate dust and its impacts with a non-hydrostatic global variable-resolution model.Lau et al. (2020) simulated dust and its climatic impacts with a non-hydrostatic global variable-resolution dynamic core implemented in the CESM (Model for Prediction Across Scale [MPAS]-Community Atmosphere Model version 5 [CAM5]) (Zhao et al., 2016).Their modeling framework is designed for simulating dust climatic impacts instead of simulating specific dust events, for example, how dust plays a role in weather and air quality forecasting.In addition, their modeling framework does not have the capability for mesh refinement to the convection-permitting scale (a few kilometers) due to some model infrastructure limits, which is needed for investigating dust impacts on convective systems.Finally, their modeling framework uses the modal size representation for dust.Zhao et al. (2013a) found that models with modal size FENG ET AL.

10.1029/2023MS003636
3 of 32 representation have relatively large biases in simulating dust mass loading and radiative forcing, even with constrained AOD by observations.They found a difference of 25% in dust mass loading between the simulations using modal and sectional size representations.Therefore, in this study, based on a global variable-resolution non-hydrostatic atmospheric model, MPAS (Model for Prediction Across Scale, v7.0) (Du et al., 1999;Ringler et al., 2008;Skamarock et al., 2012), a sectional aerosol modeling framework is first time established and evaluated.The global variable-resolution simulation of dust and its radiative forcing with regional refinement to the convection-permitting scale is also realized for the first time.
The modeling sensitivities of dust and its radiative forcing to mesh refinement are also investigated because model spatial resolution is a key factor that can significantly affect the simulated regional and global distributions of dust and its impacts on weather, air quality, and climate systems (e.g., Y. Feng et al., 2022).Global variable-resolution dust modeling may better simulate dust impacts on regional weather and climate with regional refinement because some key meteorological processes, such as clouds and precipitation, can be better simulated at higher horizontal resolutions (e.g., Lau et al., 2020;Michaelis et al., 2019;Sakaguchi et al., 2015;Xu et al., 2021;Zhao et al., 2019).In addition, the complex topography can also be better resolved over mountainous regions such as the Himalayas, which is important for modeling the transport of air particles and moisture on a regional scale (e.g., Li et al., 2022;Zhang et al., 2020).This paper is organized as follows.Section 2.1 briefly introduces the MPAS-atmosphere, and Section 2.2 describes the dust modeling framework in detail.The model configuration and data used for this study are described in Sections 2.3 and 2.4, respectively.The evaluation of the new model is shown in Section 3.1.The series of global uniform-and variable-resolution experiments are analyzed in Section 3.2.The conclusion and summary are in Section 4.

MPAS-A
The model developed in this study is based on the framework of the Model for Prediction Across Scales for the atmosphere (MPAS-A, v7.0).MPAS-A is a global variable-resolution atmospheric model that uses unstructured spherical centroidal Voronoi tessellations (Du et al., 1999;Ringler et al., 2008;Skamarock et al., 2012).One can use a single scalar density function to generate global variable-resolution meshes with regional refinement in MPAS-A (Ju et al., 2011).The atmospheric solver in MPAS-A (Skamarock et al., 2012) integrates the non-hydrostatic equations, and as such, it is suitable for both weather and climate simulation, that is, for both non-hydrostatic and hydrostatic flow simulation.MPAS-A uses an unstructured centroidal Voronoi mesh (grid, or tessellation) and C-grid staggering of the state variables as the basis for the horizontal discretization in the fluid-flow solver.There are a few physics packages in the current version of MPAS-A (v7.0).Since this study focuses on the development and application of dust modeling schemes in MPAS-A, more details about physics schemes in MPAS-A can be found in previous studies (e.g., Xu et al., 2021;Zhao et al., 2019).Recently, the MPAS-A framework has been widely used to address scientific questions related to clouds and precipitation (e.g., O'Brien et al., 2013;Xu et al., 2021;Zhao et al., 2016Zhao et al., , 2019)), the structure of the intertropical convergence zone (e.g., Landu et al., 2014), atmospheric river frequency (e.g., Hagos et al., 2015), eddy-driven jets (e.g., Lu et al., 2015), regional climate (e.g., Sakaguchi et al., 2015Sakaguchi et al., , 2016)), tropical cyclones (e.g., Davis et al., 2016), and global atmospheric predictability at the convection-permitting scale (e.g., Judt, 2018).Therefore, MPAS-A has been proved to be able to provide reliable meteorological fields, such as wind and precipitation, which are required for modeling dust.

Dust Modeling Schemes in MPAS-A
In this study, a sectional aerosol modeling framework is established based on MPAS-A.Dust is one aerosol species that can be simulated in this modeling framework.To be focused, in this study, dust is the only aerosol species simulated and all other aerosol species are excluded in the experiments.All physical processes related to dust modeling are described below in detail.

Dust Emission
The dust emission scheme implemented in the model is based on the method of Ginoux et al. (2001) and Zhao et al. (2010).The dust emission flux where C is an empirical proportionality constant and is highly tunable, which is used to make the left and right sides of the equation have the same dimension.In this study, we follow Zhao et al. (2010) to set C as 1.0 × 10 −9 kg • s 2 • m −5 , which is reasonable as shown in the discussion later that the simulated global mean dust AOD is consistent with the observational estimate.S is a source function that defines the fraction of alluvium available for wind erosion.The source function S, shown in Figure 1, is prescribed as Ginoux et al. (2001).The original S has a spatial resolution of 1° × 1°.It is interpolated into different horizontal resolutions as needed using the "bilinear" interpolation method.s p is the mass fraction of dust particles at each size during the emission process.The total dust emission flux is the sum of the dust emission flux of each size.This study takes the assumption in Ginoux et al. (2001) that (a) the mass fraction of dust particles with diameter smaller than 2 μm equals 1/10 of the total mass of emitted dust particles with diameter larger than 2 μm and (b) the mass fraction of each subclass for dust particles with diameter larger than 2 μm is the same, that is, the s p values are 0.1 for the size 0.2-2 μm and 1/3 for the sizes 2-3.6 μm, 3.6-6 μm, and 6-12 μm in diameter.Zhao et al. (2010) added another subclass with diameters of 12-20 μm, and the s p values were then set as 1/4 for each size of dust particles with diameters larger than 2 μm.This study further extends the diameter range of emitted dust particles up to 40 μm, and the s p values are thus 0.1 for the size 0.2-2 μm and 1/5 for the sizes 2-3.6 μm, 3.6-6 μm, 6-12 μm 12-20 μm, and 20-40 μm.Please note that the formula listed above and the values of s p are only for calculating the total dust emission flux.The mass of emitted dust will be redistributed following the dust size distribution as introduced in Section 2.2.2.In addition, there are very few studies and measurements of dust particles larger than 40 μm in diameter.Therefore, although emitted dust particles with diameters of 20-40 μm are included experimentally in the model, this study focuses on the analysis of dust with diameters less than 10 μm, as in most previous models.More observations are needed to constrain and optimize the modeling processes associated with particles with diameters larger than 10 μm in the future.
u 10m is the wind velocity at 10 m, and u ts is the threshold wind velocity.If u 10m is smaller than u ts , dust emissions will not occur.The threshold wind velocity u ts is a function of the density and size of dust particles, air density, and surface moisture, which is calculated as where w is the relative wetness of soil that is calculated through dividing soil moisture by the saturated value of soil moisture (the saturated value of soil moisture is calculated in the land-surface scheme (F.Chen & Dudhia, 2001)), and u ts0 is the threshold velocity for dry soil (while w = 0.1), which is calculated following Marticorena and Bergametti (1995): where ρ dust is the density of soil particle that is 2.5 g/m 3 for the particle with diameter smaller than 2 μm and 2.65 g/m 3 for others following Ginoux et al. (2001), ρ air is the density of air, d is the diameter of dust particles, and B is a dimensionless number called the friction Reynolds number with an approximate expression as  =   dust +  , in which a = 1,331 cm −x , b = 0.38, and x = 1.56.It is worth noting that this formula for u ts0 is only valid for the range of 0.03 < B < 10.In this study, the diameters of dust particles range from 0.04 to 40 μm, so the values of B range from 0.38 to 0.622.Therefore, this formula holds in this study.The formula for larger dust particles can be found in Marticorena and Bergametti (1995) and is not considered in this study.

Size Parameterization
With total dust emission, the emitted dust particles are redistributed following the log-normal size distribution as follows: where N total is the total number concentration, d is the particle diameter, d n is the number median diameter, and σ n is the standard deviation of the number distribution.Equation 4 can be written as n(ln d)/N total ∼ N(ln d n , ln σ n ), where N(μ, σ) is the normal distribution function.Combining  =  6  3 and Equation 4, the volume distribution of dust particles is Consider d 3 = exp(3 ln d); then, so: where and ln d v = 3 ln 2 σ n + ln d n .Equation 7 can be written as v(ln d)/V total ∼ N(ln d v , ln σ n ).Following Zhao et al. (2010), two log-normal distributions are used to fit dust size distributions.In emissions, 15% of the total emitted dust mass follows the first log-normal distribution, and 85% follows the second one.The specific parameters for the two log-normal distributions were suggested by Zhao et al. (2010) with d v = 2.91 μm and σ n = 2.20 for the first one and d v = 6.91 μm and σ n = 1.73 for the second one.Following most global models, the size distributions of emitted dust are treated consistently globally because there is still insufficient information to treat them differently globally, which may affect the modeling results of dust and deserves further investigation in the future.
In this study, the sectional approach is selected to represent the size distributions of dust particles because Zhao et al. (2013a) found that the sectional approach with enough bins is more accurate than the modal approach with two parameters (mean and standard deviation) in simulating the evolution of dust particles in the atmosphere.Some previous studies have used the sectional approach for dust modeling, and their dust diam-eters generally range from nanometers to 10 μm (e.g., Ginoux et al., 2001;Mahowald et al., 2006;Takemura et al., 2009;Zender et al., 2003;Zhao et al., 2013a).In addition, most of them represented dust size distribution with a small number of bins (e.g., 4 and 5), which may result in some biases in number concentration and radiative forcing (Zhao et al., 2013a).Therefore, in this study, 8 bins (PSD8) are used to represent the dust size distribution from approximately 0.04-10 μm, as listed in Table 1 (Bin 1∼Bin 8), for simulating the dust life cycle and radiative forcing.Given the number of bins (N) and the upper and lower bounds (i.e., D U and D L ) of the entire size range (i.e., 10 μm and 40 nm in this study), the upper and lower size bounds of the ith bin (i.e., d U (i) and d L (i)) are calculated following Zhao et al. (2010): Some studies have emphasized the importance of dust with diameters up to 20 μm (e.g., Adebiyi & Kok, 2020;Kok, 2011;Van der Does et al., 2018).Therefore, to establish the capability of studying coarse dust particles in the future, the upper bound of the diameter in the model is further extended to 40 μm by adding two more bins ranging from 10 to 20 μm (ninth bin) and 20-40 μm (tenth bin) (Table 1).Please note that there are very few studies on dust particles larger than 10 μm in diameter using global models, and the physical processes for these large dust particles globally, such as emission and dry deposition, have not yet been fully assessed, which need to be constrained and optimized with Therefore, the ninth and tenth bins (10-40 μm) are only implemented for future model development and are not focused on in this study.Each size bin is assumed to be internally mixed so that all particles within a size bin have the same properties.
Given the lower and upper bounds of the diameter and the log-normal distribution, the volume fraction of each size bin of emitted dust can be calculated as: where v is the dust volume fraction of each size bin, d L and d U are the lower and upper bounds of the diameter, respectively, and erf(x) is the "error function" encountered in integrating normal distribution, which is defined as The volume fractions for each size bin are listed in Table 1.Please note that following the two log-normal distributions described above, the sum of the dust volume fraction for each size bin is close to 1 for 10 bins and is less than 1 for other size ranges (98% for 9 bins, 80% for 8 bins).The volume size distributions of emitted dust are demonstrated in Figure 2. The dust particles within each size bin are assumed to be uniformly distributed in this study, and the number of emitted dust particles within each size bin can be estimated with the effective diameter of each size bin (   = √  () * () ) and the spherical assumption of dust particles.

Dust Transport
All transport of dust aerosol for different sizes is performed online, including up to 30 aerosol-related variables that are fully prognostic, such as dust mass and number.The transport scheme of dust aerosols follows the same rules as the scheme of moisture scalars in MPAS-A as described in Skamarock and Gassmann (2011).The scheme employs a third-order Runge-Kutta method on the irregular Voronoi (hexagonal) meshes of MPAS-A to solve the fully  1 for emitted dust.The red curve represents the distribution described by Equation 7. The markers represent the volume fractions for each size bin.
compressible non-hydrostatic equations of motion.The prognostic equations integrated are cast in conservative form for conserved variables.These equations are discretized in a finite volume formulation, and as a result, the model exactly conserves mass.The turbulent transport of dust particles within the planetary boundary layer (PBL) is parameterized based on the exchange coefficient that is diagnosed from the PBL schemes.In this study, cumulus convective transport is not included.Although cumulus convection activity is often relatively weak over dry areas with dust sources, it may affect the long-range transport of dust over remote regions, which deserves further development and investigation in the future.

Dry Deposition and Gravitational Settling
First, dry deposition is parameterized based on Peters and Eiden (1992) in this model and is active at the bottom vertical layer of model.The formula for dry deposition velocity is given by In this function, R a is the aerodynamic resistance and R s is the surface resistance.The aerodynamic resistance R a is calculated following the method described in Wesely and Hicks (1977).
where k is the von Karman constant (approximately 0.41), u * is the friction velocity, and ψ c is a correction for atmospheric diabatic effects (see Wesely and Hicks (1977) for details).The surface resistance is calculated as: where u * is the friction velocity.The efficiency term for diffusion E dif is given by: where S c is the Schmidt number, calculated by S c = υ/D, υ is the kinematic viscosity of air, and D is the molecular diffusivity, which is calculated using the Stokes-Einstein equation: where k B is the Boltzmann constant and T is the temperature.The efficiency terms for impaction E imp and interception E int are expressed as: where R′ is the rebound correction factor used in Slinn (1982).
z r is the roughness length, and S t is the Stokes number, which is given by where d c is the collector diameter of the obstacle with a value of 0.002 in this study.
The gravitational settling velocity is calculated as: where μ is the dynamic viscosity of air, d is the particle diameter, and C c is the Cunningham correction factor that follows the method in Pruppacher and Klett (1978): C c = 1.0 + (1.257 + 0.4e −1.1d/2λ ) (2λ/d), λ is the mean free path of air.Gravitational settling not only removes dust particles from the bottom vertical layer of the model but also brings particles from the vertical layers above to below in the atmosphere.
At every vertical layer, the adjustment of mass concentrations of dust particles within each size bin due to dry deposition and gravitational settling is calculated as: , Bottom layer of model , Other vertical layers (20) where q dust (k) is the concentration of dust particles in the kth layer at time step t,   ′ dust () is at time step (t + ∆t), ∆z(k) is the height of the kth layer, and ∆t is the model time step.

Aerosol Cloud Activation and Wet Scavenging
Aerosol-cloud interaction is implemented in the model based on the method described by Gustafson et al. (2007) for calculating the activation and resuspension between dry aerosols and cloud droplets.Aerosol activation (or droplet nucleation) is based on a maximum supersaturation determined from a Gaussian spectrum of updraft velocities, similar to the methodology used in Ghan et al. (2001).The activation droplet number is coupled with the Thompson microphysics scheme (Thompson et al., 2008).In this way, aerosols can affect cloud droplet number, and clouds can also alter aerosol concentration through aqueous processes and wet scavenging.Within clouds, aerosol particles are activated into the cloud water phase in each size bin, for example, represented by Dust cw01 for size bin one.Dust particles can also be considered as IN.Within the Thompson cloud microphysics scheme, the number of IN in mixing-phase clouds from dust (N i ) is calculated following the formula proposed by DeMott et al. (2010): where a = 5.94 × 10 −5 , b = 3.33, c = 0.0264, d = 0.0033, T 0 = 273.15K, n is the number concentration of dust particles in each size bin with effective diameter (as described in Section 2.2.2) larger than 0.5 μm and T is the temperature in K.It should be noted that the calculated N i is an overall number without size distributions.
This study only considers the wet scavenging process of activated dust aerosols into cloud droplet, ignoring the conversion of dust into IN.The wet scavenging process of dust particles includes in-cloud (rainout) and below-cloud scavenging.Wet scavenging related to parameterized cumulus convection is not considered in this version of the model and deserves further development.Within clouds, dust particles that are in the cloud water phase are collected by rain, graupel, and snow.The mass concentrations of Dust cw in each size bin are removed by the same loss rate following the formula: where dt is the simulation time step and qlsink is the loss rate of cloud water calculated by the Thompson cloud microphysics scheme (Thompson et al., 2008) as follows: where q represents concentration; r is for rain, i is for cloud ice, s is for snow, w is for cloud water, and g is for graupel.The subscripts represent different processes: au is for autoconversion, fz is for freezing and cw is for collection by cloud water.For example, qr wau means the changing concentration of cloud water into rain through the process of autoconversion.qc represents the total concentration of cloud water.Below clouds, wet scavenging follows Slinn (1984).The scavenging rate is calculated as follows: where r and p represent the radius (in cm) of rain drops and particles, respectively; ϵ and E are the hydrometeor retention and collision efficiencies, respectively; the product ϵE is the collection efficiency; v t and v g are the terminal speed of rain drops and gravitational settling speed of the particles, respectively; N is the number density of rain drops; and A x = π(r + p) 2 is the effective cross-sectional area.To simplify the calculation, here, the gravitational settling speed of the particles is ignored, and the retention efficiency is taken to be unity (ϵ = 1).The radius range of rain drops from 0.005 to 0.25 cm is considered.The final formula for the scavenging rate is as follows: The particle collision efficiency is calculated as: where R e is the Reynolds number, S c is the Schmidt number, κ = p/r, f = v I (max)/v t , and v I is the internal circulation.S t is the Stokes number in Equation 18 with v t replacing u * and is the critical Stokes number.In addition, convective transport and wet scavenging of dust by cumulus clouds are parameterized based on the schemes of H. Liu et al. (2001) and Zhao et al. (2009Zhao et al. ( , 2013a)).

Optical Properties
Aerosol optical properties such as extinction, single scattering albedo, and asymmetry factor for scattering are computed as a function of wavelength and three-dimensional position.As the computation of Mie scattering requires large computing resources, to compute the extinction efficiency Q e and the scattering efficiency Q s more efficiently, the methodology described by Ghan et al. (2001) is applied in this study.In this way, the full Mie calculation is performed only once to obtain a table of seven sets of Chebyshev expansion coefficients, and later, the full Mie calculations are skipped, and Q e and Q s are calculated using bilinear interpolation over the Chebyshev coefficients stored in the table.The detailed method of the computation of aerosol optical properties in the model is similar to the description in Barnard et al. (2010), Fast et al. (2006), andZhao et al. (2013b).The Optical Properties of Aerosols and Clouds data set (Hess et al., 1998) is used for the shortwave (SW) and longwave (LW) refractive indices of dust aerosols.The radiative feedback of dust is coupled with the rapid radiative transfer model (RRTMG) (Iacono et al., 2000;Mlawer et al., 1997) for both LW and SW radiation, as in Zhao et al. (2011).The AOD and direct radiative forcing of dust aerosols are diagnosed following the methodology by Zhao et al. (2013b).In this methodology, the calculation of aerosol optical properties and the radiative transfer scheme is performed twice with and without considering the masses of dust and its associated water aerosol.
After this diagnostic procedure, the optical properties (e.g., AOD) and direct radiative effects for dust can be estimated by subtracting the optical properties and direct radiative effects calculated by the without-dust procedure from those calculated by the with-dust procedure.It can be described as:

Numerical Experiments
In this study, the model is configured with 55 vertical layers, and the model top is set at 30 km.Please note that, in this version of model, the vertical grid cannot be nested.Five sets of experiments are conducted with different mesh structures: two quasi-uniform-resolution meshes and three variable-resolution meshes.The meshes are shown in Figure 3.The quasi-uniform mesh has essentially the same mesh spacing globally, while the variable-resolution meshes have finer mesh spacing in the refined regions with a transition zone between the fineand coarse-resolution meshes.More details about the mesh generation method can be found in Ju et al. (2011).
In this study, the quasi-uniform meshes (U120 km and U60 km) have grid spacings of approximately 120 and 60 km globally.The first variable-resolution mesh features a circular refined high-resolution region with 16 km grid spacing centered over East Asia (hereafter referred to as V16 km-EA), which covers the major dust source region of East Asia and the mountainous region with complex topography around the Tibetan Plateau (TP) (Figure 3c).The resolution outside this region is approximately 120 km (Figure 3c).This region is also largely influenced by the Asian monsoon.Lau et al. ( 2020) highlighted the impacts of dust-cloud-radiation-precipita- tion interactions on the Asian summer monsoon over this region and found that the interactions and impacts can be better simulated at higher horizontal resolutions.The complex topography over this region needs to be resolved with high resolution to better simulate moisture and aerosol transport toward the TP (e.g., Li et al., 2022;Zhang et al., 2020).The configuration of this global variable-resolution simulation lays the foundation for further investigation of dust impacts on the Asian monsoon over the TP region.The second variable-resolution configuration has a similar mesh to the first one, but the circular refined high-resolution region is centered over North Africa (referred to as V16 km-NA hereafter).The refined regions cover the largest dust source area of the globe.The topography over this region is less complex than that over East Asia.These two configurations may help to investigate the impacts of refinement on dust modeling over regions with different topographic complexity.One more global variable-resolution mesh features a circular refined region at 4 km resolution centered around the TP (referred to as V4 km hereafter) to cover the areas with East Asian dust sources and complex topography (Figure 3g).The resolution outside this region is approximately 60 km.This configuration is also used to demonstrate the capability for modeling dust and its impacts at the convection-permitting scale, one unique feature of this modeling framework.The detailed information of five meshes of the experiments in this study is summarized in Table 2.
Three experiments with the U120 km and V16 km meshes are conducted from 25 December 2014 to 31 December 2015.The numbers of mesh cells in U120 km and V16 km are ∼0.04 million and ∼0.17 million, respectively, and it takes ∼19.7 CPU hr and ∼262.3CPU hr to finish a 1-day simulation for U120 km and V16 km, respectively.The difference in computational time between U120 km and V16 km comes from not only their difference in mesh size but also in time step (90 s for V16 km and 600 s for U120 km) for the non-hydrostatic dynamic core.It may be noted that these simulations are much computationally slower than those with climate models because the model is implemented with a sectional aerosol module that is more precise but expensive (Zhao et al., 2013a).Please also note that although the V16 km experiment is relatively computationally expensive, it is much more efficient than conducting global simulation at uniform 16 km horizontal resolution (e.g., U16 km with ∼2.5 million mesh cells).To demonstrate the model capability of regional mesh refinement to a few kilometers, the V4 km experiment is conducted but only for 1 month of July 2015, when convection process is most active during the South Asian monsoon season, because of its large computational and storage demand with a mesh cell number of ∼0.78 million and a much shorter dynamic time step (20 s for V4 km).The U60 km experiment is conducted for the same time period as V4 km.The V4 km and U60 km experiments take ∼3,250.5 and ∼72.7 CPU hr to finish a 1-day simulation, respectively.Here, it can be noted that time step is an important factor affecting the simulation speed.Although V4 km has variable resolution from 4 to 60 km, the time step for V4 km needs to be 20 s for all mesh cells.However, please note that, in general, the number of mesh cells outside the refined region is much less than that within the refined region, for example, less than 10% in V4 km of this study.Therefore, the computational cost for mesh cells outside the refined region is small.
The meteorological initial conditions are derived from the European Center for Medium-Range Weather Forecast (ERA-Interim data set) reanalysis data set at approximately 80 km and 6-hr temporal intervals (Dee et al., 2011).The initial meteorological conditions and topography are interpolated from the same data set to each mesh.Only the results for 2015 (referred to as the simulation period hereafter) are analyzed to minimize the impacts of the chemical initial condition.To guarantee the reliability of simulated meteorological fields, the simulation is reinitialized every 5 days with the reanalysis data.In addition, to avoid instability from the initial condition, a 1-day spin-up is considered for each 5-day simulation.Therefore, the 1-year simulation consists of a few 4-day simulations (e.g., 2-5 January of the simulation for 1-5 January and 6-9 January of the simulation for 5-9 January).The chemical initial condition (i.e., dust initial condition) of each 5-day simulation is taken from the corresponding time in the last block of simulation to ensure that the concentrations are continuous.The sea surface temperature is prescribed with the reanalysis data set and updated to the simulations every 6 hr during the simulation period.The Mellor-Yamada-Nakanishi-Niino PBL scheme (Nakanishi & Niino, 2006, 2009), NOAH land-surface scheme (F.Chen & Dudhia, 2001), Thompson microphysics scheme (Thompson et al., 2008), Grell-Freitas cumulus convection scheme (Grell & Freitas, 2014), and RRTMG longwave and shortwave radiation schemes (Iacono Experiment Center of refinement Diameter of refinement (transition)

Data Sets
To evaluate the modeling results of dust, several data sets are used in this study.The climatological means of the ground-based observations of the near-surface dust mass concentration, the deposition flux, and the AERONET AOD are used.In brief, the compiled data set of annual mean surface dust mass concentration is taken from Mahowald et al. (2009), and the detailed information of each station can be found in supplemental Table 2 of Mahowald et al. (2009).The dust deposition flux data set was compiled based on merging and revising preexisting data sets for the modern climate.More details can be found in Albani et al. (2014).The retrieved total AOD is from the AERONET network (Holben et al., 1998).In this study, the monthly mean AOD from the AERONET Version 3 Direct Sun Algorithm, Level 2.0 data set is used.A subset of stations dominated by dust is selected.
The selected stations need to meet the following conditions: (a) to reduce the impact of oceanic aerosols, only the sites that are located on land are used; (b) to reflect the annual cycle of dust AOD, the sites must contain data for the entire 12 months; and (c) the monthly average Angström Exponent (AE, 500-870 nm) at each site should be less than 0.8 because the lower the AE is, the larger the dust fraction (Dubovik et al., 2002).Considering these criteria, the number of AERONET sites becomes very limited in a single year (such as the simulation period in this study).Therefore, the multiyear (2010-2019) annual average AOD at 500 nm at each station is used with a total station number of 16.

Evaluation of Modeling Framework
Please note that the evaluation here is to convince us that the performance of this new modeling framework is comparable to other dust modeling studies on a global scale (e.g., Albani et al., 2014;Y. Feng et al., 2022;Huneeus et al., 2011;Zhao et al., 2013a).The evaluation lays the foundation for further investigating the impacts of mesh refinement on dust modeling.The evaluation is conducted through the comparison between the observations and the U120 km simulation at global scale.For regional validation, the available observations are limited in terms of both space and time.For example, only two AOD observational stations and few surface concentration observations are near the dust source region of East Asia.Therefore, it is difficult to evaluate the performance inside the refined domain, which deserves further investigation in future for specific events with more observations.The impacts of mesh refinement are assessed in Section 3.2.

Surface Mass Concentration
Figure 4 shows the dust surface mass concentrations from the observations and the U120 km simulations.Clearly, according to the observations, the surface dust concentration is higher in several major dust sources and their nearby areas (such as northern and southern Africa, the Middle East, the southwestern United States, etc.) (Figure 4a).Additionally, the surface dust mass concentration over the Northern Hemisphere is larger than that over the Southern Hemisphere.However, as dust is mainly transported from source regions, the surface dust concentration is low in remote regions (e.g., Antarctica).Furthermore, the simulated dust concentrations are evaluated with the observations at all 47 sites.The simulated mean dust concentration with a value of 21.52 μg m −3 is larger than the observation with a value of 14.31 μg m −3 (Figure 4b).Overall, the simulation can reproduce the characteristics of the spatial distribution of the observed dust surface concentration.
Figure 4c shows that the difference between the observation and simulation at almost all the observation sites is within one order of magnitude over the nine regions, as shown in Figure 4a.In general, the model performance is comparable to other dust modeling studies on a global scale (e.g., Albani et al., 2014;Huneeus et al., 2011).
For the regions of N. Africa and N. Atlantic (including six sites), the mean values from the observation and simulation are 37.20 and 67.33 μg m −3 , respectively.For the Europe region (eight sites), the mean values are 10.73 and 12.01 μg m −3 , respectively.For the E. Asia and N. Pacific regions (17 sites), the mean values are 10.61 and 2.59 μg m −3 , respectively.For the regions of S. Indian, Australia, and S. Pacific (8 sites), the mean values are 0.31 and 0.16 μg m −3 , respectively.In general, the comparison shows that the modeling results over the dust source regions are higher than the observations, while over the remote regions, the modeling results are lower than the observations.This may indicate that the removal processes in the model have some positive biases, which can lead to fewer dust particles being transported to remote regions.More in situ observations of both concentration and removal flux are needed to further evaluate and optimize the model.

AOD
Since the observed total AOD includes the contribution from all aerosol compositions, only the observed AOD at the sites near the dust source region or dominated by dust AOD can be taken to evaluate the modeling results of dust.
Figure 5 shows the comparisons of AOD at 500 nm from the observation and the simulation.The overall spatial distribution of AOD dominated by dust is consistent with that of surface dust mass concentration.The mean AOD at all 16 observation sites is 0.16, while the corresponding simulated value is 0.13.Over the nine subregions defined above in Figure 4a, there are two subregions with more than five observation sites: the region of N. Africa and N. Atlantic (five sites) and the region of C. Asia and Arabian Sea (five sites).The mean AOD from the observation and simulation in the region of N. Africa and N. Atlantic are 0.24 and 0.26, respectively; the mean AOD are 0.19 and 0.13 in the region of C. Asia and the Arabian Sea, respectively.As shown in Figure 5, the difference between the observation and simulation at most of the observation sites is within the 1:2 and 2:1 lines.Although the AOD at the selected sites should be dominated by dust, the observed AOD still includes the contributions of other aerosol compositions, such as biomass burning aerosols over North Africa, and therefore is expected to be larger than the simulated dust AOD (0.16 vs. 0.13).The averaged dust AOD from the simulation on a global scale is 0.03, which is consistent with the value of 0.03 ± 0.005 based on retrievals from multiple satellite platforms (Ridley et al., 2016).The purpose of evaluation against different observations in this study is to demonstrate that the dust modeling framework can simulate a comparable result to other global dust models (e.g., Albani et al., 2014;Y. Feng et al., 2022;Huneeus et al., 2011).More in situ observations for specific dust events may be needed to further examine and optimize the modeling schemes.

Deposition
The comparisons of total dust deposition fluxes from the observation and the simulation are shown in Figure 6.The spatial distributions of total deposition fluxes from the observation and the simulation are generally similar to that of surface dust mass concentration (Figure 4).The observation shows a large range of deposition fluxes in a few magnitudes, which is well captured by the simulation.The mean value of total dust deposition fluxes at all 110 observation sites is 8.17 g m −2 year −1 , while the corresponding simulated value is 11.20 g m −2 year −1 .Over the nine subregions, the simulated results are smaller than the observed values in remote regions such as the region of S. Atlantic, S. America, S. Africa (observation vs. simulation, and hereafter; 8.32 vs. 3.88 g m −2 year −1 ), the region of Europe (6.46 vs. 4.75 g m −2 year −1 ), the region of E. Asia, N. Pacific (3.66 vs. 2.13 g m −2 year −1 ), and the region of S. Indian, Australia, S. Pacific (1.54 vs. 0.71 g m −2 year −1 ).The simulated results are higher than or similar to the observed values near the main dust source regions.In the region of N. Africa, N. Atlantic (29 sites), the mean values of observation and simulation are 11.42 and 25.63 g m −2 year −1 , respectively.In the region of C. Asia, Arabian Sea (14 sites), the mean values of observation and simulation are 25.84 and 25.16 g m −2 year −1 , respectively.The calculation of deposition fluxes in the model is strongly affected by the dust mass.Therefore, the biases of deposition fluxes are spatially consistent with the biases of surface concentration, both of which show overestimation over the dust source regions and underestimation over the remote regions.

Dust Radiative Forcing
The dust radiative forcing here is defined as the perturbation of radiative fluxes at the top of atmosphere (TOA) and the bottom of atmosphere (BOT) and radiative heating/cooling in the atmosphere (ATM) if atmospheric dust is removed without the change of meteorological fields.The dust radiative forcing is estimated with the method introduced by Zhao et al. (2013a) in the model (details in Section 2.2).The positive radiative forcing denotes net downward radiation fluxes for TOA/BOT and warming for ATM. Figure 7 shows the spatial distributions of dust radiative forcing from the U120 km experiment at TOA and BOT and in ATM.The detailed global mean values of SW, LW, and net (NET) radiative forcing are also summarized in Table 3.The U120 km experiment shows negative dust radiative forcing at TOA (−0.27 W m −2 ).The model predicts positive radiative forcing of dust at TOA over major source regions such as the Sahara, Arabian Peninsula, and East Asia.The strongest positive forcing can exceed 10 W m −2 due to the radiative absorption of dust when located over highly reflective surfaces.In comparison, negative forcings are found over oceans or lands with relatively dark surfaces.The U120 km experiment simulates negative dust radiative forcing at BOT (−1.46 W m −2 ) and positive dust warming in ATM (1.18 W m −2 ) globally with a similar spatial pattern as that of dust mass loading.
Figure 8 shows the comparison of simulated dust radiative forcing at TOA from U120 km with the estimates in previous studies (Albani et al., 2014;Y. Feng et al., 2022;Kok et al., 2017;Scanza et al., 2015).Among previous studies, Kok et al. ( 2017) estimated the dust radiative forcing at TOA to be −0.2W m −2 with an uncertain range  (Kok et al., 2017).The bias mainly comes from weaker positive LW forcing, although the simulated negative SW forcing is also weaker.Kok et al. (2017) suggested that including coarser dust particles (10∼20 μm) might lead to larger positive LW forcing (∼0.1 W m −2 ), which could reduce the modeling bias against the observation-constrained estimate.According to a previous study (Adebiyi & Kok, 2020), coarse dust particles (with diameters of 5-20 μm) may account for 50%-69% of the total atmospheric dust mass, which is higher than the ∼26% estimated in this study, although the simulated value is within the range (6%-31%) reported by other global model simulations (Adebiyi & Kok, 2020).This may indicate that the current model deposits coarse dust out of the atmosphere too quickly and (or) underestimates the fraction of coarse dust emission, which deserves further investigation in the future.

Impacts of Regional Refinement on Dust Emissions
One of the purposes of establishing dust simulations based on the global variable-resolution modeling framework is to conduct experiments with regional refinement at a higher resolution to better resolve regional processes.As discussed in Section 2 (Figure 3), two global variable-resolution simulations (V16 km-NA and V16 km-EA) are conducted in this study and are analyzed against the uniform-resolution simulation (U120 km).First, the global and regional emissions from the three experiments are summarized in Table 4.The results show that the two V16 km experiments simulate slightly smaller global emissions (∼4% for V16 km-NA and ∼1% for V16 km-EA) than the U120 km experiment.Regionally, the difference in dust emissions is mainly confined over the refined regions of North Africa and East Asia when comparing the V16 km experiments and the U120 km experiment (Figure 9).Over the refined regions, dust emissions are generally weakened, which is consistent with some of previous studies.For example, M. Liu and Westphal (2001) found that the 80 km resolution produced nearly 20% higher dust emissions compared to the 20 km resolution using the Navy's operational mesoscale meteorological model.V16 km-NA and V16 km-EA simulate 5% and 7% less dust emissions over North Africa and East Asia, respectively, than U120 km.
Dust emission flux is generally controlled by soil wetness, source function, the threshold wind velocity, and the near-surface wind speed (see details in Section 2.2).The source functions in these three experiments are from the same data set, as shown in Figure 1, and the difference over the refined regions due to interpolation among the experiments is very small (Figure S1 in Supporting Information S1).The relative soil wetness is also nearly the same in the U120 km and V16 km experiments (Figure S2 in Supporting Information S1), which may be due to that the experiments are reinitialized every 5-day and the relative soil wetness does not change much during the simulation from the same initial condition.Therefore, the difference in dust emission flux among the experiments is mainly from their difference in near-surface winds and the frequency of dust emission events (near-surface wind velocity bigger than the threshold wind velocity).Over North Africa, the regional average near-surface wind speed from V16 km-NA is 5.0 m/s smaller than 5.2 m/s in U120 km and is closer to 4.7 m/s in the reanalysis.Over East Asia, the regional average near-surface wind speed from V16 km-EA is 4.6 m/s smaller than 5.3 m/s in U120 km and is also closer to 4 m/s in the reanalysis.The difference of threshold wind velocity for dust particles with the diameter range of 2-3.6 μm in Equation 1 between U120 km (2.92 m/s over East Asia, 2.32 m/s over North Africa) and V16 km (2.95 m/s over East Asia from V16 km-EA, 2.34 m/s over North Africa from V16 km-NA) is relatively small (Figure S3 in Supporting Information S1).The difference for particles with other sizes are also small.Based on near-surface wind and threshold wind velocity, the frequency of emission events for dust particles with the diameter range of 2-3.6 μm in Equation 1 is shown in Figure 10.The emission event is determined every hour when near-surface wind velocity is greater than threshold wind velocity.The frequency is calculated as the ratio between the hours with emission event and total hours in a year.In general, dust emission event occurs less often in V16 km experiments than in U120 km experiment inside the refined regions, particularly for V16 km-EA.The frequencies for emitted dust particles of other sizes are similar (not shown here).Both less dust emission events and weaker near-surface wind velocity contribute to the lower dust emissions in V16 km experiments inside refined regions.
The reduction in wind speed over the refined regions can be explained by more topographic complexity resolved by the refined mesh in V16 than in U120 km.Figures 11a and 11d shows the spatial distributions of topography over East Asia and North Africa from U120 km.Figures 11c and 11f shows the difference in terrain variances over the refined regions between V16 km-EA/NA and U120 km, which can represent the complexity of topography.Compared to U120 m, V16 km with higher horizontal resolution resolves more topographical complexity over the refined regions.
Previous studies have suggested that the complex topography would result in an overall weakening effect on the near-surface wind speed due to the orographic drag process and more heterogeneous spatial distribution due to resolved mountains and valleys (e.g., Li et al., 2022;Lin et al., 2018;Y. Wang et al., 2020;Zhang et al., 2020;Zhou et al., 2018).Figures 11b and 11e shows the difference in near-surface wind speed between the V16 km-EA/NA and  In addition to the difference in the refined regions between V16 km and U120 km, it is also noted that the emissions over the areas outside the refined regions can also be different (Figure 9).For example, the difference in dust emissions between V16 km-NA and U120 km is also evident over the East Asian dust source region, and the difference between V16 km-EA and U120 km is also evident over the North African dust source region.In fact, the mesh resolutions are almost identical between V16 km and U120 km over the areas outside the refined regions (Figure 3).The difference between them over the nonrefined regions is unlikely due to the difference in mesh resolution.The other two factors that may affect the simulation results from V16 km and U120 km are the dynamic time step and the coefficient for horizontal diffusion, both of which are functions of the finest resolution of the global mesh.For example, the dynamic time step is set to 600 and 90 s for U120 km and V16 km, respectively, and the horizontal diffusion coefficient is a function of the length scale of 120 and 16 km for U120 km and V16 km, respectively.Therefore, even in nonrefined areas, the simulations at similar resolutions can be affected by the two factors and produce different results.To demonstrate the dependence of the simulation results on the dynamic time step and horizontal diffusion, the difference between V16 km-NA and V16 km-EA is also shown (Figure 9b).In these two global variable-resolution experiments, their dynamic time-step and the coefficient for horizontal diffusion are the same.Over their refined regions, the difference can reflect the impact of higher resolution well and is consistent with the results over the refined regions in V16 km-EA and V16 km-NA, as shown in Figures 9c and 9d, respectively.The difference over North Africa is almost identical to that in Figure 9d, and the difference over East Asia is almost identical to the opposite of that in Figure 9c.This indicates that their difference over the areas outside the refined regions is very small because the impact of refinement over East Asia (North Africa) on North African (East Asian) dust emissions is negligible.Although the simulation results outside the refined region can be affected by dynamic time step and horizontal diffusion, please note that the cells with the difference exceeding 95% statistical significance are mainly within the refined regions (Figure 9).

Impacts of Regional Refinement on Atmospheric Dust
Figure 12 shows the spatial distributions of dust mass loading from U120 km and the difference in mass loading between these simulations.On the global average, V16 km-EA and V16 km-NA simulate similar global mean dust mass loadings of 0.041 g m −2 , ∼5% lower than 0.043 g m −2 from U120 km.Similar to the spatial distributions of impacts on dust emissions, the impacts on dust mass loading between V16 km and U120 km are mainly over the source regions but also exist over the areas outside the refined regions.As discussed above, the other two factors besides the mesh resolution, dynamic time-step and the coefficient for horizontal diffusion can also affect the simulation results.Again, the cells with the difference exceeding 95% statistical significance are mainly within the refined regions (Figure 12). Figure 12b shows that the difference between V16 km-NA and V16 km-EA is negligible in the areas outside the refined regions.Over North Africa, V16 km-NA simulates slightly (∼2%) lower dust mass loading than U120 km.Over East Asia, V16 km-EA simulates ∼10% higher dust mass loading than U120 km.In V16 km-EA, not only over the East Asian dust source region, dust mass loading is evidently higher than U120 km over the Arabian Sea, which is mainly due to more dust emissions over the Thar Desert (Figure 9).As shown in Figure 11b, the near-surface wind is higher in V16 km than in U120 km over the Thar Desert.Figure 11a shows that the Thar Desert lies between the mountains, and higher-resolution in V16 km resolves deeper valley and higher mountains, which can lead to stronger valley wind along the Thar Desert (Zhang et al., 2020).
Unlike emissions, which are mainly controlled by near-surface wind, dust mass loading can be affected by other processes in addition to emissions, such as transport and deposition.Although dust emissions over the refined regions are generally lower in V16 km than in U120 km, V16 km also simulates a lower dry deposition rate (Figure S4 in Supporting Information S1), which is also partly related to the weaker near-surface wind.Therefore, dust mass loading over North Africa is close to each other between V16 km-NA and U120 km, while over East Asia, V16 km-EA simulates higher dust mass loading due to its weaker dust dry deposition rate.Since the experiments in this study are constrained toward the reanalysis data through re-initialization every 5 days (detailed in Section 2.3) to maintain reliable meteorological fields, the impacts of mesh refinement on large-scale circulation outside the refined regions are relatively small.The impacts on wind fields and precipitation are also confined within the refined regions (Figures S5 and S6 in Supporting Information S1), which may affect regional transport and wet deposition and hence the spatial variability of atmospheric dust within refined regions, but the impacts on the regional average are small, as discussed above.In particular, the difference in precipitation between U120 km and V16 km is mainly over areas with very small amounts of dust.The global and regional budgets of dust from the three experiments are summarized in Table 4.
The simulated dust surface concentration, AOD, and deposition fluxes from the V16 km experiments are also compared with the observations.The root mean square errors (RMSE) and the correlations (R) of V16 km against the stational observations are also listed in Table 5.Both V16 km-EA and V16 km-NA significantly improve the modeling results of dust surface concentration and show slight impacts on the comparison of AOD and deposition.As discussed in Section 3.1.1,the U120 km experiment overestimates the dust surface concentration over dust source regions.The improvement in dust surface concentration mainly comes from the reduction in dust emission fluxes over refined regions (E.Asia and N. Africa, which are the main dust source regions), as discussed above.More in situ observations of dust mass concentrations and dry deposition fluxes near dust source regions will be needed to further optimize and improve dust modeling results on a regional scale with mesh refinement at higher resolution in the future.Please note that the mesh refinement is not only for improving dust modeling performance but also for better simulating the climatic effects of dust at higher resolution.At higher resolution, the complex topography can be better resolved, which is important for modeling near-surface wind speed and the transport of air particles and moisture on a regional scale (e.g., Li et al., 2022;Zhang et al., 2020) (see Section 3.2.3).In addition, the interaction between dust and cloud-precipitation-radiation and its impacts on regional climate can be better simulated at higher horizontal resolutions (e.g., Lau et al., 2020;Sakaguchi et al., 2015;Xu et al., 2021;Zhao et al., 2019).Figure 13 shows the spatial distributions of dust radiative forcing at the TOA from U120 km and the difference between the experiments.The difference in radiative forcing between V16 km and U120 km occurs mainly over the refined regions.Again, the difference over the areas  1.The emission event is determined every hour when near-surface wind velocity is greater than threshold wind velocity.The frequency is calculated as the ratio between the hours with emission event and total hours in a year.
outside the refined regions is mainly due to the different dynamic time steps and the coefficient for horizontal diffusion between V16 km and U120 km, which follows the spatial distributions of the difference in dust mass loading.The smaller dust mass loading over the refined regions leads to smaller (less negative) dust radiative forcing at TOA.The global and regional averaged dust radiative forcing in the ATM and at the BOT is also summarized in Table 3.Both V16 km-EA and V16 km-NA simulate ∼5% to ∼14% smaller dust direct radiative forcing at TOA, ATM, and BOT than U120 km, which is consistent with their smaller dust mass loading.

Impacts of Regional Refinement at Convection-Permitting Resolution
To demonstrate the capability of mesh refinement to the resolution of the convection-permitting scale (i.e., 4 km in this study, see details in Section 2.3), Figures 14 and 15 show the spatial distributions of dust emission fluxes and mass loading from U60 km and the difference between V4 km and U60 km.As mentioned above, due to the high computational cost, the experiments are only conducted for 1 month (July) of 2015.Similar to the difference between V16 km and U120 km, the impacts on dust emission and mass loading are mainly over the source regions but also exist over the areas outside the refined regions, which is mainly due to the large difference in their dynamic time-step and horizontal diffusion coefficient between the experiments.Overall, V4 km produces almost identical dust emissions as U60 km on a global average.Over the East Asian deserts and the Thar desert (marked in Figure 15), V4 km produces ∼1.4% less and ∼5.5% higher dust emissions, respectively, which corre- sponds well to the changes in near-surface wind over the areas in V4 km against U60 km (Figure S7 in Supporting Information S1).As discussed above, this is due to that higher-resolution in V4 km resolves deeper valley and higher mountains, which can lead to stronger valley wind along the Thar Desert (Zhang et al., 2020).
As shown in Figure 15c, within the refined and transition region, V4 km produces higher dust mass loading over the areas around the Thar Desert (∼20%) and its downwind area than U60 km, which is mainly affected by emission and dry deposition.The emission rates in V4 km (median value: 9.9 μg/m 2 /s) is higher than that in U60 km (median value: 8.6 μg/m 2 /s).The dry deposition rates are similar in V4 km (median value: 0.33 cm/s) and U60 km (median value: 0.34 cm/s) (Figure S8 in Supporting Information S1) in the Thar Desert.In addition to these areas with higher dust mass loading from V4 km, there is also a large area to the south of TP with lower dust mass loading in V4 km compared to U60 km, which cannot be explained by their difference in emission and dry deposition.The lower dust mass loading in V4 km is mainly within the refined regions and corresponds well to high precipitation in this area (Figure 16). Figure 16 shows the spatial distributions of total/convective parameterized/grid-scale resolved precipitation in V4 km and U60 km and their differences, in which the convec-  tive parameterized precipitation is calculated from the convective parameterization and the grid-scale resolved precipitation is calculated during the cloud-microphysics processes.The total precipitation is the sum of two.The total precipitation over the Thar Desert and its downwind area, as mentioned above, is relatively small.Precipitation is mainly over the southern part of the refined region, and total precipitation shifts southward in V4 km compared to U60 km, which results in more precipitation over North India and less precipitation over the areas to the south of the TP in V4 km than in U60 km.Compared with the observations, V4 km reduces the biases in simulated precipitation in U60 km over North India (Figure S9 in Supporting Information S1).Combined with the spatial distribution of dust mass loading over this region (Figure 16a), this difference in total precipitation contributes partly to the lower dust mass loading over North India in V4 km.
In addition, V4 km simulates much larger grid-scale resolved precipitation within the refined regions at convection-permitting resolution (i.e., 4 km in this study), while U60 km simulates much larger convective parameterized precipitation within the refined regions.More precipitation over North India from V4 km is mainly from its grid-scale resolved precipitation.As with most atmospheric models, the wet scavenging scheme is mainly based on grid-scale resolved precipitation in this model (details in Section 2.2.5).Although parameterized convective transport and scavenging is also treated in this model following Zhao et al. (2013a), it is highly uncertain, as are many other models.The difference in wet scavenging processes with grid-scale resolved and convective parameterized precipitation may also contribute partly to the difference in dust mass loading between U60 km and V4 km.Considering the high uncertainty of convective parameterization, simulation at tens of kilometers with convective parameterization may underestimate dust wet scavenging compared to that at convection-permitting resolution, which deserves further investigation with more observations of dust mass concentrations and deposition fluxes over regions with active convection systems.
Previous studies (e.g., Zhang et al., 2020) also pointed out that a modeling resolution of a few kilometers could better resolve the complex topography over the TP and thus affect aerosol transport toward the TP compared to simulations at resolutions of tens of kilometers.Over the TP (the region with terrain height larger than 4,000 m), V4 km simulates much lower dust mass loading throughout the entire TP (0.0053 g/m 2 on average within the black box) than does U60 km (0.011 g/m 2 ) (Figure S10 in Supporting Information S1).This lower dust mass loading over the TP in V4 km should result partly from its slightly lower emissions over East Asian deserts in the north and its more efficient wet scavenging in the south.Its resolution of more complex topography may also play a role.This large difference in dust mass loading around the TP may result in large differences in modeling dust-radiation interactions in the atmosphere and snowpack (Zhang et al., 2020) and dust-cloud-precipitation interactions (e.g., Lau et al., 2020).Particularly, at a resolution of a few kilometers, dust impacts on cloud microphysics can be better resolved.Investigation of dust-radiation and dust-cloud interactions and hence feedbacks to regional weather and climate systems around the TP during the South Asian monsoon season at resolutions of a few kilometers using this modeling framework deserves in future.

Summary and Discussion
In this study, a sectional dust modeling framework based on a global variable-resolution non-hydrostatic atmospheric model (MPAS) is established.The global variable-resolution modeling of dust and its radiative forcing with regional refinement to the convection-permitting scale is also introduced for the first time.In this model, atmospheric dust is simulated simultaneously with the meteorological fields, and dust-radiation interaction is also implemented in the model.Five configurations of global mesh (U120 km and V16 km with refinement over the dust source regions of East Asia and North Africa; U60 km and V4 km with refinement over the TP) are used to understand the impacts of regional refinement on the dust lifecycle at regional and global scales.The model shows the capability to reasonably produce the overall magnitudes and spatial variabilities of global dust metrics such as surface mass concentration, total deposition, AOD, and radiative forcing compared to the observations and the results from previous studies.The model simulates a global mean dust AOD of 0.03, which is consistent with the value of 0.03 ± 0.005 based on retrievals from multiple satellite instruments (Ridley et al., 2016).The simulated dust radiative forcing at TOA is −0.27W m −2 , which is within the range of −0.48∼+0.2W m −2 estimated by Kok et al. (2017) constrained by observations.It is noteworthy that here, the dust radiative forcing is defined as the direct impact of dust on radiation, which does not account for the perturbation of radiative fluxes due to dust-induced change of meteorological fields.Besides, this study only simulates dust particles without considering the complex mixture of dust and other aerosols, which may also affect the simulation of dust and its radiative impacts.The dust radiative impacts along with other aerosols in this modeling framework deserves further investigation in future.The modeling results also demonstrate some positive biases of dust surface mass concentrations near the source regions and some negative biases in the remote regions.This may reflect that the removal processes in the model have some positive biases, which leads to fewer dust particles being transported to remote regions.
Two global variable-resolution simulations with refined regions over two major deserts of North Africa (V16 km-NA) and East Asia (V16 km-EA) at a horizontal resolution of 16 km are conducted and analyzed against the uniform-resolution simulation (U120 km).V16 km-NA/EA simulates 5%/7% less dust emissions inside the refined regions due to the weakened near-surface wind speed, which is due to more topographic complexity resolved by high resolution in V16 km than in U120 km.Although dust emissions over the refined regions are generally lower in V16 km than in U120 km, V16 km also simulates a lower dry deposition rate, which is also related to weaker near-surface wind.Therefore, dust mass loading over North Africa is close to each other between V16 km-NA and U120 km, while over East Asia, V16 km-EA simulates higher dust mass loading.V16 km-EA and V16 km-NA generally have better performance in dust surface mass concentration (smaller RMSE and higher spatial correlations) compared to the observations than U120 km.The improvement of the simulated AOD in V16 km against U120 km is relatively small.
One more global variable-resolution experiment (V4 km) with refinement at 4 km resolution around the TP covering the areas of East Asian dust deserts is also conducted to demonstrate the capability for modeling dust and its impacts at the convection-permitting scale.Although similar to the difference between V16 km and U120 km, emission and dry deposition are slightly different between V4 km and U60 km, which contributes partly to the difference in their dust mass loading.The large difference in dust mass loading between V4 km and U60 km exists within the refined region to the south of the TP, which is mainly due to their difference in simulating the spatial distribution of precipitation and partitioning grid-scale resolved and convective parameterized precipitation.V4 km with convection-permitting resolution within the refined region shifts precipitation southward and simulates much larger (less) grid-scale resolved (convective parameterized) precipitation than does U60 km.The results indicate that wet scavenging with grid-scale resolved precipitation may be more efficient than that associated with convective parameterized precipitation.Therefore, simulation at tens of kilometers with convective parameterization may underestimate dust wet scavenging compared to that at convection-permitting resolution, which deserves further investigation with more observations of dust mass concentrations and deposi- tion fluxes over regions with active convection systems.As a result, V4 km also simulates much lower dust mass loading over the TP than U60 km.
This study demonstrates the capability of simulating atmospheric dust with a global variable-resolution modeling framework.Please note that the ground-based observations of dust on a global scale are only available on climatological mean, while the modeling results in this study are for one single year (2015) limited by computer resources.Therefore, the evaluation in this study is not for apple-to-apple comparison.Although the evaluation proves that the model performance is reliable and comparable to other models on a global scale (e.g., Albani et al., 2014;Y. Feng et al., 2022;Huneeus et al., 2011;Zhao et al., 2013a), the model may be applied for specific events with more in situ observations of both concentration and removal flux to further demonstrate its advantage for modeling the dust lifecycle regionally with regional refinement at high resolution.In addition to dust mass and radiative forcing, the results also show significant differences in some key meteorological fields, such as wind and precipitation, within the refinement region at high resolution.The topography is also better resolved at high resolution in the refinement region.Therefore, at resolutions of a few kilometers, dust impacts on radiation and cloud microphysics and hence the feedbacks to regional weather and climate systems may be quite different from previous studies with relatively coarse resolutions of tens of kilometers, which deserves further investigation with this modeling framework in the future.Moreover, this study mainly focuses on the impacts of horizontal resolutions.Previous studies (e.g., Menut et al., 2013;Teixeira et al., 2016) have shown that vertical grid resolution is also an important factor that might affect the dust modeling results, which deserves further exploration.
What's more, the difference between global variable-resolution and uniform-resolution experiments may also exist over the nonrefined regions even though the resolution there may be the same, which is partly related to their difference in dynamic time-step and the coefficient for horizontal diffusion, both of which are functions of the finest resolution of global mesh.This needs to be considered when comparing global variable-resolution and uniform-resolution experiments and deserves further investigation in the future.Finally, although this study demonstrates that, in terms of achieving regional high resolution, global variable-resolution simulation is much computationally cheaper than global high-resolution simulation (e.g., U16 km with ∼2.5 million mesh cells vs. V16 km with ∼0.17 million mesh cells), the simulation with regional refinement to the convection-permitting scale (a few kilometers) is still quite computationally expensive, as illustrated.To conduct numerical experiments efficiently, more effort will be made to improve the computational efficiency of this modeling framework (e.g., Gu et al., 2022).

Figure 1 .
Figure 1.Spatial distribution of the dust source function for the dust emission scheme in the model.The major dust source regions are marked on the map.

Figure 3 .
Figure 3. (a) Global quasiuniform resolution mesh of 120 km (U120 km); (b) Global variable-resolution mesh of V16-128 km with the refinement region over the East Asian dust source region (V16 km-EA); (c) Global variable-resolution mesh of V16-128 km with the refinement region over the North Africa dust source region (V16 km-NA); (d, e) Spatial distributions of mesh size for V16 km-EA and V16 km-NA.(f) Global quasiuniform resolution mesh of 60 km (U60 km); (g) Global variable-resolution mesh of V4-60 km with the refinement region over the East Asian dust source region (V4 km); (h) Spatial distribution of mesh size for V4 km.

Figure 4 .
Figure 4. Comparison of dust surface concentrations (μg m −3 ) from U120 km and observations; (a) Observations; (b) Modeling results associated with observations, the contour plot shows the modeling results and the colored circles represent the values of observed dust mass concentrations; (c) Scatter plot of simulations versus observations.The value of the Y-axis comes from the model grid nearest to the observational site.The black solid line is the 1:1 line, whereas the black dotted lines correspond to the 10:1 and 1:10 lines.Different types of markers represent different observational data sources: Crosses = University of Miami sites; Triangles = sites from Mahowald et al. (2009).Colors of markers are based on the geographical regions delineated in panel (a).

Figure 5 .Figure 6 .
Figure 5. Same as Figure4but for total aerosol optical depth (AOD) at 500 nm in dust-dominant areas.Panel (c) uses the normal coordinate system instead of the logarithmic coordinate system.The black solid line is the 1:1 line, whereas the black dotted lines correspond to the 2:1 and 1:2 lines.

Figure 7 .
Figure 7. Spatial distributions of net dust direct radiative forcing from the U120 km experiment.(a) At the top of the atmosphere (TOA); (b) in the atmosphere (ATM); (c) at the bottom of the atmosphere (BOT).

Figure 9 .
Figure 9. (a) Spatial distribution of dust emission flux from the U120 km simulation; (b) The difference in dust emission flux between V16 km-NA and V16 km-EA.The two dashed circles denote the refinement and transition zones.The region outside the thin dashed circle is at coarse resolution.(c) The difference in dust emission flux between V16 km-EA and U120 km.(d) The difference in dust emission flux between V16 km-NA and U120 km.The statistical significance of the difference in dust emission flux is estimated using the Student's test with 365 samples of daily mean value at each cell.The cell with significance value exceeding 95% is shaded by the purple symbol "+."

Figure 10 .
Figure10.Same as Figure9but for the frequency of emission events for dust particles with the diameter range of 2-3.6 μm in Equation 1.The emission event is determined every hour when near-surface wind velocity is greater than threshold wind velocity.The frequency is calculated as the ratio between the hours with emission event and total hours in a year.

Figure 11 .
Figure 11.(a, d) Terrain heights around East Asia and North Africa in the U120 km case; (b, e) Differences in near-surface wind speed over East Asia and North Africa between the V16 km-EA/NA and U120 km cases; (c, f) Differences in variance of terrain height over East Asia and North Africa between the V16 km-EA/NA and U120 km cases.The variance at each mesh cell is calculated based on all the neighboring mesh cells within 2° × 2°.

Figure 12 .R
Figure 12.Same as Figure 9 but for dust mass loading.

Figure 13 .
Figure 13.Same as Figure 9 but for dust direct radiative forcing at top of the atmosphere (TOA).A positive value refers to the downward direction.

Figure 14 .
Figure 14.(a) Spatial distribution of dust emission flux from the U60 km simulation for 1 month (July) of 2015; (b) The difference in dust emission flux between V4 km and U60 km.(c) Same as (b) but over the refined region.The Thar Desert is marked within the black box.The two dashed circles denote the refinement and transition zones.The region outside the thin dashed circle is at coarse resolution.

Figure 15 .
Figure 15.Same as Figure 14 but for dust mass loading.

Figure 16 .
Figure 16.(a-c) Spatial distributions of total/convective parameterized/grid-scale resolved precipitation from the U60 km simulation for 1 month (July) of 2015; (d-f) The difference in total/convective parameterized/grid-scale resolved precipitation between V4 km and U60 km.

Table 1
Lower/Upper Bound and Volume Fraction for 10 Size Bins of Dust in the Model Figure 2. Size distributions of volume fractions for 10 size bins as listed in Table

Table 2
Summary of the Characteristics of Five Meshes in the Experiments in This Study et al., 2000; Mlawer et al., 1997) are used.Please note that the experiments with different meshes use the same parameterizations without being tuned for each experiment.

Table 3
Annual Mean Global and Regional Dust Direct Radiative Forcing at TOA, ATM, and BOT in Numerical Experiments

Table 4
Global and Regional Budgets of Dust in Numerical Experiments