Saltation‐Induced Dust Emission of Dust Devils in the Convective Boundary Layer—An LES Study on the Meter Scale

Dust devils are vertically oriented, columnar vortices that form within the atmospheric convective boundary layer (CBL) of dry regions. They are able to lift a sufficient amount of soil particles including dust to become visible and are considered as a potentially important dust source for the atmosphere. Mineral dust, a key component of atmospheric aerosols, influences the climate by affecting the radiation budget and cloud formation. Current estimates of the contribution of dust devils to the global, regional, and local dust release vary considerably from less than 1% to more than 50%. To address this uncertainty, we perform the highest resolved large‐eddy simulation (LES) study on dust emission in the CBL to date, using the PALM model system and the saltation‐based Air Force Weather Agency (AFWA) dust emission scheme. Our results show that under desert‐like conditions, dust devils are responsible for an average of 5% of regional dust emissions, with temporary maxima of up to 15%. This contrasts with previous measurement‐based (>35%) and LES‐based estimates (∼0.1%). Local emissions of dust devils (up to 10 mg m−2 s−1) are 1–3 orders of magnitude higher than the emission in the surroundings. This makes dust devils important for air quality and visibility. Additionally, our study reveals previously unknown large‐scale convective dust emission patterns. These patterns are tied to the CBL's cellular flow structure and are the main cause of dust release. Contrary to other studies, our findings clarify the important role of saltation‐induced dust emission.


Introduction
Dust devils are atmospheric vortices with a vertical axis of rotation that frequently occur under convective conditions when the surface is heated by insulation, causing strong superadiabatic temperature gradients near the ground.During the last 80 years, dust devils have been studied extensively with field measurements (e.g., Ives, 1947;Lorenz & Lanagan, 2014;Sinclair, 1964), laboratory experiments (e.g., Kaestner et al., 2023;Mullen & Maxworthy, 1977;Neakrase & Greeley, 2010) and, especially since the 21st century, with numerical simulations, utilizing large-eddy simulation (LES) and direct numerical simulation (e.g., Giersch & Raasch, 2021;Kanak et al., 2000;Raasch & Franke, 2011).All these studies showed a wide range of values for the characteristics of dust devils, covering several orders of magnitude (Balme & Greeley, 2006;Murphy et al., 2016;Spiga et al., 2016).For example, dust devils have spatial extents of 1 m to more than 100 m horizontally, range from a few meters to more than 1,000 m vertically, and show lifetimes from a few seconds to hours.They cause pressure drops of up to several hundred pascal and maximum horizontal wind speeds of 25 m s 1 , which is why they are able to lift a sufficient amount of soil particles including dust to become visible.The lifted particles are often transported to altitudes far away from the ground by the swirling upward motion of up to 15 m s 1 .The dust fluxes provoked by dust devils and the corresponding contribution to the total atmospheric dust amount are frequently discussed (Klose et al., 2016).Especially larger dust devils might contribute significantly because they potentially lift a large amount of dust-sized particles into the boundary layer.These particles can be further transported into the free atmosphere, where they remain for several days or weeks, affecting the Earth's climate system (Knippertz & Stuut, 2014).
Atmospheric dust, as a major contributor to the atmospheric aerosol content, interacts with the climate system.The aerosols modify the radiation budget through scattering and absorbing shortwave radiation as well as absorbing and re-emitting longwave radiation (Miller et al., 2014).Additionally, they modify micro-physical processes of cloud formation by acting as ice nuclei and, thereby, influencing the cloud's feedback on the climate (Nenes et al., 2014).Moreover, dust contains a variety of organic and inorganic substances, which might serve as nutrients for the local ecology after deposition (Barcan et al., 2023), but which can also cause environmental and health issues (Morman & Plumlee, 2014).Despite this important role of atmospheric dust, its lifting mechanisms are inadequately assessed and estimates of the total global dust emission show large uncertainties, for example, Huneeus et al. (2011) proposed a range of 0.5-4 × 10 9 t yr 1 .These uncertainties can be partly explained by the contributions of small-scale phenomena such as dust devils, which are insufficiently quantified.Previous studies on the global and regional contribution by dust devils to the total dust emission did not show consistent results.Estimations based on data from the European Center for Medium-Range Weather Forecasts (ECMWF), theoretical considerations, and observational data like in-situ measurements of dust fluxes suggest a global contribution between 3.4% (Jemmett-Smith et al., 2015) and 35% (Koch & Renno, 2005).However, both studies presented large uncertainties of roughly 15%-50% (Koch & Renno, 2005) and 1%-30% (Jemmett-Smith et al., 2015) in their estimated global contributions, even though they did not take into account the variability in the estimates of the total global dust emissions, which would further increase the contribution uncertainty.Regional estimates based on the numerical Weather Research and Forecast (WRF) model, thermodynamic theory, and measurements vary between 38% for North Africa (Pan et al., 2021) and up to 53% for Western China (Han et al., 2016).Numerical simulations with the WRF Chemistry (WRF-Chem) model coupled with a new parameterization scheme for dust devils revealed a contribution in East Asia between 17.4% and 43.4% (Tang et al., 2018).Employing LES instead of large-scale weather prediction models, Klose and Shao (2016) estimated a regional contribution for Australia in the range of 0.03%-0.19%.The primary challenge in estimating the contribution of dust devils to the overall dust release is the quantification of typical dust fluxes.For dust devils, laboratory investigations, in-situ measurements and numerical simulations have not been able to adequately quantify dust fluxes of similar magnitude so far (Klose et al., 2016).While laboratory studies usually require the artificial genesis of convective vortices in a vortex chamber (e.g., Mullen & Maxworthy, 1977), measurements of dust devils suffer from the limited area that can be reliably monitored (e.g., Lorenz, 2014), and numerical simulations are mainly constrained by limited computing power (e.g., Giersch & Raasch, 2023).In addition, field studies are restricted to arid or semi-arid regions and are further complicated by the sporadic genesis of dust devils.Nevertheless, spatially fixed and portable measurement techniques were able to quantify at least basic dust devil characteristics like wind speeds and pressure drops (Balme & Greeley, 2006;Murphy et al., 2016).However, measurements of dust fluxes by dust devils are particularly difficult.Dust fluxes are determined by the product of the mass (or particle) concentration and the vertical velocity.Therefore, two quantities instead of one must be measured simultaneously (Klose et al., 2016).There have been multiple attempts to quantify dust fluxes in the field, for example, the aircraft measurements of Gillette and Sinclair (1990), LIDAR measurements of Renno et al. (2004), or measurements derived from instrumented vehicles (e.g., Mason et al., 2014).One of the most extensive field campaign to date was conducted by Metzger et al. (2011), who estimated PM10 particle fluxes in 33 dust devils to be in the range of 4 × 10 1 to 1.1 × 10 2 mg m 2 s 1 .The studies of Koch and Renno (2005) and Jemmett-Smith et al. (2015) on the global contribution of dust devils and the studies of Han et al. (2016), Tang et al. (2018), and Pan et al. (2021) on the regional contribution assumed a dust flux per dust devil of 7 × 10 2 mg m 2 s 1 , a value that can be associated with the total suspended particle flux rather than the flux of dustsized (PM10) particles (Metzger et al., 2011).This particle type dependency of the fluxes already clarifies that "typical" emission fluxes, on which the estimated contributions are based, must be carefully determined.
Assessing the statistics of dust entrainment by dust devils via measurements is challenging and very costly because dust devils of different intensities and sizes must be measured under a variety of atmospheric conditions and soil types.Nevertheless, such statistics are crucial for calculating typical fluxes and, thus, evaluating the global and regional contribution.Following Spiga et al. (2016), the numerical simulation with LES is a very promising approach to complement measurements of dust devils because it allows access to all properties of the simulated vortices such as wind speeds, temperature, and pressure as well as to local environmental conditions.Also, Neakrase et al. (2016) considers LES to be a viable option for assessing dust fluxes in dust devils.
Previous LES studies on convective vortices in the terrestrial convective boundary layer (CBL) were mostly able to reproduce the characteristic vortex structure and flow features of dust devils, similar to field measurements.However, the simulated intensities expressed through the pressure drop in the core were often too low (Ito et al., 2013;Kanak et al., 2000;Raasch & Franke, 2011).After the simulation of the first dust devils with intensities as observed (Giersch et al., 2019), Giersch and Raasch (2023) carried out a comprehensive grid sensitivity study on dust devil characteristics in LES, utilizing multiple resolutions down to a grid spacing of Δ = 0.625 m.The authors confirmed the assumption of Ito et al. (2010) that the intensity of vortices in simulations is strongly affected by the grid spacing and found that an adequate quantitative investigation on dust devils requires a resolution of at least Δ = 1 m.Apart from the required spatial resolution, a large horizontal model domain (∼10 km 2 ) is also essential to simulate dust devils of observed intensities because the occurrence of these vortices is connected to the large-scale convection pattern, which appears as polygonal cells in the vertical wind component (e.g., Giersch et al., 2019;Kanak, 2005;Kanak et al., 2000;Schmidt & Schumann, 1989).This cellular pattern is characterized by broad downward motions in the cell center and narrow upwind areas at the cell edges, also known as cell branches.It is reminiscent of a honeycomb-like pattern, or open cellular convection during cold air outbreaks.Due to flow continuity, the near-surface flow diverges beneath the downdrafts and converges beneath the updrafts.The strongest updrafts are usually found at the vertices, where several convergence lines merge.The strongest horizontal wind speeds are usually observed in regions of high horizontal gradients of the vertical velocity, that is, where up-and downdrafts alternate over short distances.In terms of the convective cells, this occurs at the transition from the broad downwind region to the narrow upwind region.Dust devils are exclusively located along the branches and vertices of the cellular pattern.
The numerical resolution requirements established by Giersch and Raasch (2023) have far-reaching consequences for studying dust fluxes with LES.As mentioned above, the grid spacing decisively influences the core pressure drops of simulated dust devils and, consequently, the horizontal and vertical wind speeds.Higher values are simulated for higher resolutions.We will show in Appendix B that also the friction velocities and dust fluxes are significantly larger for higher resolutions.Instead, the core radii of simulated dust devils decrease with better resolution (Giersch & Raasch, 2023).Therefore, dust devils cover smaller areas of higher friction velocities and higher dust fluxes at lower grid spacings, which stresses that LES-based dust fluxes are significantly affected by the resolution.The results of Giersch and Raasch (2023) suggest that studies with grid spacings much larger than 1 m lead to a significant underestimation of the dust flux.Therefore, we will focus on simulation results with a resolution of 1 m.This will be the highest-resolved LES on dust fluxes in the CBL to date.Previous investigations on dust fluxes by dust devils used grid spacings of Δ = 20 m (Ito et al., 2010) and Δ = 10 m (Klose & Shao, 2016).This is probably too coarse for a quantitative analysis and it is not surprising that the simulated fluxes and concentrations in these studies (∼10 3 -10 0 mg m 2 s 1 and ∼10 3 -10 0 mg m 3 ) are rather at the lower end of the measured values by Metzger et al. (2011) (∼10 0 -10 2 mg m 2 s 1 and ∼10-10 2 mg m 3 ).
Beside the grid spacing, we expect the choice of the dust emission scheme to strongly influence the dust flux by dust devils in LES studies.Dust emission schemes calculate the dust emission flux based on bulk properties of the atmosphere and the underlying surface.This typically includes the surface drag given by the friction velocity u * and soil properties like the particle size distribution or erodibility (see e.g., LeGrand et al., 2019;Neakrase et al., 2016).There are three different physical mechanisms of dust emission, namely direct aerodynamic entrainment, saltation bombardment, and aggregate disintegration (Shao, 2008).Emission schemes can be based on one or more of these mechanisms.Direct aerodynamic entrainment is the direct lifting of dust particles due to a strong aerodynamic drag.The most important form of direct aerodynamic entrainment is called convective turbulent dust emission (CTDE).It describes a mechanism that generates strong, localized and intermittent surface shear stresses which cause dust emission in the absence of saltation (Klose, 2014;Li et al., 2014).Schemes based on direct entrainment usually use the empirical parameterization of Loosmore and Hunt (2000) or the physics-based parameterization of Klose et al. (2014), which accounts for the stochastic behavior of inter-particle cohesive forces and the statistical distribution of momentum fluxes.However, saltation bombardment and aggregate disintegration are the most effective dust emission mechanisms on Earth (e.g., Shao et al., 1993;Tingting et al., 2018).Both require saltation as an intermediate process before the lifting of dust-sized particles can occur (Neakrase et al., 2016).Saltation describes the streamwise, hopping-motion of coarser particles or particle aggregates.As the hopping particles hit the ground, dust-sized particles are lifted and a vertical dust flux is generated (Shao, 2008).Saltation depends on several soil and surface properties like the soil moisture, the particles' density and diameter, and the distribution of vegetation and roughness elements (Bergametti et al., 2007;Shao & Lu, 2000).It is first initiated by sand-sized particles with a diameter of 80 μm as soon as the threshold friction velocity u * t of approximately 0.2 m s 1 is exceeded (Marticorena & Bergametti, 1995;White, 1979), which is usually the case in the event of dust storms (Klose et al., 2016).However, even during such strong wind erosion events, the threshold is exceeded only occasionally (Stout & Zobeck, 1997).Thus, saltation is considered as an intermittent rather than a continuous process (Shao, 2008).Intermittent saltation is also observed in the CBL when turbulent motions of air exceed the saltation threshold (Shao, 2008).The frequency of such intermittent saltation is still unclear and its statistical behavior is not well understood until today (Liu et al., 2018).Klose et al. (2016) considered it as controversial whether or not the drag in dust devils is sufficient to initiate saltation.For example, the study of Klose and Shao (2016) stated that the saltation threshold is often not reached in dust devils.We will disprove this statement in the following.
While saltation bombardment, once initiated, is a dominant dust emission process for nearly all soil types, the contribution of aggregate disintegration to the total dust emission depends more on the specific type and its properties like the amount of aggregates in the surface layer and the binding strengths of the soil aggregates (Bergametti et al., 2007).Following Shao (2008), the importance of aggregate disintegration is probably similar to that of saltation bombardment.However, the vertical dust flux due to aggregate disintegration is often not parameterized independently of saltation bombardment due to its complexity.Instead, it is assumed that it scales with the horizontal saltation flux.It can be considered in the so-called sandblasting efficiency that links the horizontal saltation flux with the vertical dust flux (Marticorena & Bergametti, 1995).State-of-the-art saltationbased dust emission schemes are the mineral dust entrainment and deposition model (L.Zhang et al., 2001) and the Air Force Weather Agency (AFWA) dust emission scheme for the Goddard Chemistry Aerosol Radiation and Transport (GOCART) model, which is part of the WRF-Chem model (Jones et al., 2011(Jones et al., , 2012;;LeGrand et al., 2019).A more complex scheme which captures both saltation bombardment and aggregate disintegration is provided by Shao et al. (2011).However, all of these schemes have mainly been used for large-scale modeling and not for local investigations (Wang et al., 2012;Y. Zhang et al., 2018;Tian et al., 2021).
Both the observed near-surface sand skirts (Murphy et al., 2016) and the high concentration of sand-sized particles (Raack et al., 2018) in dust devils indicate that these vortices provoke saltation.As the saltation-induced vertical dust flux is assumed to be one order of magnitude larger than direct aerodynamic entrainment (Shao et al., 1993), Neakrase et al. (2016) suggest to use a saltation-based parameterization to estimate the dust entrainment in dust devils.However, all previous LES studies on dust fluxes caused by dust devils utilized emission schemes based on direct aerodynamic entrainment (Ito et al., 2010;Klose & Shao, 2016).To the best knowledge of the authors, local investigations of saltation-based dust emission in the CBL have never been performed with high-resolution LES.Therefore, we make use of the AFWA dust emission scheme and perform the first high-resolution LES of saltation bombardment in the CBL.The AFWA scheme recently showed good performance in simulating and forecasting dust storms (Yuan et al., 2019) and is rather simple (LeGrand et al., 2019).The paper's structure is as follows: Section 2 gives an overview of the methodology including the PALM model system, the implemented dust scheme, the simulated setup, the detection and tracking of dust devils, and how the calculation of the contribution of the dust devils to the total dust emission is realized.The results are introduced and discussed in Section 3 with a focus on the spatial distribution of dust emission due to saltation bombardment and the contribution by dust devils to the total dust emission.A summary and conclusion completes our study.

Methodology
In the following, the PALM model system is used for the numerical simulations (Maronga, Banzhaf, et al., 2020).By default, it does not contain a physics-based parameterization of the local particle release and transport.
Therefore, the model must be coupled with such a scheme.We will focus on dust-sized particles.The coupling enables the simulation and investigation of dust fluxes, patterns of dust emission and the contribution of dust devil-like vortices to this emission.We start this section with a brief introduction to PALM, followed by an overview of the newly implemented dust physics.The simulated setup and the detection and tracking of convective vortices are described afterward.Finally, it is shown how the contribution of dust devil-like vortices to the total dust emission is determined.Note, as in other studies (e.g., Giersch et al., 2019;Kanak, 2005;Raasch & Franke, 2011), the term dust devil, dust devil-like vortex, and (convective) vortex are used as synonyms.A differentiation between visible and non-visible vortices is not made.

The PALM Model System
All numerical simulations are carried out with the PALM model system in LES mode (revision 4732).PALM is a Fortran-based code, which has been developed for studying a variety of atmospheric and oceanic flows (Maronga, Banzhaf, et al., 2020;Maronga et al., 2015;Raasch & Schröter, 2001).By default, PALM solves the nonhydrostatic, spatially filtered, incompressible Navier-Stokes equations in Boussinesq-approximated form, assuming a constant air density ρ a .Prognostic equations for up to seven variables are solved on a staggered Cartesian Arakawa-C grid: the velocity components u, v, w, the potential temperature θ, the subgrid-scale turbulence kinetic energy e, the water vapor mixing ratio q v and the passive scalar s.Dry conditions are assumed in this study (q v = 0), and the dust mass concentration is implemented via the passive scalar s.To guarantee incompressibility of the flow, a Poisson equation for the so-called perturbation pressure p* is solved by applying a predictor-corrector method after Patrinos and Kistler (1977).As in other LES studies of dust devils (e.g., Giersch & Raasch, 2023;Kanak, 2005;Kanak et al., 2000;Raasch & Franke, 2011), we determine the dust devil pressure drop with respect to the surroundings from the total dynamic pressure perturbation π* = p* + 2/3ρ a e.
For the resolved-scale advection, PALM employs the fifth-order scheme of Wicker and Skamarock (2002) together with a third-order Runge-Kutta-time-stepping scheme (Williamson, 1980).For the subgrid-scale transport, PALM follows the gradient approach, which assumes that the transport is proportional to the local gradients of the mean resolved quantities, and utilizes a 1.5th-order closure after Deardorff (1980) in the formulation by Moeng and Wyngaard (1988) and Saiki et al. (2000).
As outlined in Section 1, the friction velocity is an important quantity for dust emission.In PALM, the friction velocity is computed at each horizontal grid point through the local application of Monin-Obukhov similarity theory (MOST, Monin & Obukhov, 1954).A constant flux layer between the surface (z = 0 m) and the first computational grid level (z = 0.5Δz) is assumed, with z being the height above ground and Δz the vertical grid spacing.The calculation of u * at every grid point in the surface layer requires the knowledge of the local resolved horizontal velocity components.Thus, the near-surface horizontal wind speeds of dust devils control the magnitude of the friction velocities.For more details about PALM, the reader is referred to Maronga, Banzhaf, et al. (2020).

Implemented Dust Physics
The dust physics in PALM shall consider five processes: dust emission from the surface, passive advection with the resolved-scale turbulent wind, subgrid-scale turbulent transport, gravitational settling, and dry deposition.The individual parameterizations and calculations, which are implemented in addition to PALM's standard treatment of a passive scalar, are presented in the following paragraphs.
The dust emission parameterization follows the AFWA dust scheme (LeGrand et al., 2019), which calculates the vertical dust emission flux F e caused by saltation bombardment (Kawamura, 1951;Marticorena & Bergametti, 1995).At first, the total vertically-integrated streamwise (horizontal) saltation flux G in kg m 1 s 1 is calculated by where H(D p ) denotes the partial vertically-integrated streamwise saltation flux of the saltation size bin p with the effective particle diameter D p and dS rel (D p ) describes a bin-specific weighting factor.The bin-specific weighting factor is calculated from the mass distribution of particles in the surface soil dM(D p ) = s ssc (p) × s frac (D p ), where s scc is the mass fraction of the soil separate class (ssc: sand, silt, clay) the bin p is assigned to and s frac (D p ) is the bin-specific mass fraction in the corresponding soil separate class.For more details, the reader is referred to LeGrand et al. (2019).We assume a homogeneous surface with mass fractions based on the soil type sand of the STATSGO-FAO database with s Sand = 0.92, s Silt = 0.05, s Clay = 0.03 (Pérez et al., 2011).For the saltation size bin configuration defining s frac (D p ), we follow the recommendation by LeGrand et al. (2019).All values are summarized in Table 1.Each bin-specific vertically-integrated streamwise flux H(D p ) is calculated according to Kawamura (1951): where g is the gravitational acceleration and C mb is an empirical constant.Here, we use C mb = 1 as suggested by Marticorena et al. (1997) and Darmenova et al. (2009) instead of the original value of 2.61 according to White (1979) and Marticorena and Bergametti (1995).The threshold friction velocity, as a function of the particle diameter D p , is calculated via the semi-empirical equation of Marticorena and Bergametti (1995): with the empirical parameters a = 1.75 × 10 6 m x , b = 0.38, c = 6 × 10 5 kg m 0.5 s 2 and x = 1.56.
In a second step, the vertical bulk emission flux of emitted dust-sized particles F e in kg m 2 s 1 is determined by the product of G and the sandblasting efficiency α in m 1 through The sandblasting efficiency (LeGrand et al., 2019;Marticorena & Bergametti, 1995) is calculated via where the factor 100 results from the conversion of cm 1 , as used in the formulation of LeGrand et al. (2019), to m 1 .The variable % Clay = 3 is the mass fraction of clay in percentage.In this study, F e represents the flux of dust particles that are assumed to be uniform and spherical with a diameter of 10 μm and a density of 2,650 kg m 3 .
Note, most of the in situ measurements evaluate the dust emission via the product of the dust mass concentration c and vertical velocity w.To enable a comparison between measurements and simulations, we also determine this dust flux, which we call the vertical dust transport F t from now on.It is calculated via Mass contribution a (dM) 0.03 0.0125 0.0125 0.0125 0.0125 0.0189 0.0377 0.0330 0.3585 0.4718 a Assuming the soil category 1 ("sand") of the STATSGO-FAO database (Pérez et al., 2011).

Journal of Geophysical Research: Atmospheres
10.1029/2023JD040058 with F g and v g being the gravitational settling flux and the gravitational settling velocity, respectively.The latter can be calculated by using Stokes' law (see also Farrell & Sherman, 2015;L. Zhang et al., 2001).Field observations usually do not use a standardized altitude to evaluate dust fluxes.Measurement heights vary from less than a meter (Metzger et al., 2011;Raack et al., 2018) to several 100 m (Gillette & Sinclair, 1990;Renno et al., 2004).
In this study, we choose a height of 10 m for the assessment of the vertical dust transport F t .This height level corresponds to the vortex detection height explained in Section 2.4.
Dry deposition is implemented for the land use category "desert" based on a scheme proposed by L. Zhang et al. (2001).It estimates the dry deposition flux in a bulk-transfer formulation via where is the dry deposition velocity (here in the formulation by Zeng et al., 2020), c 1 denotes the bulk mass concentration of dust at the first computational grid layer, and R a and R s describe the aerodynamic resistances above the canopy and the surface, respectively.For more information, the reader is referred to L. Zhang et al. (2001).
As previously mentioned, the dust mass concentration field is equal to the passive scalar s in PALM.Dust emission and deposition are combined to the net surface flux F n = F e + F d , which represents the surface scalar flux in the model.F n modifies the concentrations at the first computational grid level above the surface as an additional source or sink term, depending on its sign.Gravitational settling is implemented for all heights above the surface layer and alters the local concentration as soon as a divergence of F g occurs.

The Simulation Setup
The main simulation of this study follows the setup R5N1 of Giersch and Raasch (2023).All boundary conditions, the initialization, and the numerical schemes are the same.In the following, only the most relevant settings and the discrepancies to the original setup are explained.For more details, the reader is referred to Giersch and Raasch (2023).R5N1 features a temporally and spatially constant vertical sensible heat flux of 0.24 K m s 1 at the surface to force convection.The roughness length, which needs to be prescribed for the application of MOST at the lower boundary, is set to 0.1 m.During model initialization, vertical wind and potential temperature profiles are prescribed.The velocities are set to zero because no background wind is considered and free convection is simulated.The initial potential temperature is constant (300 K) up to a height of 1,000 m and increases with 0.02 K m 1 until the top of the domain.During the beginning of the simulation, random perturbations are imposed on the horizontal wind to accelerate the development of convection and to reach a quasi-stationary state of the CBL more quickly.In addition, PALM's nesting technique (Hellsten et al., 2021) is applied in vertical direction, that is, two domains with the same horizontal but different vertical extensions are simultaneously simulated, utilizing different resolutions.The inner domain, also called child domain, spans 4,000 × 4,000 × 240 m 3 .The outer (parent) domain has a vertical extent of 2,248 m.The spatial resolutions are 5 and 1 m for the parent and child domain, respectively.To compare the results with the findings of Klose and Shao (2016), another setup called R20N10 is used, which has a parent resolution of 20 m and a child resolution of 10 m.Both setups are summarized in Table 2.The domain extents along the Cartesian coordinates x-, y-, and z are indicated by L x , L y , and L z , respectively.At the top boundary, a Neumann (zerogradient) condition is used for the dust mass concentration (∂s/∂z = 0 kg m 4 ).Note that the surface scalar (or dust) flux is not explicitly set but dynamically calculated as described in Section 2.2.The initial concentration is set to 0.0 g m 3 .For both setups, the simulation time t s is 4 hr.Following Giersch and Raasch (2023), the first 45 min are considered as model spin-up time t su , which is why the actual analysis time is defined as t a = t s t su .If not otherwise stated, all results in Section 3 refer to simulation R5N1.R20N10 is discussed in Appendix B. At this point we want to note that the R5N1 setup demands substantial computational resources, with a single simulation requiring approximately 10 days of wall-clock time on 6,900 cores of an Atos/Bull system equipped with Intel Xeon Platinum 9242 processors.

Detection and Tracking of Vortices
In the following, details of the detection and tracking algorithm of convective vortices are introduced.The algorithm is principally designed as in Giersch and Raasch (2023) with minor changes.Here, a more generalized version is presented that is not explicitly developed for a grid sensitivity study.The algorithm can be split into two parts.The first part takes care of the detection of vortex centers during the simulation.The second part filters and combines the detected centers in a post-processing step, which finally results in dust devil-like vortices that are analyzed.
During the simulation, vortex centers are identified via criteria for the modified perturbation pressure π* and vertical vorticity ζ at or slightly above an altitude of 10 m (e.g., scalar quantities are defined at 10.5 m in R5N1).
The criteria read as follows: 1.A local minimum of π * < π * th = 3 std(π * ) must be given.Its position defines the location of the vortex core.2. A local extremum of |ζ| > ζ th = 5 std(ζ) must be reached, which is located somewhere within a square of 20 × 20 m 2 around the π*-minimum.
The thresholds are based on the standard deviation (std) approach by Nishizawa et al. (2016) and set to π * th = 3.4 Pa and ζ th = 1.08 s 1 in accordance with the Δ = 1 m simulation of Giersch and Raasch (2023).The square of 20 × 20 m 2 , which limits the spatial offset of the pressure minimum and vorticity extremum, mimics typical extents of intense dust devil-like vortices in the 1 m simulation of Giersch and Raasch (2023).It allows the maximum absolute values of pressure drop and vorticity to be slightly displaced and also ensures that they belong to the same vortex center.For each detected center, the core radius R c is calculated as the distance at which the tangentially averaged modified perturbation pressure is less than 50% of its peak value at the center for the first time (see also Giersch & Raasch, 2023;Giersch et al., 2019;Kanak, 2005;Raasch & Franke, 2011).This method agrees well with empirical and analytical models of dust devils (Lorenz, 2014), where the core radius defined as above matches the location of the highest tangential velocity.
Regarding the tracking of vortices, the first step is to filter all detected centers when the simulation is finished by the following three criteria (for an explanation see further below): A. Vortex centers with core radii R c larger than 50 m are deleted.B. Centers are neglected if a stronger center (rated by π*) is found within a radius of 20 m at the same time step in order to consider the merging of vortices and to omit counting the same vortex structure with several subcenters twice or more.
Second, the remaining vortex centers are sequentially processed to generate so-called dust devil tracks, having a certain duration.Centers are assigned to the same track if the following criteria are satisfied: 1.The maximum allowed displacement between two consecutive detections is limited.It is determined by the larger value of (a) 20 m or (b) the distance calculated by a translation speed of 10 m s 1 times the time difference Δt from the previous detection of the track.2. The area-averaged vorticity ζ av (in a square of 20 × 20 m 2 around the center) must have the same sign.3. The change in π* and ζ av must be less than 10% between two consecutive centers.4. A new vortex track is initiated if no center, satisfying Criteria 1-3, is found within 3 s of simulated time.
Note, we follow the suggestion of Klose and Shao (2016) and remove all short vortex tracks with a duration of less than 30 s to increase the comparability of the data with field measurements (short-lived dust devils are hard to detect in the field) and to eliminate strong, non-coherent turbulent fluctuations that do not correspond to fully developed vortices.The theoretical and technical foundations for the criteria above (A, B, and 1-4) are well explained in Giersch and Raasch (2023), however, with a special focus on grid sensitivity.For the more generalized algorithm here, we decoupled the algorithm from its focus on the comparability for different grid spacings.The maximum core radius of Criterion A is reduced from 100 to 50 m based on a comprehensive investigation of the data from the R5N1 simulation of Giersch and Raasch (2023).This investigation showed that more than 99% of the dust devil tracks with lifetimes exceeding 30 s have mean core radii of less than 50 m, and that the remaining dust devil tracks with mean radii of more than 50 m show only weak intensities, accumulating at less than 10 Pa (not shown).Thus, we decided to neglect centers with a radius larger than 50 m.Criterion 3 is extended by the plane-averaged vorticity to better ensure that two subsequent detections belong to the same vortex track.Criterion 4 enables gaps in the vortex tracks of up to 3 s that might occur from values of the pressure drop or vorticity, which are temporarily lower than the absolute values of the applied detection thresholds.For time intervals of more than 2 s, Criterion 1 allows a larger displacement between two consecutive centers than 20 m based on the assumed maximum translation speed of 10 m s 1 , which is a reasonable value in accordance to measurements (e.g., Murphy et al., 2016).

Contribution of Dust Devils to the Dust Emission
To estimate the contribution of dust devil-like vortices to the overall dust release, an area must be defined that delimits the vortices' dust emission from the background emission.While the core area can be considered as the visible dust column (e.g., Balme & Greeley, 2006;Luan et al., 2017), the area of the total dust emission by a dust devil-like vortex is not necessarily restricted to its core.We follow the approach of Klose and Shao (2016) and assume that the relevant area for dust emission is equal to a circle of twice the core radius.This ensures that potentially high dust emission fluxes just outside the core region are also assigned to the emission fluxes by dust devils.Hereinafter, these circular regions are termed as dust devil flux areas and denoted by σ, that is, σ(n, t) stands for the flux area of the nth dust devil center of the whole sample N dds (t) detected at time t.The union of all individual dust flux areas at t is denoted as Thus, Ω(t) accounts for all emission flux relevant areas assigned to dust devils.Areas that are covered by more than one vortex (overlapping dust devil flux areas) are counted only once.In the following, ω(t) = Ω(t)/D and ω t a describe the instantaneous fractional area covered by dust devils and its time-averaged value, respectively.The horizontal domain is denoted as D and spans 4,000 × 4,000 m 2 .
The instantaneous mass flow rate Ṁ A x (t) due to a given mass flux F x is defined as the mass of lifted dust per unit time (kg s 1 ).The subscript x = {e, t} refers either to the dust emission flux F e or the vertical dust transport at 10 m altitude F t .The instantaneous mass flow rate is calculated via spatial integration of the respective flux over a certain area A with the surface elements dA: The total amount of (emitted/transported) dust mass M A x is calculated by the temporal integration of ṀA x over the analysis period t a , that is, from the model's spin-up time t su until the end of the simulation t s : Both, Ṁ A x (t) and M A x can be considered for the nth dust devil (A = σ), for all dust devils (A = Ω), or for the whole domain (A = D).If the vertical dust transport is considered instead of the emission flux at the surface, we neglect the negative values of F t and consider only the positive values of the vertical dust transport.With this restriction, the values related to F t can directly be related to the results from field measurements, which consider the positive vertical dust flux as the product of the concentration and positive vertical velocity (e.g., Metzger et al., 2011;Renno et al., 2004).With Ṁ A x (t) and M A x , the contribution by dust devils to the overall dust emission can be calculated.We distinguish between an instantaneous contribution r x (t), which is calculated via and a time-integrated contribution R x defined as: Lastly, we apply the concept of spectral frequency analysis to investigate the spectral distribution of the friction velocity, as the main simulation parameter for saltation-induced dust emission.Bins of size 10 3 m s 1 within the interval from 0 to 3.0 m s 1 are chosen.During the analysis period, u * -values at each grid point of a considered region are assigned to the corresponding bin.In this way, instantaneous frequency distributions are generated.The time-integrated frequency spectra are finally determined by an accumulation of the instantaneous distributions over the whole simulation time.We calculate both global spectra and dust devil spectra, which show the frequencies of u * over the whole simulated domain D and over the union of all dust devil flux areas Ω, respectively.

Results and Discussion
This Section clarifies the question about the contribution of dust devils to the total dust release and transport.For this purpose, we start with a domain-wide analysis of the friction velocity because of its large influence on the simulated surface dust flux.Later on, friction velocities are investigated within the dust devil flux areas.In Section 3.2, saltation-induced dust emission is analyzed in the whole simulated domain.This emission will be found to be caused by large-scale convective patterns and several exceptionally high dust fluxes associated with dust devils.Both phenomena are separately studied in Sections 3.3 and 3.4, respectively.Finally, the dust devils' contribution to the overall dust emission and vertical transport is estimated in Section 3.5.

Friction Velocities and the Saltation Threshold
Figure 1 illustrates a snapshot of the horizontal cross-section of the friction velocity for the whole simulation domain at an arbitrary time step of the simulation.To separate areas with saltation from areas without saltation, we chose a threshold friction velocity of u * t = 0.21 m s 1 , corresponding to the sixth saltation size bin of the AFWA scheme with an effective diameter of 70 μm (see Table 1).This size bin provides the minimum threshold friction velocity above which saltation of particles is possible and can be considered as AFWA's saltation threshold.
It can be seen that intermittent saltation occurs frequently during daytime convection.This conflicts with the opinion that saltation contributes only slightly to the background dust loading and that saltation thresholds are only exceeded during strong wind events like dust storms (e.g., Klose & Shao, 2016;Klose et al., 2016).In our simulation, saltation is organized and arranged along large-scale meandering patterns.These patterns are comprehensively addressed in Section 3.3 in the context of the dust emission field.The threshold friction velocity of 0.21 m s 1 is exceeded in roughly half of the horizontal area.The temporally averaged area fraction at which saltation occurs is a t a u * > u * t ≈ 52 %.In comparison, dust devils, which are visible as small light spots in Figure 1, occupy a much smaller fraction of the total horizontal area.This temporally averaged fraction is determined as ω t a = 0.16 %, utilizing the dust devil flux areas as described in Section 2.5.Combining field observations with the thermodynamic theory about natural convective as a heat engine (Rennó & Ingersoll, 1996), Koch and Renno (2005) estimated the fractional area covered by dust devils to be ω obs = 0.003% ± 0.002%, which is even smaller than ω t a .Their fractional area is defined as the region where dust devils are strong enough to produce saltation through perturbations in surface velocity.Thus, we expect more than 99% of the area where saltation is present to be outside of dust devils.This shows that saltation-induced dust emission might not only be important for the dust release of dust devils but also for the continuous, ambient dust emission during convective conditions.Local mechanisms such as strong electric fields or the Δp-effect, which are especially prevalent in dust devils, could cause a significant decrease of the threshold friction velocity (Balme & Hagermann, 2006;Esposito et al., 2016), resulting in even higher dust emissions than simulated (see Section 4).Contrary, for non-idealized surfaces, soil crusting can increase the threshold friction velocity by a factor of 2 (Pi & Sharratt, 2019), and soil moisture can lead to a further increase (Yang et al., 2019).The soil may also not contain abundant sand particles with diameters of about 70 μm or they may be shielded by larger particles.Consequently, saltation would not be initiated at a friction velocity of 0.21 m s 1 but higher.
Considering the variable values of the saltation threshold for different soils and atmospheric conditions, Table 3 summarizes the temporally averaged area fraction at which a certain friction velocity is exceeded.The value of u * t ≈ 0.2 m s 1 corresponds to the minimum value of Equation 3, u * t ≈ 0.21 m s 1 corresponds to the minimum threshold friction velocity of the AFWA scheme, and u * t = 0.4 m s 1 as well as u * t = 0.6 m s 1 follow the suggestions of Li et al. (2014) and Ju et al. (2018), respectively, to clearly separate saltation from CTDE.
Under the assumption of a higher friction velocity threshold of u * t = 0.4 m s 1 to clearly separate intermittent saltation from CTDE (e.g., Li et al., 2014), the area occupied by saltation is a t a u * > u * t ≈ 1.8%.In this scenario, roughly 90% (following ω t a ) or 99.8% (following ω obs ) of the saltation area is found outside dust devils.Here, we have assumed that saltation is active throughout the whole dust devil flux areas, which is mostly the case.Otherwise, the fraction of the saltation area outside dust devils (following ω t a ) would be even higher.Even for an area fraction of 1.8%, the ambient saltation might cause a decisive contribution to the total dust emission because saltation bombardment is considered to produce fluxes an order of magnitude larger than direct entrainment (Shao, 2008).To investigate the friction velocity distribution in more detail, we carry out a frequency analysis of u * (as described in Section 2.5).By this analysis, a better understanding of saltation-induced dust emission in the CBL is achieved.The results are illustrated in Figure 2. The frequency spectrum follows a Gaussian curve with a mean value of approximately 0.22 m s 1 and a standard deviation of 0.079 m s 1 .The tail of the distribution extends to 2.59 m s 1 .Considering the log-scale distribution from Figure 2b, the frequencies for u * > 1 m s 1 are mainly located within the dust devil flux areas.However, there is a small portion of very high friction velocity counts which is not assigned to dust devils.This is visible through the slight offset between the global (blue) and the dust devil spectra (green).The offset is caused by the algorithm for vortex identification (see Section 2.4) and the definition of the dust devil flux areas (see Section 2.5).It is discussed in detail in Appendix A. We conclude that almost all strong saltation events are an exclusive feature of dust devils.This would also explain their pronounced visibility in the field.Outside the dust devils, u * regularly exceeds the saltation threshold of 0.21 m s 1 .In some cases, values up to 0.8 m s 1 are reached.This once again stresses the important role of saltation for dust emission in the CBL.We therefore expect that observed dusty plumes, as mentioned by Koch and Renno (2005), are related to the saltation-induced dust emission caused by daytime convection and are not solely related to CTDE.
The statistical analysis of all dust devils shows that the instantaneous peak friction velocity (found during the individual dust devil lifetimes and within the dust devil flux areas σ), has a mean value (averaged over all dust devil tracks) of 0.89 m s 1 and a maximum of 2.59 m s 1 .Both values are in very good agreement with field observations by Balme et al. (2003), who derived near-surface peak friction velocities within 10 dust devils between 0.9 and 2.4 m s 1 .The friction velocity, averaged over both the individual dust devil lifetimes and over σ, has a mean of 0.32 m s 1 with a maximum of 1.28 m s 1 .The relatively low value of 0.32 m s 1 is explained by the pressure threshold |π*| ≥ 3.5 Pa used for the detection of vortex centers (see Section 2.4).If only the most intense dust devils are considered that would be able to lift a sufficient amount of dust to become visible in nature (|π*| ≥ 30 Pa, see Lorenz, 2014), the corresponding average value is 0.60 m s 1 .Apart from dust devils, temporal and spatial averaging over the analysis period and the remaining regions (D\Ω) lead to a mean friction velocity of 0.21 m s 1 , which roughly corresponds to the mean value of the Gaussian frequency spectrum discussed above.This highlights again that dust devils cover only small areas in the simulated domain and that they do not determine the overall frequency distribution of u * (except the right tail).
Note that the simulated setup considers a homogeneous roughness length of 0.1 m in agreement with Giersch and Raasch (2023).However, roughness lengths for flat sandy surfaces are generally lower by about one to three orders of magnitude (Chapman et al., 2017;Kurgansky, 2018).Therefore, our simulations might overestimate the friction velocity, whose magnitude is controlled by the roughness length.We suggest performing further simulations to investigate the impact of lower roughness lengths on the simulation results, especially concerning the friction velocity and how often the saltation threshold is exceeded.

Saltation-Based Dust Emission in the CBL
Figure 3 illustrates a domain-wide horizontal cross-section of the saltation-based dust emission flux computed at the same time as the friction velocities of Figure 1.A comparison of both Figures reveals that the friction velocity mainly controls the dust emission flux, which is the case in almost every existing emission parameterization, regardless of the considered emission mechanism (e.g., Kawamura, 1951;Klose et al., 2014;LeGrand et al., 2019;Loosmore & Hunt, 2000;Shao et al., 2011;Zender et al., 2004).In Figure 3, a very strong dust devil is visible with a peak pressure drop of about |π*| = 256 Pa.The dust devil is highlighted with the left black rectangle.The right rectangle marks an area with very strong, large-scale dust emission that is not connected to any intense vortex.More detailed illustrations of both areas can be found in Figures 5 and 6a.It can be seen that saltation-induced dust emission is organized along cellular, large-scale patterns distributed all over the domain, similar to the patterns observed for u * .Because dust emission is directly connected to the flow field, this pattern confirms that the CBL is determined on a large scale by polygonal convection cells as described in Section 1. Averaged over the analysis period, the mean dust emission flux over all locations with F e > 0 and outside of dust devils (A F e > 0 \ Ω), is 1.06 × 10 2 mg m 2 s 1 .We will refer to 10 2 mg m 2 s 1 as a typical background emission flux.Local peaks along the large-scale emission patterns are in the order of 10 1 to 10 0 mg m 2 s 1 , corresponding to a friction velocity of roughly 0.46-0.82m s 1 Figure 3 also shows that the highest dust emission fluxes are limited to very small areas.These areas can always be assigned to intense dust devils.Instantaneous peak emission fluxes of intense dust devils reach up to 46.7 mg m 2 s 1 .Therefore, dust devils clearly distinguish themselves from their surroundings with dust fluxes that are one (10 1 compared to 10 2 ) to three (10 compared to 10 2 ) orders of magnitude larger compared to the typical ambient dust emission.
Note that the application of MOST at the lower boundary is theoretically founded only for horizontally averaged quantities, but the local application has become standard in most of today's LES codes (Maronga, Banzhaf, et al., 2020).However, it is known that the local application of MOST between the surface and the first gird level causes a systematical overestimation of the averaged wind shear near the surface.Following Maronga, Knigge, and Raasch (2020), this leads to a systematical underestimation of the surface shear stress and surface friction velocity.Consequently, the general level of dust emission in the whole domain might be too small.Unfortunately, no meaningful measurement data exist for the mean background emission in the CBL that would allow a direct comparison to our values.

Large-Scale Convective Dust Emission
The large-scale patterns of dust emission are closely connected to convective motions of air in the CBL. Figure 4 illustrates both a snapshot of the horizontal cross-section of the vertical velocity w at 100 m altitude in (a), and the dust emission field together with w in (b).It is evident that the dust emission bands are located between adjacent updrafts and downdrafts, where, due to the continuity of the flow, high horizontal velocities occur.Along these regions of high horizontal velocities, we regularly find friction velocities of up to 0.5 m s 1 , which is significantly  Figure 5 displays the dust flux for the very pronounced large-scale emission structure that is highlighted by the right rectangles in Figures 3 and 4. Due to the high resolution of 1 m, details of this band-like emission pattern are well captured, revealing a large variation of F e from 0 to more than 1 mg s 1 m 2 even on short distances.The large-scale bands generally extend over several hundred meters and are composed of many, partly parallel, line-like structures of high dust emission, which follow the regions of high near-surface horizontal velocity.These structures are reminiscent of (elongated) streaks in the surface layer that have often been reported in the literature for the velocity field (e.g., Asmuth et al., 2021;Leonardi et al., 2004;Moeng & Sullivan, 1994).Analogous patterns of large-scale dust emission caused by horizontal winds due to turbulent convection were also examined in studies based on direct aerodynamic entrainment and termed as CTDE events (e.g., Ju et al., 2018;Klose & Shao, 2012;Klose et al., 2014;Li et al., 2014).However, the parameterization of CTDE usually assumes the absence of saltation and, consequently, that the respective dust emission is solely governed by direct aerodynamic entrainment.In contrast, our study reveals for the first time that convective motions of air can cause significant intermittent saltation on a large scale.Therefore, saltation should be considered in future CTDE studies.It might be decisive for the overall daytime ambient dust emission in arid and semi-arid regions.The consideration of large-scale convective dust emissions in the calculation of total global dust emissions could also potentially reduce its existing uncertainties.Note that the AFWA dust emission scheme is based on the assumption of quasi-stationary saltation, which may not always be applicable to these large-scale patterns.As this study reports these large-scale emission patterns for the first time, further research on the corresponding dust fluxes and suitable parameterizations is required.

Dust Fluxes and Concentrations Within Dust Devils
Similar to the large-scale dust emission patterns, dust devils are closely connected to the convective motions in the CBL. Figure 4a illustrates the locations of detected dust devil centers, exceeding |π * | ≥ 10 Pa at that time.We excluded weaker detections (3.5 Pa < |π * | < 10 Pa) from the illustration to make sure that the detected vortex centers do not overlap and are visually distinguishable from each other.For illustrations including all dust devil centers at a specific time, the reader is referred to Giersch and Raasch (2023).By comparing Figures 4a and 4b, we find a high correlation between large local emission fluxes (cyan color) and relatively strong vertical vortices in terms of the pressure drop (yellow dots).In addition, vortices are exclusively found at or very close to the updraft regions of the cellular pattern, which is in agreement with previous findings (e.g., Giersch et al., 2019;Kanak, 2005;Raasch & Franke, 2011).The reason is that dust devil-like vortices require strong updrafts and sufficient wind shear for their genesis and maintenance.As stated by Willis and Deardorff (1979) and Raasch and Franke (2011) both requirements are fulfilled at the vertices and branches of the convective cells.Although largescale dust emission bands might enclose weaker dust devils, very intense ones show a clear spatial offset from these bands.The large-scale emission mainly occurs directly adjacent to the updraft regions, while dust devils are preferentially located within them.The mean value and standard deviation of the area fraction occupied by dust devils is ω t a = 0.164 ± 0.039%.As already noted in Section 3.1, this value exceeds previous estimates based on field observations and thermodynamics by at least one order of magnitude.For example, Koch and Renno (2005) estimated the fractional area covered by dust devils to be ω obs = 0.003% ± 0.002%.An extensive statistical analysis by Lorenz and Jackson (2016) showed area fractions between 3 × 10 4 and 4 × 10 6 .The discrepancy to our simulation results is explained as follows: First, it is difficult to obtain good statistics on the occurrence of dust devils during field observations because only very intense dust devils are easily visible.Instead, our simulations capture the whole range of convective vortices, which agrees with the result of significantly higher detection rates in LES compared to observations (Lorenz & Jackson, 2016).Second, the area assigned to a given dust devil is not consistently defined (Klose & Shao, 2016;Koch & Renno, 2005;Lorenz & Jackson, 2016;Lorenz et al., 2021).
We recommend to revisit the definition of the dust devil flux area in future studies.Consequences of this definition are addressed in Appendix A.
Figure 6 displays snapshots of horizontal cross-sections of the surface dust emission flux for the four strongest dust devils.Vortex B features the highest pressure drop of almost 280 Pa.Instantaneous peak fluxes in the order of 10 mg m 2 s 1 at or very close to the vortex center are typical.In addition, it can be seen that the calculated dust devil flux areas capture the highest dust emission fluxes reasonably well.Note that the most intense dust devil in terms of the absolute pressure drop does not necessarily cause the highest dust emission fluxes.Instead, the highest flux is caused by a rather concentrated vortex with a core pressure drop of roughly 150 Pa and a welldeveloped central downdraft (not shown).We therefore speculate that other factors beside the vortex's intensity, like the strength of the central downdraft and the radius, result in particularly high near-surface horizontal velocities and, consequently, dust emission fluxes.If only intense dust devils are considered that would probably be visible in nature (|π*| ≥ 30 Pa, see Lorenz, 2014), we found typical peak dust emission fluxes during the vortices' lifetimes between 7.80 × 10 1 and 46.7 mg m 2 s 1 .Thus, our peak dust emission fluxes exceed the LES results from Klose and Shao (2016) by 1-2 orders of magnitude.They determined peak fluxes in the order of 10 3 to 10 0 mg m 2 s 1 .In laboratory experiments, Neakrase and Greeley (2010) determined sediment fluxes in the range of 4 × 10 0 -10 8 mg m 2 s 1 .Our peak values fit into this interval, but are much closer to the lower end than the upper end.The fluxes from the laboratory represent the bulk ranges including all sediment types (dust and sand-sized particles).The experiments further guaranteed that sufficient surface material was available for a continuous particle lifting.In addition, Neakrase and Greeley (2010) conducted their terrestrial experiments with steady horizontal wind speeds of up to 10 m s 1 and pressure drops of up to 10 hPa (1% of the Earth's ambient pressure of 1,000 hPa), which corresponds to very intense dust devils only.All this might have caused the large upper limit of 10 8 mg m 2 s 1 .Focusing on dust-sized particles, Neakrase and Greeley (2010) determined the relationship 5.68 × 10 6 |π*| 2.24 for calculating the flux.Assuming a typical pressure drop of 100 Pa for our intense dust devils, this relationship results in a dust flux of 10 5 mg m 2 s 1 , which is still several orders of magnitude larger than our maximum values.Metzger (1999) measured sediment fluxes of approximately 6 × 10 2 -5 × 10 3 mg m 2 s 1 in the field.Because these fluxes were determined for mixed sediment, including both sand and dust-sized particles, Metzger's flux estimates can be interpreted as an upper bound for our dust fluxes.Newer measurements by Metzger et al. (2011) indicate dust fluxes of 10 0 -10 2 mg m 2 s 1 , which shows again that Neakrase's lower limit of 4 × 10 0 mg m 2 s 1 from the laboratory and, thus, our determined fluxes are much more realistic.Averaged over all detected dust devils, the lifetime-and spatially-averaged (over σ) dust emission flux is 9.35 × 10 2 mg m 2 s 1 .If we apply the same averaging to the most intense dust devils that would probably be visible in nature (|π*| ≥ 30 Pa, see Lorenz, 2014), the corresponding mean value is 5.90 × 10 1 mg m 2 s 1 .
Figure 7a displays the positive vertical dust transport of Vortex A at 10 m height for the same time as in Figure 6a.
A comparison of the dust transport and surface dust emission reveals that the area of F t > 1 mg m 2 s 1 is significantly smaller compared to the area with F e > 1 mg m 2 s 1 .The peak value of vortex A is F t = 118.3mg m 2 s 1 , which is almost seven times the corresponding peak flux of F e = 17.6 mg m 2 s 1 .Averaged over its dust flux area σ, vortex A has an instantaneous mean vertical dust transport of F σ t = 14.7 mg m 2 s 1 , which is again roughly four times larger than the mean dust emission flux of F σ e = 4.0 mg m 2 .These observations highlight that F e and F t are not directly comparable in terms of amplitude and shape and that observed dust fluxes by dust devils are significantly influenced by the considered height.
The large discrepancy between F e and F t is further confirmed by a more profound statistical analysis.Averaged over the domain and analysis period, F t has a mean value of 1.25 × 10 2 mg m 2 s 1 , which is twice the mean of F e = 5.77 × 10 3 mg m 2 s 1 .Averaged over all dust devils and their lifetimes, the spatially-averaged (over σ) value of F t is 4.69 × 10 1 mg m 2 s 1 .This is five times the corresponding value of F e = 9.35 × 10 2 mg m 2 s 1 .Applying the same averaging procedure to the peak value within σ, we derive a dust transport of 2.77 × 10 0 mg m 2 s 1 , which is about 35% more than the corresponding value for F e = 2.06 × 10 0 mg m 2 s 1 .The total maxima of F t = 1.61 × 10 2 mg m 2 s 1 and F e = 46.7 mg m 2 s 1 differ by a factor of three.
The vertical structure of the dust column, which is, based on field observations, defined as the visible column of a dust devil (e.g., Balme & Greeley, 2006;Luan et al., 2017), can be related to the dust concentration field in numerical setups.Figure 7b shows an instantaneous yz-cross-section of the dust mass concentration through the center of vortex A. The results support the findings of Morton (1966) and Hess and Spillane (1990) that the observed maximum height to maximum width ratio is of order 10 for a wide range of sizes.For vortex A, the width in terms of the diameter is about 5-10 m, which would suggest a (visible) height of 50-100 m.We also find that the dust concentration field significantly tapers from the surface to a height of a few meters (3.5 m for vortex A), where the minimum horizontal extent is reached.This height interval is called Region I in the following.The contraction of the dust column in the first meters above ground agrees well with the observation of a near-surface radial inflow of dust particles (e.g., see Balme & Greeley, 2006;Sinclair, 1966).Above this first region, the dust column is sharply confined with a small, almost constant radius (depicted by Region II).For intense dust devils, Region II usually reaches heights between 10 and 50 m (approximately 17.5 m for vortex A), that is, it includes the detection height of 10 m, where the previously mentioned transports F t were evaluated.Above Region II, the horizontal extent of the concentration field slightly increase and the dust devils begin to blur.We term this area as Region III.At a certain height, the dust column is fully blurred and the horizontal extension of the visible column would be almost constant (if visible at all).We depict this height interval by Region IV, which often extends beyond elevations of 100 m and can potentially reach the top of the boundary layer.Thus, Figure 7 only captures the lowermost part of this region.The blurring effect agrees well with observations of Renno et al. (2004), who stated that dust devils at 100 m above the surface have no clear core but rather a uniform dust content.All in all, the previously defined regions match well with those established by Sinclair (1966) and revisited by Murphy  2016): Region 1 describes the near-surface, radial inflow zone that is heavily particle-loaded and often has a v-shaped form.Region 2, at an intermediate height, is characterized by strong rotation and uplift.It includes the near-vertical column of rotating dust.In the upper-most Region 3, the structure dissipates, that is, the rotation decays and the dust devil "fades" into the ambient atmosphere.Our classification can be regarded as an extension of these three regions by a fourth one as explained above.
The tapering from Region I to II is also visible in the flux fields.As previously mentioned, the area of high dust fluxes is significant narrower for F t than for F e .Due to continuity of the mass flow, the narrowing of the area of dust fluxes causes the fluxes to increase significantly, which is why F t is larger for any statistical measure that we have calculated above.Table 4 shows the peak fluxes at different altitudes for the four strongest dust devils.Within the first 20 m, fluxes for each dust devil vary up to one order of magnitude.Particularly the lowest four m show a significant increase.For example, the vertical dust transport of Vortex A increases by a factor of two from 1 to 2 m altitude.Therefore, the height at which the dust flux is determined is critical for both numerical simulations and field observations.This conclusion is also true for the σ-averaged fluxes (not shown).We strongly recommend to use a uniform height in future studies.Note, values referring to an altitude above 10 m might not be located in the dust devil flux areas due to the commonly observed tilting of dust devils (e.g., Kaimal & Businger, 1970).
Lastly, we want to address the influence of dust devils on the dust concentration within the boundary layer.The vertical dust transport by intense dust devils is visible even at altitudes of several hundred meters.Figure 8 illustrates a xz-cross-section of the concentration field averaged along the y-direction.Through the entire vertical extent of the child domain (240 m), the most intense dust devil-like vortices cause a significant increase of the yaveraged concentration compared to the ambient value.This stresses the important role of dust devils for the dust transport into higher heights.
A statistical analysis of the dust mass concentration at the detection height of 10.5 m shows that if only intense and visible dust devils are considered (|π*| ≥ 30 Pa, see Lorenz, 2014), the mean value (over all detected dust devils) of the instantaneous peak dust mass concentration within σ and during the vortex lifetime is 8.92 × 10 1 mg m 3 .The total maximum is 16.2 mg m 3 .The mean value (over all dust devils) of the temporally and spatially averaged dust mass concentration (over the individual lifetimes and σ) is 6.02 × 10 1 mg m 3 .Due to different definitions of the dust column diameter utilized in field measurement and the dust devil flux areas in simulations, mean concentrations can hardly be compared with observational data.We limit the comparison to the observed peak values.Table 5 summarizes the peak dust mass concentrations within the first 9.5 m of the four strongest dust devils.Simulated peak dust concentrations decrease with altitude, and the highest values are always found at the first grid level.In contrast to Klose and Shao (2016), who determined a mass concentration of 10 3 -10 0 mg m 3 at an altitude of 2 m, our dust devils show values between 11.3 and 15.3 mg m 3 , respectively.Thus, our values are closer to the measurements of Metzger et al. (2011), who determined PM10 peaks between 1.3 and 162.0 mg m 3 for 21 dust devils in altitudes of 1.14-4.5 m.Renno et al. (2004) measured mass concentrations inside strong dust devils and dusty plumes of 100 mg m 3 .However, they did not consider the grain-size distribution of the lifted material.As the dust mass contributes about 10% to the total suspended particles (Metzger et al., 2011), our values are again in very good agreement.Note, because the dust emission flux F e controls c at the surface and vortex A-D do not cause the simulation's peak dust emission fluxes, we expect the simulation's concentration peaks at 2 m to be even higher.

Contribution by Dust Devils
Figures 9a and 9b illustrate the instantaneous dust mass flow rates defined in Section 2.5.Both the mass flow rates caused by the dust emission flux and by the vertical dust transport at 10 m height are shown, each for the union of all dust devil flux areas and for the total domain.The instantaneous contribution by all dust devils is also displayed in (c) and (d).A statistical summary of the different mass flow rates and the corresponding contribution of dust devils to these rates is given in Table 6.According to our knowledge, this table provides the most precise and extensive overview of saltation-induced dust mass flow rates and total emitted/transported dust masses in the CBL to date.
The dust emission flux of the total domain fluctuates around 100 g s 1 and is approximately 20 times larger than the dust emission flux caused by all dust devils (Figure 9a).Therefore, the (regional) contribution of dust devils to the total dust release amounts to an average of 5%, with instantaneous fluctuations approximately between 1 and 15% (Figure 9c).Note that it remains unclear whether the AFWA dust emission scheme is well-suited for the ambient dust emission flux along the large-scale patterns.AFWA assumes quasi-stationary saltation, which might not always be satisfied along the patterns.In addition, the assumed roughness length of 0.1 m might cause too high friction velocities for flat sandy surfaces and, hence, too high values for the ambient dust emission.Considering a possible overestimation of ambient dust emission and the potential to improve the considered region which are relevant for dust emission by dust devils (see Appendix A), the contribution could increase further.A contribution of 5% is significantly less than previous regional estimates of about 30%-50% that were based on observational data and/or large-scale modeling (Han et al., 2016;Pan et al., 2021;Tang et al., 2018), but it is significantly higher    2021) assumed that the total amount of lifted dust aerosols is caused by dust storms, dust devils and dusty plumes.Thus, they completely neglected the background dust emission during daytime convection, which we found to be a main dust emission source but is most likely invisible in the field.Second, all three studies utilized an emission flux of 7.0 × 10 2 mg m 2 s 1 , as suggested by Koch and Renno (2005).This value is one order of magnitude larger than our highest dust emission flux of 46.7 mg m 2 s 1 , probably because they did not differentiate between dust and sand-sized particles and simply considered all lifted particles as dust.Regarding the fact that roughly 90% (Metzger et al., 2011) to 99% (Raack et al., 2018) of the lifted material is larger than dust-size, Koch and Renno's value is roughly one order of magnitude too large if it is applied to describe dust fluxes.Third, Klose and Shao (2016) determined a mean dust flux based on a 10 m resolution LES.This is too coarse to generate dust devils of observed intensity (see Appendix B).Therefore, they underestimated the dust fluxes.Moreover, their dust emission scheme is not based on saltation, which is active in dust devils and the dominant dust emission process (Shao, 2008).Both explain their low estimate of contribution compared to our results.
As indicated by Table 6, the dust transport values at 10 m height and the contributions by dust devils to these transports are similar compared to Ṁe .However, during the simulation, a significant accumulation of dust in the atmosphere can be observed.At 10.5 m altitude the domain-averaged concentration increases from 1.82 × 10 2 mg m 3 at 2,700 s to 5.96 × 10 2 mg m 3 at 14,400 s.This is caused by the net dust surface flux which remains positive because total emission exceeds deposition for the simulated period.Consequently, the (positive) vertical dust transport outside dust devils is enhanced on average because it scales with the concentration c (see Equation 6).This enhancement finally causes a continuous increase of Ṁ D t , as visible in Figure 9b.Dust devils, however, are not that affected by the background dust concentration increase (see ṀΩ t in Figure 9b).Their dust content is mostly governed by the local dust emission rather than the advection of ambient dust (e.g., Gu et al., 2006;Zhao et al., 2004).This in combination with the enhanced vertical dust transport outside dust devils results in an overall negative trend of the dust devils' contribution to the total positive vertical transport with simulated time (see Figure 9d).We observe that the mean contribution averaged over 1 hr is 5.71% and 3.57% for the first and last hour of the analysis period, respectively.For that reason, conclusions regarding the contribution of dust devils to the total vertical dust transport crucially depend on the background concentration of atmospheric dust.For future studies, we suggest to either analyze both positive and negative transports at certain altitudes or to follow a quadrant analysis approach that focuses on the turbulent transports as a deviation from the mean (Y.Zhang et al., 2018).Note, the calculation of the dust devil flux areas is a critical step for the determination of the contribution of dust devils to the total dust emission and transport.Therefore, we address this issue in Appendix A in more detail.

Summary and Conclusion
In this study, saltation-induced dust emission fluxes in the dry atmospheric CBL were simulated, focusing on terrestrial dust devil-like vortices.The local and regional contribution of dust devils to the overall dust release of PM10 particles was estimated, which might have strong consequences for the consideration of dust devils in the dust budget of the climate system (Klose et al., 2016).
The numerical simulations were performed with the PALM model system, utilizing the LES approach.The model core was extended with a dust physics scheme to consider the emission, passive advection, gravitational settling, and dry deposition of dust.With the help of PALM's nesting technique, very high resolution LES were performed.For the first time, grid spacings down to the meter scale were used to simulate the saltation-induced dust emission in a simultaneously large domain of about 4 × 4 × 2 km 3 .Such grid spacings follow the resolution guidelines of Giersch and Raasch (2023) for simulating dust devils in the CBL.So far, similar studies (e.g., Ito et al., 2010;Klose & Shao, 2016) have used too coarse grid spacings of 10 m or more and have not applied a dust emission a Notation as in Figure 9.
Journal of Geophysical Research: Atmospheres 10.1029/2023JD040058 parameterization based on saltation bombardment, which is one of the key processes for the release of soil particles into the atmosphere (Shao, 2008).
The simulated friction velocity, as the main simulation parameter that controls saltation and, consequently, dust entrainment into the atmosphere, agrees well with measurements.Balme et al. (2003) derived peak near-surface friction velocities in dust devils of 0.9-2.4m s 1 .We showed that peak values of up to 2.59 m s 1 can occur.Typical maxima of the friction velocity during the vortices' lifetimes amounted to 0.89 m s 1 .However, the threshold friction velocity above which saltation can occur is much more often exceeded outside dust devils.Assuming a threshold of 0.21 m s 1 , we showed that more than 99% of the area where saltation was present was not covered by dust devils.This relatively high proportion could partly be caused by the friction velocity assumed in the simulation, which is rather at the upper limit for typical desert-like conditions.
The simulated dust emission fluxes of dust devils fit very well to the most extensive terrestrial field measurements of dust fluxes to date, which indicated values in the range of 10 0 -10 2 mg m 2 s 1 .Our fluxes were simulated in a range between 10 1 and 10 mg m 2 s 1 , while in the surroundings mean emission fluxes of 10 2 mg m 2 s 1 occurred.Thus, the local contribution of dust devils to the dust release can be significant but also varies strongly, which is in agreement with the conclusions of Klose et al. (2016).The maximum flux of 46.7 mg m 2 s 1 was caused by a dust devil with an instantaneous pressure drop of approximately 150 Pa and a maximum tangential wind velocity of 7.4 m s 1 .Averaged over all dust devils, the mean dust emission maxima during the vortices' lifetimes amounted to 2.06 mg m 2 s 1 and the temporally as well as spatially averaged (over the lifetime and the horizontal vortex sphere) dust emission fluxes showed a mean value of 9.35 × 10 2 mg m 2 s 1 .To the best knowledge of the authors, this was the first comprehensive statistical evaluation of dust emission fluxes by dust devils.Moreover, the values above indicate that previous LES studies significantly underestimated dust fluxes of dust devils, like the study from Klose and Shao (2016) (10 3 -10 0 mg m 2 s 1 ) or Ito et al. ( 2010) (10 3 -10 1 mg m 2 s 1 ).For future studies that rely on dust fluxes by dust devils, we suggest to use the value of 5.90 × 10 1 mg m 2 s 1 as the typical dust emission flux for intense dust devil-like vortices that would also be visible in the field.This value corresponds to the mean temporal and spatial average (over the vortex lifetime and horizontal sphere) of all detected dust devils with |π*| ≥ 30 Pa and a minimum duration of 30 s.
Finally, we estimated the mean contribution of dust devils to the total dust emission for desert-like regions on Earth to be approximately 5%.This is much less than previous estimates from regional studies for North Africa (38%; Pan et al., 2021), Western China (53% including dusty plumes; Han et al., 2016), or East Asia (30,4%;Tang et al., 2018), but much larger than the estimate of Klose and Shao (2016) for Australia (0.03%-0.19%).The resolution in numerical simulations and the considered dust emission phenomena in the individual studies are the main uncertainty factors that cause this variety.The rather low contribution in our case could be attributed to large-scale patterns of relatively strong dust emission, which are tightly connected to the cellular convection pattern of the CBL and dominate the overall dust release.As we did not investigate cases with prescribed background winds caused by, for example, dust storms triggered by synoptic lows, the simulated ambient dust emission might be underestimated.Similarly, the consideration of direct aerodynamic entrainment as a further dust emission process might also enhance the ambient dust emission.Both background winds and the inclusion of direct aerodynamic entrainment could reduce the relative contribution of dust devils to the overall dust release.On the other hand, our setup utilized a roughness length of 0.1 m, which is relatively high for desert-like conditions.Consequently, the resulting ambient friction velocity and, hence, the background dust emissions might be overestimated, potentially resulting in a higher relative contribution of dust devils.
In future work, all relevant dust release processes should be implemented into the simulation.Apart from the saltation-based emission, this includes aggregate disintegration and direct aerodynamic entrainment (Shao, 2008).More advanced environmental conditions should also be incorporated.Especially effects of surface properties, heterogeneities, and background winds on dust emission and transport need to be investigated because these parameters strongly influence the simulated intensity of dust devils (Giersch et al., 2019) and the ambient dust emission.Besides, the consideration of the so-called Δp-effect (Balme & Hagermann, 2006) and electrical fields (Esposito et al., 2016;Franzese et al., 2018;Kok & Renno, 2006) in dust devils could significantly modify the simulated dust fluxes.A more technical study about the dust emission area attributed to dust devils is also missing.We followed the procedure by Klose and Shao (2016) to define this area as a circle of twice the core radius.However, there is neither a theoretical foundation for this approach nor measurements or simulations that support the application of such a circle.Therefore, we highly recommend to investigate the emission area relevant for dust devils for different vortex features in more detail because it significantly determines the emission contribution by dust devils.
simulation performs best in replicating the entire range of observed peak values.Similarly, the 2.5 m simulation produces peak friction velocities within the observed range, albeit primarily toward the lower end.Coarser resolutions than 2.5 m significantly underestimate the peak values.Particularly the resolution of 10 m, as applied by Klose and Shao (2016), underestimates the maximum observed friction velocities by a factor of about 3.5 (2.4/ 0.69).We assume that the significantly lower values for the coarser resolutions are cause by the reduced peak dust devil intensities with increasing grid spacing, as discussed in Giersch and Raasch (2023).Overall, the above findings confirm the grid requirements of Giersch and Raasch (2023) to use a resolution on the meter scale if quantitative dust devil analysis shall be performed with LES.et al., 2018) to 10% (Metzger et al., 2011) of the total lifted particle mass to be in the regime of fine dustsized particles.According to Raack et al. (2018), the mass distribution within dust devils is dominated by coarser sand-sized particles.However, we plan to repeat the grain-size-resolved simulation utilizing a grid spacing of 1 m.

Figure 1 .
Figure 1.Instantaneous horizontal cross-section of the friction velocity for the whole domain.

Figure 2 .
Figure 2. Frequency spectra of the friction velocity based on approximately 10 12 counts.(a) displays the frequency N by using a linear y-axis, whereas in (b) a logarithmic y-axis is used.Frequencies of the whole domain (D) are displayed in blue, while the frequencies within the dust devil flux areas (Ω) are given in green.The vertical orange line marks the saltation threshold of the Air Force Weather Agency scheme.The red area indicates the interval (0.9 m s 1 , 2.4 m s 1 ) measured by Balme et al. (2003).

Figure 3 .
Figure 3. Instantaneous horizontal cross-section of the saltation-based dust emission flux.Areas with a vanishing dust flux are displayed in white.The color scale changes at F e = 10 1 mg m 2 s 1 (as marked by the white line).The left black rectangle contains the most intense dust devil at this time step and the right rectangle contains a very strong structure of largescale ambient dust emission.

Figure 5 .
Figure 5. Instantaneous horizontal cross-section of the surface dust emission, focusing on the strong large-scale dust emission pattern, corresponding to the right rectangle of Figure 3.The red circle depicts the dust devil flux area of a vortex with |π*| ≥ 10 Pa.

Figure 4 .
Figure 4. Instantaneous horizontal cross-sections of the vertical velocity at an altitude of 100 m in combination with detected vortex centers (yellow dots) in (a) and with the surface dust emission flux in (b).Only centers with |π*| ≥ 10 Pa are considered so that dust devils visible in (b) can be clearly assigned to detections in (a).The rectangles are the same as in Figure 3.

Figure 6 .
Figure 6.Instantaneous horizontal cross-sections of the surface dust emission flux for the four most intense dust devils at the time of their pressure minima.The red circles indicate the dust devil flux area (sphere with a radius of two times the core radius).

Figure 7 .
Figure 7. (a) Instantaneous horizontal cross-section of the (positive) vertical dust transport at 10 m height around Vortex A. The red circle depicts the dust devil flux area.The blue dotted line marks the location of the vertical cross-section (yz-plane) of the dust mass concentration field in (b).The orange dashed lines separate different height intervals (see text).

Figure 8 .
Figure 8. Vertical xz-cross-section of the dust mass concentration at 5,812 s, averaged over the whole domain length L y perpendicular to the cross-section.Positions of the strongest vortices with |π*| > 30 Pa are marked by the red arrows.

Figure 9 .
Figure 9. Instantaneous dust mass flow rates caused by (a) the dust emission flux Ṁe and by (b) the positive vertical dust transport Ṁt , each for the union of all dust devil flux areas Ω and for the total domain D. Instantaneous and time-averaged contribution of the dust devils to (c) the total emission flux r e and (d) the positive vertical dust transport r t .The areas shaded in light green represent the intervals of ± the standard deviation.

Figure B2 .
Figure B2.Friction velocity (left) and dust emission flux (right) within a square of 150 × 150 m 2 around the two strongest dust devil centers observed in P20N10 and R5N1.The black (red) circles indicate the dust devil flux area defined as a sphere with a radius of two times the core radius.

Table 1
Configuration of Saltation Size Bins and Associated Attributes for the Air Force Weather Agency Scheme

Table 2
Domain Size and Number of Grid Points for Both Simulated Setups a a In the parent domain, vertical grid stretching is applied.

Table 3
Temporally-Averaged Area Fraction a t a u * > u * t at Which a Given Threshold Friction Velocity for Saltation u * t is Exceeded

Table 4
Peak Values of F e and F t Around the Four Strongest Dust Devils at Different Altitudes aValues are given for an area of 150 × 150 m 2 at the time where the vortices reach their maximum intensity.The height of 0 m refers to F e .All other heights show the vertical transports as defined in Section 2.2.Fluxes are given in mg m 2 s 1 . a

Table 5
Peak Values of the Dust Mass Concentration c Around the Four Strongest Dust Devils at Different Altitudes a Values are given for an area of 150 × 150 m 2 at the time where the vortices reach their maximum intensity.Concentrations are given in mg m 3 .The values at 2.0 m altitude are linearly interpolated from the adjacent grid levels. a

Table 6
Statistics of the Mass Flow Rates and the Corresponding Contribution of Dust Devils a