How Do Dust Devil-Like Vortices Depend on Model Resolution? A Grid Convergence Study Using Large-Eddy Simulation

Dust devils are organized convective vortices with pressure drops of hundreds of pascals that spirally lift surface material into the air. This material modifies the radiation budget by contributing to the atmospheric aerosol concentration. Quantification of this contribution requires good knowledge of the dust devil statistics and dynamics. The latter can also help to understand vortex genesis, evolution and decay, in general. Dust devil-like vortices are numerically investigated mainly by large-eddy simulation (LES). A critical parameter in these simulations is the grid spacing, which has a great influence on the dust devil statistics. So far, it is unknown which grid size is sufficient to capture dust devils accurately. We investigate the convergence of simulated convective vertical vortices that resemble dust devils by using the LES model PALM. We use the nesting capabilities of PALM to explore grid spacings from 10 to 0.625 m. Grid spacings of 1 m or less have never been used for the analysis of dust devil-like vortices that develop in a horizontal domain of more than 10 km2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^2$$\end{document}. Our results demonstrate that a minimum resolution of 1.25 m is necessary to achieve a convergence for sample-averaged quantities like the core pressure drop. This grid spacing or smaller should be used for future quantifications of dust devil sediment fluxes. However, sample maxima of the investigated dust devil population and peak velocity values of the general flow show no convergence. If a qualitative description of the dust devil flow pattern is sufficient, we recommend a grid spacing of 2.5 m or smaller.


Introduction
Dust devils are columnar convective vortices with a vertical axis of rotation. On Earth, they typically extend to 1-10 m horizontally and 10-100 m vertically and last for several minutes B Sebastian Giersch giersch@meteo.uni-hannover.de (Balme and Greeley 2006). Dust devils mainly occur in arid regions like deserts ) and swirl up loose material from the ground, which makes them visible. The basic meteorological conditions for dust devil development are clear skies with strong insolation that favors a large super-adiabatic lapse rate near the surface and background winds below certain critical values (∼ 5 m s −1 ) (Ives 1947;Williams 1948;Horton et al. 2016;Giersch et al. 2019). These conditions cause a convective boundary layer characterized by a nearsurface cellular pattern with a size on the order of the boundary layer height (∼ 1 km). This pattern shows small cell boundaries where the horizontal flow converges and the air rises. Simultaneously, the cell centers contain broad regions of descending air (e.g., Khanna and Brasseur 1998). Dust devil-like vortices preferentially appear at the cell branches and vertices (Kanak 2006;Raasch and Franke 2011).
Typical characteristics of dust devils include high vertical vorticity as well as positive temperature and negative pressure deviations (Sinclair 1964;Mullen and Maxworthy 1977). Their mean flow features a radial inflow near the surface, with highest wind speeds just outside the dust devil core, and a spiraling upward motion with positive peak vertical velocities located near the border of the visible dust column similar to the maximum tangential velocities (Sinclair 1973;Balme and Greeley 2006). However, the short-term flow structure is often more complicated and is influenced by, for example, secondary vortices (Sinclair 1973;Zhao et al. 2004) or central downdraughts (Kaimal and Businger 1970;Sinclair 1973).
In addition to the impressive appearance and dynamical behavior described above, dust devils carry particles into the atmosphere that may effect cloud microphysics (DeMott et al. 2003) and the radiation budget (Myhre and Stordal 2001). The contribution of dust devil sediment flux to the terrestrial dust budget is still under debate (Koch and Rennó 2005;Neakrase and Greeley 2010;Jemmett-Smith et al. 2015;Klose and Shao 2016). Also, dust devils modify the vertical momentum and heat fluxes significantly (as shown by field measurements from Kaimal and Businger (1970)). All this motivates the ongoing research on dust devils.
Since the late 1990 s, direct numerical simulations (DNS, e.g., Cortese and Balachandar 1998;Giersch and Raasch 2021) and large-eddy simulations (LES) (e.g., Kanak 2005;Ohno and Takemi 2010a;Raasch and Franke 2011;Spiga et al. 2016;Giersch et al. 2019) have been increasingly used to study convective vertical vortices. One crucial parameter in numerical set-ups is the grid spacing, which often affects the simulated flow dynamics and the statistics of dust devils. However, meaningful numerical results must not depend on the model resolution. A prerequisite for an adequate simulation of dust devils is the sufficient resolution of the mean flow and its low-order moments. For daytime convective planetary boundary layers (PBLs) in which dust devils can form, Sullivan and Patton (2011) proposed a resolution of z i /Δ f > 60 so that the majority of the low-order moments (means, variances, and fluxes) become grid independent in the boundary layer interior (0.1 < z/z i < 0.9). z i and Δ f describe the boundary layer height and LES filter width, respectively. In our LES simulations, Δ f can be interpreted as the spatially uniform grid spacing Δ. For z i /Δ f > 230 (Δ ≈ 5 m), Sullivan and Patton (2011) demonstrated a coupling between large-scale thermal plumes and dust devil-like vortices. The dust devil cores tended to develop in the branches or spokes of the surface updrafts. Wurps et al. (2020) concluded that for convective simulations (z i ∼ 1000 m) a resolution of 20 m is adequate to accurately capture the profiles of the mean flow variables, the resolved turbulence kinetic energy (TKE), the resolved turbulent shear stresses, and the energy spectra. From Bopape et al. (2020), it can be inferred that a minimum grid width on the order of 25 m is sufficient to provide an accurate simulation of the convective boundary layer. However, the aforementioned studies focused on the grid-resolution requirements for low-order moments or flow structures with scales much larger than those of dust devils.
Numerical simulations by  (20 and 50 m grid spacing) and Zhao et al. (2004) (50 m grid spacing vertically and 100 m horizontally) indicated that relatively coarseresolution LES are capable of resolving basic dust devil characteristics qualitatively but  admitted that higher resolutions would affect the intensity and structures of the vortices. To simulate smaller dust devils, they suggested a grid spacing of a few meters. Resolutions of several tens of meters only capture very large dust devils, which are quite rare in nature (Kurgansky 2006). Even the largest vortices may not be accurately captured but only the large-scale thermal updraughts associated with them (Sinclair 1969). According to Kanak et al. (2000), horizontal wind speeds of the dust devil-like vortices are expected to strengthen and their circulation diameters are expected to decrease with increasing horizontal resolution. Kanak (2005) performed high-resolution LES with 2 m grid spacing that enabled the detection of smaller vortices. Despite the high resolution, vortex characteristics did not match the observations, especially the lower vortex vertical wind speed and the absolute pressure drop at the center. Raasch and Franke (2011) and Giersch et al. (2019) partly focused on the effects of grid resolution on simulated dust devils. Raasch and Franke (2011) concluded that vortices were stronger and appeared more frequently in high-resolution runs, whereas vortex diameter and height were the same. Therefore, a 2 m grid spacing should be sufficient to simulate the typical spatial structure of vortices. However, they only compared two different resolutions (1 and 2 m) and used model domains that were far too small to capture the large-scale cellular patterns to which the dust devils are tied. By comparing a 10 m to a 2 m simulation, Giersch et al. (2019) obtained a similar result as Raasch and Franke (2011). Dust devil-like vortices were more numerous and intense with finer model resolution. Giersch et al. (2019) also showed that the mean radius and lifetime of a whole dust devil sample decreased with higher resolutions and stated that a further reduction of the grid spacing (below 2 m) might result in even higher core pressure drops. Klose and Shao (2016) concluded that their number density, defined by the number of dust devils per square kilometer and per hour was smaller than in other studies due to the lower horizontal resolution they used (10 m). This was explained by an underestimation of the number of small dust devils. Interestingly, Ito et al. (2013) found no significant change in the dust devil strength (measured by vertical vorticity) between model runs with 50 and 5 m grid size. Both resolutions suggested a typical value of ∼ 10 −1 s −1 . Zhao et al. (2004) and Gu et al. (2008) simulated dust devils at ultra-fine resolution (down to 0.1 m). These studies modeled single, idealized dust devils in a cylindrical domain by applying appropriate boundary and initial conditions to force vortex development. While the results helped to understand basic vortex physics, it is unclear how applicable they are to atmospheric dust devils in more realistic environments.
For individual dust devils, some of the aforementioned studies showed that higher resolutions reveal more flow details. For example, Zhao et al. (2004) examined the flow structure of dust devil-like vortices at different grid spacings. Their vortex at 100 m resolution was characterized by local maxima in vertical vorticity, vertical velocity and temperature. However, their high-resolution vortex at 0.1 m suggested that these maxima do not necessarily occur at the same location, which is consistent with the simulations of Raasch and Franke (2011). The high-resolution vortex of Zhao et al. (2004) also revealed dynamic flow features like downdraughts at the core or several intense secondary vortices. Both phenomena were not observed in their coarse-resolution vortex. As admitted by Kanak et al. (2000), structures at the dust devil scale can not be resolved by a grid size of several tens of meters. Instead, the larger-scale circulations in which dust devils are embedded are captured. Our investigations of selected vortices will show how smaller grid spacings modify the appearance of dust devils and reveal many more flow details. However, a detailed investigation of the flow dynamics is By determining the resolution where convergence occurs, future studies can improve estimations of particle concentrations and fluxes within dust devils. Numerical simulations require emission schemes that connect the turbulent wind to the particle release and transport (e.g. Klose and Shao 2013). Therefore, a prerequisite for these schemes to work appropriately is a well-resolved flow within and around dust devils. Currently available concentration and emission data vary a lot and are subject to great uncertainties. For example, Rennó et al. (2004) found through field experiments a typical dust content and vertical dust flux of 0.1 g m −3 and 0.1-1 g m −2 s −1 , respectively. Neakrase and Greeley (2010) estimated the sediment flux of dust devils by means of laboratory investigations to be between 10 −3 and 10 3 g m −2 s −1 . Estimates from LES studies indicate dust concentrations of 10 −5 -10 −4 g m −3  and 10 −4 -10 −3 g m −3 (Klose and Shao 2016). In these studies, grid spacings of 20 and 10 m were applied, which is too high for an adequate quantification of the particle load, especially within smaller dust devils with a horizontal size of a few meters. Both  and Klose and Shao (2016) mentioned this insufficient resolution with respect to smaller vortices. Table 1 provides a summary of simulation parameters that have been used to study terrestrial dust devil-like vortices by LES while being able to capture the large-scale cellular pattern (about 1 km × 1 km horizontal domain size or more) and while using a grid spacing of 10 m or below. L i describes the spatial extents along the x-, y-, and z-axis. All grids were stretched vertically except the one of Ito et al. (2013). For all we know, a grid spacing of 2 m is the highest resolution that has ever been used so far to investigate terrestrial dust devils in the convective PBL. Higher resolutions at large domains or larger domains at high resolutions were computationally too intensive. We overcome this problem by using a nesting technique (Hellsten et al. 2021) that enables simulations of a 4 km × 4 km large horizontal model domain with grid spacings down to 0.625 m.
The above considerations indicate that there is a strong need to know which model resolution is necessary to adequately resolve terrestrial dust devils with LES and to capture their characteristics qualitatively and quantitatively correctly. Considering this, the paper structure results as follows. Section 2 gives an overview of the dust devil detection and analysis algorithm, the analyzed quantities characterizing a dust devil, the PALM model system, which is used in this study to perform the LES simulations, and the numerical set-ups. The results are introduced and discussed in Sect. 3 with a focus on the convergence behavior of general flow features, the dust devil statistics, and their three-dimensional structure. In Sect. 4, a summary and conclusions are given.

Methodology
This section introduces the numerical and analysis methods that are used in our study. First, the PALM model system with its nesting technique is presented, followed by an introduction to the numerical set-ups that are performed. We will only concentrate on those features of PALM (formerly an abbreviation for Parallelized Large-eddy Simulation Model but now an independent name) which are actually used. Second, it is clarified how a dust devil center is detected during the numerical simulation and how centers are combined to form a dust devil with a certain lifetime. To conclude this section, specific dust devil quantities are introduced, which are statistically analyzed and provide the basis for the grid convergence study. Note, (convective) vortex, dust devil, and dust devil-like vortex are used as synonyms in the following and do not distinguish between a dust-laden vortex and an invisible one. Neither dust lifting nor dust transport processes are included in our simulations and, consequently, no potential effects of the particles on the turbulent flow (e.g., Richter and Sullivan 2013) are considered.

The PALM Model System and its Nesting Technique
The numerical simulations are carried out with the PALM model system (revision 4732, Raasch and Schröter 2001;Maronga et al. 2015Maronga et al. , 2020, which is an open source model code written in Fortran designed for atmospheric and oceanic boundary-layer flows. It is publicly available at http://palm.muk.uni-hannover.de/trac/browser/?rev=4732 and designed for massively parallel computer architectures with distributed memory, utilizing the message passing interface (e.g., Gropp et al. 1999). In its default state, PALM calculates the flow solution based on the non-hydrostatic, filtered, Navier-Stokes equations in a Boussinesqapproximated form by solving the conservation equations for momentum, mass, and internal energy on a staggered Cartesian Arakawa-C grid (Arakawa and Lamb 1977). On such a grid, scalars like the potential temperature θ are defined in the grid box center, whereas the velocity components u, v, and w are shifted by half of the grid spacing in x-, y-, and z-direction, respectively. Hence, the velocities are defined in the middle of the side walls of the grid box. The filtering of the smallest eddies is done implicitly based on a spatial scale separation approach after Schumann (1975). The spatial discretization of the model domain is realized through finite differences with an equidistant horizontal and a variable vertical grid width. A third-order Runge-Kutta scheme (Williamson 1980) is used for the integration in time. Advection is discretized by a fifth-order scheme of Wicker and Skamarock (2002). The parameterization of the SGS turbulence uses a 1.5-order closure after Deardorff (1980) but in the revised formulation of Moeng and Wyngaard (1988) and Saiki et al. (2000). PALM's data output is based on the [UC] 2 data standard (Scherer et al. 2020) and is performed with the open, self-describing netCDF format.
One main assumption of the Boussinesq-approximated system of equations is the incompressibility of the flow. To reach this, a Poisson equation for the so-called perturbation pressure p * is solved by applying a predictor-corrector method (Patrinos and Kistler 1977) together with a fast Fourier transform after every Runge-Kutta sub-time step. The total dynamic pressure perturbation π * = p * + 2/3 ρe is interpreted as the dust devil pressure drop with respect to the surroundings, similar to other LES studies of dust devils (e.g., Kanak et al. 2000;Raasch and Franke 2011;Giersch et al. 2019). The second term includes the air density ρ and the SGS TKE e. It describes the isotropic part of the SGS tensor that arises from the filtering of the model equations.
In order to allow for a high resolution of small-scale turbulent processes near the surface and to simulate a sufficiently large model domain to capture the large-scale cellular flow pattern in the convective PBL, PALM's self-nesting capabilities are used (Hellsten et al. 2021). In our study, we apply the pure vertical (one-dimensional) nesting, where the high-resolution nested and coarse-resolution domains (also called child and parent domains) have identical horizontal model extents and where the children only obtain their boundary conditions from the parent at the top boundary instead of all boundary surfaces. Identical horizontal model extents are mandatory to capture the large-scale cellular pattern in all domains and enough high-resolved dust devil-like structures for statistical analysis. Regions above the dust devils can be simulated with the coarser parent resolutions to spare computational costs. The domain with the coarsest resolution is called the root domain.

Numerical Set-ups
To study the effects of model resolution on the dust devil statistics, we performed differently resolved and mostly nested simulations while maintaining the horizontal model extent. In the following, an overview about the simulation initialization is given (see also Fig. 1). To force convection, a homogeneous heating is prescribed at the surface with a temporally and spatially constant vertical sensible heat flux wθ 0 of 0.24 K m s −1 (approximately 285 W m −2 ). The subscript 0 indicates a surface value. The vertical potential temperature profile at the beginning features a constant value of 300 K up to a height of 1000 m. Above, a capping inversion with a gradient of 0.02 K m −1 is prescribed. For heights of 1300 m and more, a sponge layer (Rayleigh damping) is applied for all prognostic variables to reduce spurious reflections of vertically propagating waves from the model top. Wind velocities are initially set to zero everywhere because no background wind is considered. To accelerate the development of convection, random perturbations with a maximum amplitude of 0.25 m s −1 are imposed on the horizontal velocity field during the beginning of the simulation until a prescribed domain-averaged perturbation energy (or resolved-scale TKE) limit of 0.01 m 2 s −2 is exceeded. Besides, large-scale subsidence with a magnitude of up to 0.023 m s −1 guarantees a constant boundary layer height during the simulation. Otherwise, the dust devil statistics might change in time because of the dependency between the boundary layer height and dust devils (e.g., Hess and Spillane 1990;Fenton and Lorenz 2015). With these initial conditions, a quasi-stationary state of the convective PBL is simulated that mimics the reality during the afternoon of a sunny day.
The bottom and top boundary are regarded as impermeable (w = 0m s −1 ) with no-slip (Dirichlet) conditions (u = v = 0m s −1 ) at the ground and free-slip (Neumann) conditions (∂u/∂z = ∂v/∂z = 0s −1 ) at the top. For the perturbation pressure and the SGS TKE, a Neumann condition (∂ p * /∂z = 0 Pa m −1 , ∂e/∂z = 0 m s −2 ) is applied at the bottom surface. At the top surface, a fixed value of p * = 0 Pa is set, except for child domains, where a Neumann condition (∂ p * /∂z = 0 Pa m −1 ) is used. The potential temperature at the model top is derived locally by linear interpolation during each time step, utilizing the initial gradient at that height (horizontally homogeneous and constant during the run) and the temperature one grid point below, which is calculated through the prognostic model equation for the thermal internal energy. A further Neumann condition is set for e (∂e/∂z = 0 m s −2 ) at the top boundary. The nested domains obtain their top boundary conditions for the prognostic variables (θ , u, v, w) from their respective parent solution through a zero-order interpolation as described in Hellsten et al. (2021). In addition, a constant flux layer is assumed between the surface and the first computational grid point above. In this layer, unknown fluxes are calculated using  Monin-Obukhov similarity theory (MOST), which requires to specify a roughness length of the surface. We choose a value of 0.1 m (typical for rural areas). In the horizontal directions, cyclic boundary conditions are applied. Finally, the Coriolis force, although probably not that important for dust devils (e.g., Balme and Greeley 2006), is considered by setting the earth's angular velocity to 7.29 × 10 −5 rad s −1 and the latitude to 52 • .
The computational domain has a horizontal extent of 4 km × 4 km in order to resolve the large-scale polygonal convective cells. This is necessary because vertical vortices in the convective PBL are strongly tied to these cells (see Sect. 1). By the application of the vertical grid stretching above 1.2 km height, a vertical extent of the model domain of approximately 2 km is reached, which is well above the inversion layer. The child domains are about 240 m high, except for the simulation with the highest resolution (0.625 m), where a two-stage nesting (three domains) is used. Here, the second child domain has a height of 120 m. The selected heights are a compromise between computational costs and the demand to resolve the whole vertical thickness of a dust devil with fine resolution. The grid spacings in the root domains are 10 or 5 m along each spatial direction (except where the vertical grid stretching is applied) to capture the general flow statistics of the convective PBL well enough (see Sect. 1) and to limit the grid spacing ratio between parent and child domain to a maximum value of 5. This limitation is applied due to observations made by Hellsten et al. (2021) who state that the child solutions are almost independent of the chosen grid spacing ratio and that they all match to the non-nested case with the child resolution everywhere. Besides, our own test simulations show no significant difference between the dust devil statistics derived from an ensemble of 10 nested simulations with grid spacing ratio 5 (parent: 10 m, child: 2 m grid spacing) and the non-nested case, where a fine resolution of 2 m is used everywhere. In the nested set-ups of the convergence study, the child grid spacing is uniform in each direction and gradually reduced from 5 to 0.625 m (5, 2.5, 2, 1.25, 1, and 0.625 m), which results in 7 simulations each with a total simulated time of 4 h. An overview off all performed simulations is given in Table 2, stating the simulation name, the grid spacing, the domain size, and the number of grid points. The name is used as a reference in this study and is selected as follows: The first part indicates the root domain's grid spacing (e.g., R10). The second part, if any, describes the nested domain's grid spacing (e.g., N5) and simultaneously clarifies that a single-stage nesting (one child) is applied. If a two-stage nesting is used, a third part illustrates the resolution of the second nest and the first nest acts also as a parent domain. In the same row of the simulation name, the root domain characteristics are given. The first and second line below show the characteristics of the first and second child domain, respectively. Due to the grid stretching, the vertical extent of the root domains is rounded. Note that the dust devil data of a simulation are only taken from the nested domain with the finest resolution.

Detection of Vortex Centers and the Generation of Vortex Tracks
In the following, we will comprehensively clarify how a vortex center is identified and how various centers compose a vortex track, which we interpret as a dust devil-like vortex with a certain lifetime. The different tracks form the basis for our statistical analysis of dust devillike vortices (see Sect. 2.4). Our experiences show that the dust devil statistics can strongly depend on details of the detection and analysis algorithm, especially if a quantitative analysis is performed. In principle, the algorithm is based on the vortex detection and analysis approach from Raasch and Franke (2011) and Giersch et al. (2019), which we abbreviate with VDA11 (Vortex Detection and Analysis approach 2011) in the following. However, their approach has some drawbacks, which would heavily affect the comparison of results from runs with different resolution. For example, their vortices were detected on a horizontal plane taken from the first computational grid level above the surface at Δ/2. If the grid spacing changes, the physical detection height also changes, which will result in different statistics. Such a procedure would complicate or even prevent a quantitative comparison between differently resolved simulations. Therefore, a requirement for the algorithm in this study is a constant detection/analysis height, which we fix to 10 m or the next higher level, where a scalar grid point is located (e.g., 12.5 m height in the 5 m simulation). We select this height level to exclude potential effects of the SGS model on the statistics at least for high-resolution runs (only the first grid points above the surface are significantly affected by the SGS model (e.g., Schmidt and Schumann 1989;Gadde et al. 2021) and to enable comparisons with field measurements. Besides, dust devils are ground-level vortices and their vertical extent is often limited to several tens of meters (Balme and Greeley 2006), which is why the detection height is chosen close to the ground. Vortex centers are identified during the simulation after each model time step and after the model spin-up time of 45 min (see Sect. 3) by local minima of the dynamic perturbation pressure and local maxima of the absolute value of the vertical vorticity once certain thresholds are exceeded (see below). The vertical vorticity (hereafter only "vorticity") is defined as the vertical component of the rotation of the velocity field: The local minima/maxima refer to the lowest/highest value within a resolution-independent square of size 20 × 20 m 2 instead of just considering the eight neighbouring grid points as in VDA11. This square size reproduces typical sizes of large dust devils (e.g., Kurgansky 2006) and simultaneously guarantees comparability between all performed simulations of different resolution. For example, an area of 20 × 20 m 2 is the smallest possible square that can be defined for a local minima/maxima in R10. The exact position of the center is exclusively defined by the location of the pressure minimum and not by the location of the absolute vorticity maximum, which is a further difference to VDA11. This definition is more consistent with the calculation of the vortex radius that is determined from the tangentially (in circular direction) averaged pressure drop distribution around each center. The radial distance at which the absolute pressure drop is 50% of that at the center is interpreted as the (core) radius (see also Kanak 2005;Raasch and Franke 2011;Giersch et al. 2019;Giersch and Raasch 2021). This procedure agrees well with analytical and empirical models of dust devils, where the radius also often describes the location of the maximum tangential velocity (Lorenz 2014). The numerically calculated core radius can be used as an estimate for the visible dust column radius of real vortices.
Simultaneously to the pressure minimum, an absolute vorticity maximum must be located within a square of size 20 × 20 m 2 around a center. In this way, the positions of the maximum rotation and the pressure minimum can be slightly shifted, which can especially happen for high resolutions (see discussion about Fig. 3). The size of the square should not be too large. Otherwise, rotations in the flow field that do not belong to the same vortex structure might be assigned to it by mistake.
The vortex center is only considered for the generation of vortex tracks if certain pressure and vorticity thresholds are exceeded. Different values have been proposed in former studies that mainly depend on the model resolution or the detection height (Ohno and Takemi 2010a;Raasch and Franke 2011;Klose and Shao 2016;Nishizawa et al. 2016). These studies have in common that certain pressure and vorticity deviations from the horizontal mean are chosen for the thresholds. In our simulations, 3 times the standard deviation of the dynamic perturbation pressure and 5 times the standard deviation of the vorticity are chosen. The standard deviation is calculated as the arithmetic mean of the standard deviations from 4 instantaneous horizontal cross-sections at analysis height after 1, 2, 3, and 4 h simulated time. It must be taken into account that all simulations from Table 2 (except R10N2.5N0.625) are executed at least twice, once for calculating the correct thresholds and once for identifying the dust devillike vortices. Two simulations of R10N2.5N0.625 with 4 h simulated time are too expensive for us in terms of time and available computational resources (one simulation needs about 20 days wall-clock time on 6688 cores of an Atos/Bull system equipped with Intel Xeon Platinum 9242 processors) and thresholds are based on the first cross-section after the model spin-up time of 45 min. For the vorticity, the values 0.1, 0.24, 0.51, 0.63, 0.92, 1.08, and 1.56 s −1 are used, which indicates a continuous increase in local vorticity fluctuations with decreasing grid spacing, similar to other studies Giersch et al. 2019). The thresholds for π * show no clear trend and a grid-independent value of 3.5 Pa is taken. Note, the threshold values are derived empirically to capture most of the vortex lifetime, to eliminate the background noise of non-coherent turbulence, and to limit the data for the post processing, which especially includes the formation of vortex tracks (see next paragraph).
To form dust devil tracks, the detected centers are sequentially processed after the simulation in a grid-independent way. The processing steps are first described technically. Details for each step and the main differences to VDA11 are explained further below.
1. Centers of vortices with a radius r ≥ 100 m are neglected. 2. For finer resolutions than 10 m (all nested simulations), only the strongest center (rated by π * ) is considered if two or more centers are closer than 20 m to each other at the same time. 3. Two centers at different times are combined to the same vortex track if the displacement speed from the location of the center detected first to the center detected at later times (second center) is ≤ √ 2 × 10 m (corresponds to the distance of one R10 grid spacing in each horizontal direction) per 1.4 s (this time is defined by the mean time step of R10). If no second center is found after 4.2 s (three mean time steps in R10), a new dust devil track is initiated. 4. For simulations with Δ < 10 (all nested runs), centers of the same vortex can move up to 20 m in just one time step (see detailed explanation below), which would exceed the allowed maximum displacement speed by criterion 3. Therefore, two centers from a nested simulation that are detected at different time steps will be connected regardless of the third criterion if they have a distance to each other of 20 m or less even for very short time periods between the detections. 5. A further condition for connecting two vortex centers to a track is that their difference in π * must be less than 10% of the value of the first center. 6. Finally, the area-averaged (over 20 × 20 m 2 around the center's location) local vorticities must have the same sign to connect two vortex centers.
The first action, which was not part of VDA11, guarantees that only well-developed dust devil-like centers are connected. An investigation of the events with radii above 100 m has shown that such centers are weakly pronounced with values close to the detection thresholds so that they often vanish again quite quickly.
The second step avoids counting a single dust devil with several sub-vortices (see e.g., Bluestein et al. 2004) twice or even more and it takes care of the merging of centers, which is an important process for vortex intensification (Ohno and Takemi 2010a). Additionally, the comparability between differently resolved simulations is increased. Due to the definition of a center as a local extremum within a square of 20 × 20 m 2 , the 10 m run can not resolve situations where two or more centers are closer than 20 m (see Fig. 2, which shows the smallest possible distance between two centers A and B that might occur in R10). However, all nested simulations can principally resolve distances below 20 m as it is shown in Fig. 2 for R10N5 (center A and C). With the second criterion, we prohibit this technical discrepancy. The old algorithm VDA11 disregards all weaker vortex centers as soon as a center is within the radius of another (sub)vortex at the same time step. Taking the radius instead of a fixed value of 20 m is critical. Our results will show that the radius is strongly affected by the resolution. As a consequence, the detection and analysis method would indirectly depend on the grid, which we want to avoid as much as possible.
The third step allows the dust devil to occasionally have absolute vorticity or pressure drop values less than the thresholds. For example, consecutive centers of a dust devil track in R10 might have a temporal distance of three mean R10 time steps (3 × 1.4 s) if thresholds were not exceeded after the first and second time step and, thus, no suitable center was stored at these times. For the nested runs, the physical scales regarding time and space are maintained, i.e., the maximum allowed covered distance in 1.4 s is always √ 2 × 10 m, which corresponds to a displacement speed of about 10 m s −1 , and a period of 4.2 s is always scanned for a potentially subsequent center independent of the actual number of time steps that are needed to cover this time frame. For the highest resolution (R10N2.5N0.625), it follows that consecutive centers of a dust devil track are allowed to have temporal distances of more than 100 time steps (the mean time step in R10N2.5N0.625 is 0.038 s). In VDA11, two vortex centers were only connected if the position of the second center is not more than two grid points away from the first one after a maximum of three time steps. This procedure links the detection of dust devils to the model resolution. Therefore, we reject it.
Criterion 3 is completed by step 4 to allow a center displacement speed of much more than 10 m s −1 , which might be physically unrealistic but can happen technically anyway. Imagine a situation in the child of simulation R10N5 (mean time step of 0.75 s) where two centers of the same dust devil-like vortex with similar strength are closer than 20 m to each other (see Fig. 2 center A and C, distance is approximately 18 m). It might happen that at one time step center A is preferred and C is sorted out (based on criterion 2) but at the subsequent time step center C is preferred because it becomes stronger than A, which is then rejected. For this situation, the vortex center's displacement speed is 24 m s −1 (18 m/0.75 s) and, thus, it exceeds the maximum allowed value of 10 m s −1 according to point 3. Therefore, a new track would be generated by mistake instead of connecting center A (remaining center after first time step) with center C (remaining center after second time step). We enable this connection by step 4.
Finally, point 5 and 6 describe additional criteria that intend to avoid counting two centers from different dust devils to the same track and to prohibit unrealistically high short-term fluctuations of the dust devil features during its lifetime. The area-averaged local vorticity is taken instead of the local vorticity as used in VDA11because the area average is more robust with respect to the overall rotation direction of the dust devil-like vortex, especially for finer resolutions. To illustrate this, Fig. 3 shows two snapshots of the local vorticity in and around a cyclonic (a) and anticyclonic vortex (b) resolved with 10 and 1 m grid spacing, respectively. For 1 m, the structure is quite diverse and although the overall rotation is in clockwise direction some locations show a pronounced positive vorticity (anticlockwise spin), even in the central region. Therefore, it can happen that the value of the local vorticity attributed to a vortex center can change its sign within just two consecutive time steps. This would initiate a new dust devil track and centers would not be combined if local vorticity values are taken. However, the area averages of the local vorticities around the two vortex centers possess a clear positive and negative value (0.4 s −1 and −0.58 s −1 ).
The high positive vorticities in the clockwise rotating vortex can be better understood by transforming the vorticity into natural coordinates: with V the horizontal wind speed, n the normal direction to the velocity vector, and R the curvature radius (R > 0 for ⤿ and R < 0 for ⟳). Accordingly, the first term describes the . Center C is located at a grid point that only exist in the child of simulation R10N5, which is why it is visualized by a diamond. The line with two arrows highlights the distance of about 18 m between center A and C. Note, point A is used in the text as a reference location for a center that occurs in R10 as well as in the child of R10N5 rate of change of V normal to the direction of the flow (shear vorticity) and the second term describes the turning of the wind along a streamline (curvature vorticity). By analyzing highdensity streamlines and vector fields (see Fig. 3c, d), it turns out that the high positive vorticity values in the vortex center are a combination of both the shear and curvature vorticity. At these locations (e.g., at around x = -11m, y = 4 m, see white rectangle), the flow shows a cyclonic curvature (Fig. 3c) and also partly a decelerating flow along the normal direction ( Fig. 3d), resulting in an overall positive vorticity.

Statistical Vortex Analysis
The resulting vortex tracks from the formation algorithm above create a subset of an unknown dust devil population for which statistics shall be derived. We call this subset a dust devil sample. Each simulation creates such a sample but with a different resolution the underlying dust devil population changes. For example, high-resolution runs produce very short-lived vortices of a few seconds, while coarse-resolution runs do not capture them at all because of their higher time steps. This complicates a comparison between the statistics of differently resolved simulations. Therefore, it is mandatory to restrict the statistical analysis only to those vortices with a lifetime above a certain threshold. In this way, the dust devil population can be compared much better among the differently resolved simulations because all dust devil lifetimes that might occur in the population can principally be captured with every resolution. Additionally, we assume that very short-lived vertical vortices of several tens of seconds do not appear in nature as well-developed dust devils because of the small time frame for reaching a significant strength that initiates strong dust lifting. Therefore, they are usually not considered in any dust devil statistics which is based on observations. To increase the comparability with data from field studies and between simulations with different grid spacings, we restrict our dust devil population to those vortices with a lifetime τ of at least 120 s. The horizontal size of the vortices within the regarded population must also be constrained to allow for a good comparability between the differently resolved simulations. If the grid spacing is reduced, more and more smaller vortices occur (e.g., Giersch et al. 2019). For example, R10 allows a minimum radius of 10 m, while in R5N1 also vortices with radii of several meters can be simulated. Therefore, we limit the spatial vortex scale at the lower bound, i.e., the vortex must have a certain size and all other vortices smaller than this size are not considered. If the lower bound is selected high enough, a detection of vortices with spatial scales close to this bound would become possible for each applied grid spacing. In addition, limiting the statistics to larger vortices enables a proper resolution of their dynamics, which typically requires a minimum of 5-10 grid points along the dust devil's axes. In our study, the most obvious selection for the lower bound would be a minimum radius of 10 m because such a vortex could principally occur in each simulation, even though poorly resolved in the coarse-resolution runs. However, the minimum vortex radius must not be too large either. Otherwise, each sample contains only a few dust devil-like vortices from which a statistically profound derivation of vortex characteristics is impossible. Observational data suggest that a core radius of 10 m and more is already too large to get a sufficient number of dust devils during a simulated time of 4 h (Houser et al. 2003;Balme and Greeley 2006;Kurgansky 2006). They usually show radii of around 5 m. That is why we exclude all vortices with a lifetime-averaged radius lower than this value from our statistical investigation as a compromise between the above requirements for comparability and statistical analysis.
As explained above, our regarded dust devil population only considers vortex tracks with lifetimes of at least 120 s and lifetime-averaged radii of ≥ 5 m. We mainly calculate doubleaveraged quantities at the detection height of 10 m to describe this population statistically. The first average is executed over the vortex lifetime. Subsequently, an average over the sample is performed. In this way, each model run creates one value for one specific quantity. Instantaneous events of single dust devils or extreme values of a sample might be rare events and do not necessarily represent typical characteristics of a population. Taking only these events would make it difficult to compare results among the differently resolved simulations. That is why we focus our statistical analysis to mean values (if not otherwise specified).
The application of the sample average necessitates the first quantity that is analyzed: the sample size, which is expressed by the number of dust devil-like vortices N or the number density n in km −2 h −1 or km −2 d −1 . In addition, the statistical focus is on the sampleaveraged quantities π * , τ , r , the tangential, radial, and vertical velocity u tan , u rad , and w d ("d" for dust devil), respectively, the area-averaged vorticity ζ av and, finally, θ . Here, the overbar indicates a time-averaged value over the whole vortex lifetime. The quantities π * , ζ av and θ are defined in the vortex center. Instead, u tan , u rad , and w d represent the maximum of the tangentially averaged velocity distribution of the respective cylindrical component around the vortex center. Note that the radial and tangential components are calculated during the simulation through the transformation of the total Cartesian velocity components u and v to polar coordinates. Thus, radial and tangential velocities also contain the translational speed of the dust devil.
Based on the quantities above, the minimum lifetime of a vortex, which shall be considered in the statistics, is further restricted. An air parcel, which moves with the velocity u tan , must be able to circulate the vortex with a circumference of 2πr at least once during its lifetime. This is a further reasonable condition to focus the dust devil analysis only to well-developed vortices.

Results and Discussion
This section starts with an overview of general flow characteristics, which describe the physics of the simulated PBL in more detail. A focus is on the grid convergence of flow quantities that are typically used within studies of the PBL, like vertical profiles of the potential temperature or the friction velocity. In Sect. 3.2, dust devil statistics for variable grid spacings are analyzed with respect to their convergence. In addition, a quantitative comparison to observational data is performed. Finally, features of the three-dimensional structure of selected dust devil-like vortices are addressed without the claim to give a comprehensive description and explanation of the dust devil flow dynamics. A focus is again on how these features change with grid spacing.

General Flow Features
The general development of the flow in all root domains is very similar and can be evaluated through results from R10. A quasi-stationary state, where the turbulence statistics do not change substantially, is reached after 45 min as, for example, indicated by time series of the domain-averaged total kinetic energy E = 0.5 × (u 2 + v 2 + w 2 ) of the flow (see Fig.5a).
Here, the domain average refers to all heights up to the top of the (first) child domain (240 m) in every simulation to enable better comparability between the different model domains. The subsequent analysis only includes periods after the spin-up time of 45 min. Figure 4 shows horizontally and temporally averaged vertical profiles of the potential temperature (a), the total vertical turbulent heat flux (b), composed of the resolved-scale and subgrid-scale turbulent heat flux, and variances of v (c) and w (d). The horizontal average is marked by angular brackets, the overbar describes a temporal average over a period of 15 min before the respective output time, and the prime denotes a resolved-scale turbulent fluctuation, which is interpreted as the deviation of an instantaneous resolved-scale quantity from its horizontal domain average. The profiles reveal the typical characteristics of a convective PBL (see also Schmidt and Schumann 1989;Moeng and Sullivan 1993;Park and Baik 2014). The potential temperature indicates constant values in the so-called mixed layer and strong vertical gradients near the surface and in the entrainment zone at around 1 km, where the heat flux becomes negative (downward flux). The negative slope of < wθ > indicates negative divergences, causing a mean temperature increase in the PBL with time. In upper layers (≈ 1000-1100 m), overshooting thermals mix warmer air from the inversion layer downwards. This, together with the large-scale subsidence (see Sect. 2.2), results in a warming of the PBL and the free atmosphere during the simulation. Variances of the horizontal velocity components (only < v 2 > is shown because < u 2 > looks very similar) show stronger turbulence near the surface and in the entrainment zone generated by wind shear and buoyancy forces. Vertical velocity fluctuations reach a maximum at the lower third of the boundary layer height, which is constant over time (≈ 1 km ). In our study, the boundary layer height is defined as the point where the minimum of the total sensible heat flux profile is reached. Profiles in b, c, and d still show fluctuations around the actual mean state because the temporal average refers to a period of only 15 min (∼ one to two times the large-eddy turnover time defined by the ratio of the boundary layer height z i and the convective velocity scale w * , which is also referred to as the Deardorff velocity scale). Normally, several large-eddy turnover times are necessary to better capture the mean state. However, the shape remains similar.
The grid convergence of the general flow state is first investigated by time series of E, the horizontal average of the friction velocity < u * > calculated by means of MOST , and the maximum of the vertical velocity w max of the whole domain (see Fig. 5). The domain-averaged total kinetic energy represents the whole (resolved) flow, including a mean and turbulent contribution. However, because no mean wind is applied, E and the domain-averaged resolved-scale TKE of the flow are the same. All simulations show similar values for E that oscillate around a mean state of approximately 1.7 m 2 s −2 . Therefore, this quantity can be regarded as converged for every resolution.
The friction velocity, characterizing the turbulent momentum exchange near the surface, shows a grid spacing independent behaviour. This indicates that a resolution of 10 m is already fine enough to resolve the mean turbulent transport in the surface layer. According to Lyons et al. (2008), the ratio of the convective velocity scale and the friction velocity must amount to larger than 5 for dust devil formation. Our simulations indicate a grid-independent value of 2m s −1 (not shown) for the convective velocity scale and, thus, a ratio w * /u * of approximately 10.
In contrast to the friction velocity, the total maximum values of the velocity components, which occur somewhere in the respective domain and which represent extraordinary flow events, do not converge (the time evolution of u max and v max look very similar to w max ). The magnitude constantly increases with a decrease in the grid spacing, resulting in peak velocities of about 20 m s −1 for each component in simulation R10N2.5N0.625. This corresponds to about 10 times the convective velocity scale. In R5N1, maximum vertical velocities always occur near the center of a dust devil-like vortex after 1, 2, 3, and 4 h simulated time. The average height of the location of w max after the spin-up time is 30 m. In most cases, maxima occur below 20 m. Only in sporadic cases, the strongest updrafts in the child are simulated at heights between 100 and 240 m. Thus, we presume that maximum velocities are mostly caused by dust devils or the strong larger scale updrafts connected to them. Because w max values of up to 25 m s −1 have already been reported for dust devils (Balme and Greeley 2006), we assume that our u max , v max , and w max are still realistic and we speculate that a grid spacing of 0.625 m is just before the limit of what is required to show a convergence of maximum velocities. To clarify this point, a further simulation with even finer resolution would be required, which is currently beyond our available computational resources. However, the time series are a first indicator that extreme values in the dust devil statistics are not well suitable at this stage for evaluating a grid convergence (see Sect. 3.2).
In Fig. 6a, vertical profiles of the potential temperature show a constant increase of the surface temperature and the near-surface temperature gradients with a decrease in the grid spacing (the surface value is set to the value at the first computational grid level above but does not enter the prognostic model equations). With higher resolution, the thin super-adiabatic  2016), a critical value of 1-10 K m −1 is needed in the first few meters above the hot surface for the onset of dust devils. This has also been shown by observations (e.g., Oke et al. 2007;Ansmann et al. 2009). Such high gradients are only simulated in the simulations with child grid spacings of 1 and 0.625 m. As a consequence of increasing surface temperatures, stronger plumes are probably be generated, resulting in higher vertical velocities. This fits to the result that w max constantly increases with decreasing grid spacing (see previous paragraph). For heights of 20 m and above, all profiles of the potential temperature overlap, indicating a converged situation. This is also true for the potential temperature variance displayed in Fig. 6b. Very close to the surface, however, more intense temperature fluctuations are simulated if the grid spacing is reduced, which is directly related to stronger gradients in the mean profile. With such gradients, even small displacements of air parcels cause high temperature fluctuations.
As Fig. 6c shows, the mean variance of the u-component < u 2 > (< v 2 > looks almost the same) increases if the resolution changes from 10 to 5 m. A further reduction causes a significant increase below heights of 20 m only until a resolution of 1.25 m is reached. For 1.25, 1, and 0.625 m grid spacing, profiles almost overlap, even very close to the surface. This is confirmed by the resolved-scale TKE e res = 0.5 × (u 2 + v 2 + w 2 ) (see Fig.6d), which includes all variances in one quantity. For resolutions of 5 m and below, maximum differences among the profiles amount to several per cent only if heights above 20 m are regarded. At lower heights, only the profiles of R5N1.25, R5N1, and R10N2.5N0.625 match well.
Finally, the resolved-scale and SGS vertical turbulent heat flux is displayed in Fig. 6e, f (a double prime denotes SGS quantities). In each simulation, the first two vertical grid levels above the surface are significantly affected by non-resolved processes. At 20 m and above, however, SGS fluxes can be neglected for all resolutions because they contribute to the total flux only marginally (less than 10%). From a height of 30 m, all profiles starts to overlap and are converged. As discussed in Sect. 2.3, the detection height of dust devil-like vortices is set to 10 m or the next higher grid level. The heat flux profiles show that the turbulence is already well-resolved at that height for most of the simulations except R10. Therefore, uncertainties that might be introduced by the surface parameterization are mostly negligible in the subsequent dust devil statistics. The total vertical turbulent heat flux profiles overlap for all resolutions (not shown).
Dust devils are connected to the large-scale cellular pattern of the convective PBL (see e.g., Kanak 2006;Raasch and Franke 2011). Figure 7 reveals this pattern by horizontal crosssections of the instantaneous vertical velocity at 100 m height and it shows how structures depend on model resolution. With a resolution of 10 m, the large-scale polygonal structures are well-resolved. However, more flow details become apparent for R10N5 (Fig. 7b). A further reduction in grid spacing does not change the overall flow pattern, although smaller and smaller turbulent scales are captured but the displayed figure size in Fig. 7 does not allow to see them. The number of detected vortex centers at 4 h simulated time amounts to 33, 58, 178, 256, 254, and 261 for R10, R10N5, R10N2.5, R10N2, R5N1.25, and R5N1, respectively. Consequently, simulations with grid spacings of 2 m and below create a similar amount of centers. A more detailed discussion about the number of detected vortex centers and, thereby, dust devil-like vortices follows in Sect. 3.2.
All in all, our results support the findings from past studies about resolutions requirements for LES of the convective PBL (see Sect. 1). These studies mainly recommend a grid spacing on the order of 10 m. However, if the research focus is more on the surface layer, processes that originate from there, or details of the flow structures, 10 m is still too coarse. This is probably also true for dust devils, that are quite small-scale flow phenomena. Especially the above number of detected vortex centers and the resolution dependent variation of the near-surface temperature profile suggest a strong influence of the grid spacing on dust devils. We will investigate this in the next Sects. 3.2 and 3.3 in more detail.

Dust Devil Statistics
In the subsequent paragraphs, dust devil statistics at detection height are quantitatively analyzed for samples of a dust devil population, which only contains vortices with a lifetimeaveraged radius of 5 m or more and with a lifetime of at least 120 s. This limited population and the quantities that are analyzed have already been motivated and explained in Sect. 2.4. The statistics are usually based on one simulation with a certain grid spacing. However, to estimate statistical uncertainties, ensembles with 10 members have been created for grid spacings of 5 (R10N5) and 2.5 m (R10N2.5) by applying different random perturbations at the beginning of the respective simulation (see also Sect. 2.2). From these members, 95% confidence intervals have been derived for the mean and the standard deviation of the corresponding ensemble. Such an interval covers the true value with a probability of 95%. We take the same statistical significance interval as defined by (Appendix C Giersch and Raasch 2021) to assess if a value might be part of the same ensemble or not, i.e., as soon as a value lies outside the significance interval, differences are rated as statistically significant. A calculation of the confidence intervals for higher resolutions was not possible because the required Fig. 7 Horizontal cross-sections of the instantaneous vertical velocity at 100 m after 4 h simulated time taken from the simulations R10, R10N5, R10N2.5, R10N2, R5N1.25, and R5N1. Vortex centers at detection height are depicted as yellow dots ensemble runs were beyond our current computational resources. The statistical uncertainty ranges of R10N5 and R10N2.5 are very similar for most of the analyzed quantities (see Fig. 8), which is why we assume that they are also applicable to the other resolutions. Note, the exact detection height varies with the grid spacing due to the arrangement of the numerical grid. In case of 10 m grid spacing, dust devil centers are actually detected at 15 m and for 1 m resolution at a height of 10.5 m. This creates a systematic uncertainty to higher or lower values. However, with decreasing grid spacing the analyzed grid level approaches the physical height of 10 m and the systematic bias is reduced.
The dust devil characteristics and their dependencies on the grid width are shown in Fig. 8. The number of detected dust devil-like vortices (Fig. 8a) varies between 700 and 2400, which corresponds to number densities of 13 km −2 h −1 and 46 km −2 h −1 , respectively (the dust devil detection time during the simulation is 3 h and 15 min). Therefore, the order of magnitude for n is assumed to be 10 km −2 h −1 . The deviations between coarser resolutions (10-−2.5 m) are mainly provoked by the definition of the regarded dust devil population. The smallest vortices with radii between 5 and 10 m, which occur more frequently than larger vortices, can hardly exist in R10, R10N5, and R10N2.5, if at all. The finer the resolution the more of these vortices are resolved, which can be seen in Fig. 9, where the radius data is grouped into bins with an equal size ratio of about √ 2. The maximum number moves towards smaller radii and increases. For R10, R10N5, and R10N2.5, it is located at the bins [15;21.63), [7.21;10.4), and [5;7.21), respectively. A further reduction in grid spacing (2.5-−0.625 m) causes no further increase of N because dust devil-like vortices with radii smaller than 5 m are neglected in our population. Instead, the total number decreases. For grid spacings of 2 m or lower, Fig. 9 and observational data by Oncley et al. (2016) suggest that the maximum of the number distribution of r is not captured anymore by a population, which neglects vortices smaller than 5 m. The maximum moves out of the considered radius range if the grid spacing reduces. According to the significance intervals from simulation R10N5 and R10N2.5, we assume a converged value of about N=1000, corresponding to 19 km −2 h −1 or 77 km −2 d −1 if a typical sunny day allows for 4 h of strong dust devil activity (see also Lorenz 2014). Optical detections during field experiments state frequencies between 0.1 and 800 km −2 d −1 , depending on the survey area (Balme and Greeley 2006;Lorenz 2009;Lorenz and Jackson 2016). Lorenz (2009) proposes the formula n ∼ 50/A, indicating that the number density in km −2 d −1 is inversely proportional to the survey area A in km 2 . An application of this formula to our study results in n ≈ 3 km −2 d −1 , which is less than the simulated values. However, according to Lorenz (2014), approximately 100 dust devil counts per square kilometer and day is the most likely formation rate of visible dust devils under favorable meteorological conditions. The simulated value of 77 km −2 d −1 is quite close to this rate. It must be noted that not all of the simulated vortices would be visible in nature because pressure and vorticity detection thresholds in the numerical simulations correspond to intensities much lower than the values that would be needed for dust lifting. If we additionally assume a threshold core pressure drop of 30 Pa for the occurrence of dust lifting (Lorenz 2014), only 200-300 vortices occur in simulation R5N1.25, R5N1, and R10N2.5N0.625 that reach such a high value at least once during their lifetimes. This is approximately equal to a frequency of 20 instead of 77 km −2 d −1 , illustrating how an apparently easy statistical quantity like the number density can fluctuate strongly depending on the details of the investigation approach.
A typical lifetime for the investigated dust devil population varies between 240 and 310 s (Fig. 8b). Maximum values fluctuate between 1500 and 2600 s. However, a convergence can not be identified. The changes from R5N1 to R10N2.5N0.625 are still statistically significant if error bars with comparable sizes to R10N5 or R10N2.5 are applied to the data of R5N1 and R10N2.5N0.625. Additionally, the results indicate that there is no benefit from performing  Lorenz and Jackson (2016). The displayed values on the x-axis mark the lower and upper bound of a bin interval, e.g., [5,7.21) for the first one high-resolution and costly simulations if the mean lifetime shall be estimated. No clear trend in a certain direction is visible. Observational data show that most dust devils last for only a few minutes (Lorenz 2013). In rare cases, a duration of several tens of minutes and even several hours is possible (Balme and Greeley 2006). Consequently, the simulated lifetimes match the observed range for each simulation. The absence of dust devil-like vortices with a duration of several hours might be attributed to the rareness of such events in combination with the idealized setups (e.g., lack of large-scale vorticity or limited simulation time). Note, lifetimes derived from field studies mostly refer to the time where the vortex is visible. Because our detection thresholds are rather low, simulated vortices would probably not be visible in nature during the whole lifetime. Therefore, our values might have a bias towards larger lifetimes compared to observations. Figure 8c displays the mean radius of each sample and Fig. 9 indicates how many vortices belong to a certain size class. It is apparent that coarse resolutions strongly overestimate the horizontal dust devil size. A smaller grid spacing enables the resolution of smaller vortices, which is why the mean radius reduces with a decrease in the grid spacing. Because the defined population is limited to vortices with radii of 5 m and more, the radius reduction slows down for smaller grid widths and a converged value between 6 and 8 m arises at resolutions of 2.5 m and below. Kurgansky (2006) indicates a highly skewed size distribution with characteristic radii between 0.85 and 4.15 m, depending on the investigated field study. Balme and Greeley (2006) provide dust devil radius frequency distributions derived from six different data sets. A mean radius of 3.5 m was derived from the observations with a distribution skewed towards the smaller sizes (positive skewness). Thus, our simulated and converged vortex size of 6-8 m is larger compared to the findings in nature. This is probably because we cut off too much of the real dust devil population at the smaller sizes if we restrict our analysis only to those vortices that have at least a 5 m radius. An analysis of the population where each radius and lifetime is allowed reveals a constant but diminishing decline to a mean radius of approximately 3.2 m in simulation R10N2.5N0.625 (not shown), which fits much better to the observed range. As indicated by Fig. 9, the largest vortices in our simulations reach lifetime-averaged radii between 20 and 40 m.
The mean potential temperature at detection height (Fig. 8d) tends to increase with a reduction in grid width, especially for coarser resolutions of the investigated range. The development hypothesizes a converged value slightly above 306 K. This corresponds to a mean horizontal temperature difference between the dust devil core and the surroundings of about 3.5 K if the profile data from Fig. 6a are taken as a reference for the ambient conditions. The mean vertical temperature gradient (difference between the mean surface value and the mean value at analysis height) constantly increases with finer resolution and amounts to 0.13 K m −1 for Δ = 0.625 m. According to Balme and Greeley (2006), measured horizontal temperature excursions of less than 10 K are common. The total possible range is quite large and stated as 1-20 K. Sinclair (1973) did some measurements at 31 ft (similar to our analysis height) that indicate typical excursions of 2-4 K. Data derived from thermal image velocimetry show a temperature difference of up to 3 K for a single dust devil (Inagaki and Kanda 2022). However, measurement studies often refer to the maximum recorded temperature excursion during the dust devil's lifetime instead of lifetime-averaged values. If we consult our maximum values and average them over the sample size, they are 1-2 K higher than the sample mean of the lifetime-averaged temperatures. In R10N2.5N0.625, one dust devil even shows an instantaneous temperature of nearly 311 K, corresponding to a deviation of more than 8 K. All this suggests that our simulated values fit quite well to reality, at least for resolutions below 2 m. Coarser grids underestimate the temperature jump within dust devils significantly because the super-adiabatic layer close to the ground from which the heat is sucked into the dust devils core is poorly resolved (see also Fig. 6a). Figure 8e demonstrates the development of the mean absolute pressure drop in the dust devil's core with the grid spacing. In accordance with Giersch et al. (2019), the strength increases with decreasing grid spacing. From a resolution of 1.25 m, the changes are insignificant and a converged value of approximately 14 Pa is reached. A fair comparison to measurements is only possible with the peak pressure drop of a dust devil-like vortex because this is the quantity that matches the most to what is reported in field studies. For finer resolutions (1.25, 1, and 0.625 m), single dust devils show maximum pressure drops between 200 and 300 Pa. Sample averages of the maximum values vary between 25 and 35 Pa. The coarser resolutions (5 and 10 m) only simulate maximum values below 120 Pa (sample averages between 10 and 20 Pa). Unlike the lifetime-averaged pressure drops, maximum values constantly increase and show no convergence similar to what has been shown in Fig. 5 for the peak vertical velocity of the model domain. In reality, intense and visible dust devils show pressure excursions of several 100 Pa (Balme and Greeley 2006), which is similar to the results from the high-resolution runs. Pressure measurements recorded by Lorenz and Lanagan (2014) disclose peak pressure dips of convective vortices that may or may not be dust-laden in a range of 20-150 Pa, which suggests smaller intensities compared to Balme and Greeley (2006). This might be due to the fact that Lorenz and Lanagan (2014) derived the pressure drops from fixed stations that were not necessarily inside the dust devil core and from dust devils that were not necessarily dust-laden. All in all, the simulated values for resolutions below 2 m agree well with the reality. Otherwise, pressure drops tend to be underestimated. Our numerical simulations would probably produce even stronger vortices if heterogeneities and background winds are considered (Giersch et al. 2019).
Another measure beside the pressure drop that is used to evaluate the vortex strength is the vertical vorticity. Figure 8f visualizes the mean absolute value of the area-averaged vorticity around the center (for the definition see Sect. 2.3) and how it depends on the resolution. Similar to the lifetime, no definite trend is observable and values fluctuate around 0.23 s −1 , which we define as our best vorticity estimate of the regarded population. The sample mean of the maximum area-averaged vorticity during a dust devil's lifetime varies between 0.3 and 0.4 s −1 , depending on the resolution. Single vortices with a very strong rotation even reach instantaneous area-averaged vorticities of approximately 1 s −1 . Higher resolutions tend to produce higher peak values. References to measured vorticities within terrestrial dust devils do rarely exist. Doppler radar measurements by Bluestein et al. (2004) show local core vorticities that range from 0.5-1 s −1 , similar to observations by Oncley et al. (2016). Temporally and tangentially averaged data for one single dust devil show a peak vertical vorticity of 1.8 s −1 in the core region, which reduces nearly linearly to approximately 0 s −1 at the core radius of about 12 m (Inagaki and Kanda 2022). All these values indicate typical vorticities of ∼ 1 s −1 within the dust devil core. If the vorticity is additionally averaged in a horizontal plane of 20 × 20 m 2 around that core, similar results as our estimate of 0.23 s −1 are expected. However, for getting this estimate, high-resolution runs are not mandatory.
The tangential velocity in Fig. 8g shows a similar pattern compared to |ζ av |. This indicates that the horizontally averaged vorticity around a dust devil-like vortex strongly correlates with the maximum tangential velocity, occurring at a certain distance apart from the center. The pressure drop do not show this correlation, which could be expected from the cyclostrophic balance. It might be that the pressure drop in the center has less effects on the maximum tangential velocity but rather on the overall rotation, which could be better described by integral quantities like the circulation or the horizontally averaged tangential velocity. However, these results are contrary to observations made by ( Fig.9f Oncley et al. 2016), which showed a fairly well cyclostrophic balance for the maximum tangential velocity but a bad one for the averaged tangential velocity. This mismatch to our results might be caused by the different methodologies used to define a vortex. For example, Oncley et al. (2016) fitted the circular vortex structure by eye, whereas we use the pressure distribution to determine the vortex scale (see again Sect. 2.3). Consequently, completely different horizontal vortex surfaces could be determined for the same vortex, which significantly influences the calculation of maximum and averaged tangential velocities. Similar to Oncley et al. (2016), we additionally argue that the radial velocity causes a significant deviation from the cyclostrophic balance and, thus, disturb a potential correlation between u tan and |π * |. Fluctuations of u tan appear in a range of 2.05-−2.20 m s −1 with no distinct development to higher or lower values. Therefore, very fine resolutions are not beneficial for this measure. Instead, instantaneous peak values of 13 m s −1 are simulated in R10N2.5N0.625, whereas simulation R10 reveals much smaller peak tangential velocities (≈ 5 m s −1 ). This is also true for the sample averages of the maximum tangential velocity during a vortex lifetime that range from 2.5 to 3.5 m s −1 , depending on the resolution. Lower grid spacings simulate higher maximum velocities. Based on several measurement studies, Balme and Greeley (2006) stated that the peak tangential component of the wind speed usually reaches 5-10 m s −1 . In extreme cases, up to 20 m s −1 are possible. Likewise, observations from a fixed array of 31 turbulence sensors demonstrate a maximum value of tangential velocity of 8.9 m s −1 (Oncley et al. 2016). Stull (1988) specifies that tangential velocities are on the order of 10 m s −1 . Novel measurements with high spatial and temporal resolution show a tangential velocity component of up to 4.2 m s −1 for a single dust devil (Inagaki and Kanda 2022). Because this value represents an averaged velocity along the circular direction and over 40 s it appears to be smaller than the other measurements mentioned before. In summary, our numerical simulations with a grid spacing of ∼ 1 m resemble the tangential velocity more realistic compared to, for example, the 10 m run.
The radial velocity converges to 1.35 m s −1 (see Fig. 8h) if grid spacings below 2 m are applied. A magnitude of about 1.35 m s −1 fits very well to the temporally and tangentially averaged measurements from Inagaki and Kanda (2022), which show radial velocities between 1 and 2 m s −1 along the radial direction of a single dust devil-like vortex. A typical maximum radial velocity that is simulated during a dust devil's lifetime and for resolutions below 2 m is 2 m s −1 , while higher values are reached for higher resolutions. In extreme events and only for resolutions below 2 m, peak radial velocities up to 4 m s −1 are simulated at certain times. Based on measurements, Sinclair (1973) derived typical maximum radial velocities of ∼ 5 m s −1 , similar to Kaimal and Businger (1970) whose time series suggest maximum radial velocities between 3 and 6 m s −1 , depending on the measured height. Note, the surface roughness as the main reason for the radial velocity and the deviation from the cyclostrophic balance is assumed as 0.1 m in our simulations. A different roughness would produce different values that might fit better or worse to observations. The vertical velocity component shows the smallest significance intervals of all investigated quantities with respect to the changes that occur due to a modification in the grid spacing. Over the entire resolution range, the lifetime-averaged vertical velocities almost double (from 1 to 2 m s −1 ) with higher values at small grid spacings, whereas the relative statistical uncertainty expressed by the significance interval is just a few per cent. No convergence is reached for the vertical component. This might be related to the potential temperature, which also constantly increases for higher resolutions (see discussion above). The buoyancy and, thus, the strength of updraughts is determined by horizontal temperature differences. The profile data in Fig. 6 indicate a grid-independent reference temperature at 10 m height. If the core region of a dust devil-like vortex becomes warmer for higher resolutions at the same height, horizontal temperature differences increase on average and, consequently, updraughts become stronger. The maximum vertical velocities during the vortex lifetime show sample averages of about 4-5 m s −1 for 1.25, 1 and 0.625 m grid spacing with peak values of 20 m s −1 . A constant increase with decreasing grid spacing is also evident for the sample-averaged maximum values. The coarsest resolution run R10 reveals a value of only 1.5 m s −1 . Comparisons with field studies suggest that fine resolutions below 2 m are necessary to capture the vertical component realistically. Using a mobile instrumented tower, Sinclair (1973) measures maximum vertical motions of ∼ 10 m s −1 in all investigated dust devils and at all regarded height levels (7 (2), 17 (5), and 31 ft (9 m)). In-situ wind speed measurements by Kaimal and Businger (1970) also show roughly height independent magnitudes of the w-component. Peak values of 3 to 4 m s −1 are reported, similar to measurement data acquired by Fitzjarrald (1973) and Tratt et al. (2003). According to Balme and Greeley (2006), typical vertical wind speeds are less than 10 m s −1 .
In addition to a comparison with measurements, dust devil data can also be compared to theoretical models like the thermodynamical scaling theory of Rennó et al. (1998). In this theory, the pressure drop and the maximum tangential wind speed across a dust devil can be approximated by: and: v m ≈ γ ηc p ΔT .
The variable γ describes the fraction of the total dissipation of mechanical energy consumed by friction at the surface with typical values between 0.5 and 1 Rennó et al. (1998), η is the thermal efficiency of a heat engine, c p = 1005 J kg −1 K −1 is the specific heat capacity at constant pressure, R = 287 J kg −1 K −1 is the gas constant of dry air, p ∞ is the surface pressure (101,325 Pa in our simulations), T ∞ is the absolute temperature of the ambient air outside the dust devil determined by the profiles of Fig. 6a, and ΔT is the effective temperature perturbation. According to Souza et al. (2000) and Kurgansky et al. (2016), η Table 3 Comparison between the dust devil scaling theory proposed by Rennó et al. (1998)  can be calculated via: where g = 9.81 m 2 s −1 is the gravitational acceleration. We use the simplified formulas suggested by Kurgansky et al. (2016). Because this model describes the order of magnitude of the maximum values, we apply it to the peak values of our simulated dust devil statistics at 10 m height. Therefore, ΔT = T m − T ∞ , with T m the peak absolute core temperature of every simulated dust devil sample. Similarly, |π * | m and u tan m describe the simulated pressure drop and tangential velocity maxima for each sample. Table 3 shows the results. Only for resolutions of approximately 1 m or lower, the simulated peak values (|π * | m , u tan m ) match the range of the theoretical values (Δp, v m ). Otherwise, our simulations underestimate the magnitude, which is in agreement with the comparison to observational studies performed before.
The above discussion demonstrates that vortex properties can vary significantly if the grid spacing is changed and that the concrete resolution dependence differs between the regarded quantities. So far, no general answer can be given to the question at what resolution the overall statistics are converged. However, resolutions below 2 m show a convergence for most of the analyzed mean quantities and results fit very well to observations and measurements of real dust devils. Comparisons of the strongest dust devil events with observations and the thermodynamical theory of Rennó et al. (1998) suggest that a resolution of 0.625 m is short before the minimum grid spacing that is required to reach a convergence even of the peak values. This section also indicates that a quantitative comparison between numerical simulations and field experiments must be carefully performed. Different analysis heights might be relevant and lifetime-averaged values must not be mixed with maximum values. Additionally, huge differences of several orders of magnitude partly exist among various measurement data, which makes a comparison to numerical results quite challenging. Also, sometimes only visible, dust-laden vortices are considered in the statistics. In other cases, a pronounced pressure drop is enough to define a dust devil-like vortex. Therefore, it is mandatory to clearly define the characteristics of the regarded dust devil population. Changes in the considered population can cause significant changes in the results. Finally, different boundary conditions (e.g., roughness length, terrain type or slope) and meteorological conditions (e.g., background wind, boundary layer height, heating rate) complicate a comparison between numerical simulations and measurements. Nevertheless, it can be summarized that with grid spacings of less than 2 m dust devils with radii of more than 5 m are detected with sufficient accuracy. In the following, the three-dimensional flow structure is investigated regarding grid spacing changes to review the above resolution suggestion of less than 2 m. However, a detailed and sophisticated flow analysis is not intended.

Three-Dimensional Structure
The study of the grid convergence of the three-dimensional vortex structure is realized by the analysis of instantaneous as well as time-averaged horizontal and vertical cross-sections of selected vortices. Only data from the most durable dust devil-like vortices are considered. These dust devils can be regarded as representative for well-developed and pronounced vertical vortices, occurring in the regarded population. The most persistent vortices in R10, R10N5, R10N2.5, R10N2, R5N1.25, R5N1, and R10N2.5N0.625 are tracked for 1829, 1590, 2549, 2207, 2203, 1838, and 1628 s, respectively. The sampling of the vortices' variables ( p * , θ , ζ , u, v, and w) uses the same algorithm as presented in Raasch and Franke (2011). It is based on a three-dimensional grid defined in the vortices' centers that moves together with them during their whole lifetime. Instantaneous (at every time step) as well as time-averaged data is stored after the simulation. The averaging procedure is performed during the model run. For each simulation, it is guaranteed that the output volume has at least an extent of 140 × 140 × 100 m 3 (see cross-section size in Figs. 12 and 13). This is a compromise between available storage space, memory, and to ensure that the dust devil's main sphere of influence is recorded. Note that two identical simulations are required for the three-dimensional analysis. Based on the first run, the dust devil to be examined is identified as well as its track. The second run uses the center coordinates from the first run to perform the sampling. A rerun of R10N2.5N0.625 was not possible due to the high computational demand associated with it. Therefore, only results for grid spacings between 10 and 1 m are discussed.

Instantaneous Data
This section provides an overview of instantaneous and short-term vortex features that are not visible in time-averaged fields. The main focus is on how these features change with the resolution. Our data reveal a frequent interaction of approaching vortices, independent of the grid spacing. The interaction can result in an intensification, dilution, or maintenance of the original vortex. Figure 10 shows an example of how a dust devil-like vortex is maintained and intensified by absorbing another vortex with the same sense of rotation (positive vorticity). Between about 6833 and 6866 s simulated time (Fig. 10a, b), the main vortex in the middle shows perturbation pressures between -37 and -44 Pa, whereas the second and weaker vortex on the right has minimum pressures in the range of -20 to -30 Pa. After the merging, only one strong vortex remains with a maximum pressure drop of more than 50 Pa (Fig. 10d). In Sect. 3.2, the clustered occurrence of dust devils along convergence lines and near vertices of the cellular flow pattern has already suggested a strong interaction between different vortex centers during their lifetimes (see Fig. 7). The results are in agreement with Ohno and Takemi (2010a), who pointed out that most of the strong dust devils are first intensified through the merger of multiple vortices and subsequently maintained and more enhanced by additionally incorporating small-scale vortices. Doppler radar observations by Bluestein et al. (2004) also support a frequent interaction between different dust devil-like vortices.
The decrease of the grid spacing discloses more and more instantaneous fine-scale flow features. Some of them are displayed in Fig. 11. Panel a and b show horizontal snapshots of the vertical velocity and the horizontal wind at detection height for the analyzed vortices taken from simulation R10 and R10N2, respectively. Beside the overestimation of the vortex size in Fig. 10 Horizontal snapshots of the perturbation pressure of the most durable vortex that occurs in simulation R10N2.5 at detection height (z = 11.25 m) after approximately 6833, 6866, 6897, and 6930 s simulated time. Arrows indicate horizontal wind vectors. Every fourth vector is shown R10 (see also Sects. 3.2 and 3.3.2), downdraughts are apparent in both core regions surrounded by pronounced positive vertical velocities. Such downdraughts appear to be an inherent feature for all investigated fully developed dust devils, which agrees well with measurements (Balme and Greeley 2006). However, vortices in R10 rarely show these central downdraughts, whereas vortices resolved with finer resolutions reveal a descending motion in or close to the center for most of the time. In the time-averaged dust devil data, upward motions are dominant nearly everywhere (see Sect. 3.3.2). This is not contradicting to the instantaneous data because the location of the region of descending air changes during the vortex lifetime. Thus, positive and negative vertical velocities alternate in the center (defined by the pressure minimum), still resulting in a positive time-averaged value. Horizontal gradients of w are much more distinct for the smaller grid spacings. The vortex in R10N2 partially shows an increase in the vertical velocity from -5 to 6 m s −1 over a distance of just 10 m. Observations related to dust devils have already shown the existence of pronounced downdraughts within the core (Kaimal and Businger 1970;Sinclair 1973). The two-cell vortex concept with descending air close the axis of rotation and upward motions aside (e.g., Mullen and Maxworthy 1977) fits also very well to our numerical simulations. Figure 11c indicates an extended low pressure region southwest to the dust devil center (x= [-30;0], y= [-40;-20]), which is another distinct flow feature that arises in our simulations independent of the resolution. We refer to this region as a tail-like structure in the following. The most striking pattern, however, is visible in the vertical velocity. Strong descending air is separated from strong upward motions in the same region where the tail in the pressure data is visible. The updraughts are always located closer to the center than the downdraughts. As far as we know, this dust devil characteristic has not been reported so far, probably because of missing high-resolution three-dimensional data. Nevertheless, a comprehensive analysis is not intended here and goes far beyond the scope of this study.
Finally, Fig. 11d indicates that dust devils can momentarily have several sub-centers that appear, for example, through neighbouring local pressure minima. This supports previous findings of secondary vortices within or around a main dust devil (Bluestein et al. 2004;Zhao et al. 2004;Oncley et al. 2016). Especially Zhao et al. (2004) highlight this issue in their dust-devil-scale simulation with a resolution down to 0.1 m. Pressure contours showed up to eight sub-centers that occur along the annular zone around the center where the strongest radial shear in both the tangential and the axial velocity components exist. In the most persistent vortex from R5N1.25, up to three distinct pressure drops appear next to each other. The number of secondary vortices scales with the resolution. A maximum of two sub-centers is visible in the dust devil data of simulation R10N2.5. For a grid width of 1 m, up to 4 sary vortices are observable. However, they seem to be quite unstable and last not more than several tens of seconds, which is why they can not be seen in the time-averaged fields.
The above discussion shows again that the model resolution is one of the most relevant parameters with respect to LES of dust devils. A definite recommendation which grid spacing must at least be used to capture the instantaneous vortex characteristics realistically is difficult to derive because appropriate field data that contain instantaneous, three-dimensional vortex information with high spatial and temporal resolution rarely exist or are limited to a few dust devil properties only (e.g., Bluestein et al. 2004;Oncley et al. 2016;Inagaki and Kanda 2022). Also the exact research question of follow-up studies will determine the resolution requirements. Our simulations with grid spacings of 10 and 5 m are missing some dust devil features like sub-vortices and they overestimate the sphere of influence significantly. This is why we suggest a minimum resolution below 5 m to qualitatively capture the instantaneous dust devil behaviour realistically.

Time-Averaged Data
A temporal average of quantities within the volume around the vortices' centers reveal typical dust devil features that have already been reported in other field or LES studies (for a comprehensive overview see Reiss et al. 2017). As indicated by horizontal cross-sections in Fig. 12, the surrounding region of a dust devil is characterized by negative pressure excursions with maximum values of the normalized perturbation pressure p * / p * min close to the center. The variable p * min describes the minimum pressure drop that occurs in the whole region of 140 × 140 m 2 . Arrows of the horizontal wind vector show that the horizontal convergent flow rotates the strongest just outside the central region. Directly around the center, a flow with a weak radial component pointing away from the core is partly visible. This was already noticed by Sinclair (1973), who stated that there might be a radial outflow within the dust column. Also, Balme and Greeley (2006) illustrate this issue in their Fig. 8. The potential temperature, the vertical vorticity and velocity show distinct maximum values in or close to the center with a pattern similar to the one in Fig. 12. The local vertical vorticity is even Fig. 11 a-c Display horizontal snapshots of the vertical velocity for the most durable vortices in simulation R10, R10N2, and R10N2.5. d Perturbation pressure for the most persistent vortex in simulation R5N1.25. Isolines with 8 Pa distance (c) and 5 Pa distance (d) highlight the distribution of the perturbation pressure. Arrows indicate horizontal wind vectors. Every fifth vector is plotted in (b). In the (c, d), every third vector is shown. Note the different color scales, reference arrows, and spatial ranges that are plotted for improving readability more concentrated in the vortex core than the pressure drop and starts to randomly fluctuate around zero if a distance of more than 20-30 m to the center is reached (not shown). The mean vertical velocity maxima are mostly located several grid points adjacent to the vortex core similar to (Fig.8 Raasch and Franke (2011)). Independent of the regarded quantity, all cross-sections demonstrate that coarser grid resolutions overestimate the horizontal vortex size and, thereby, the vortex sphere of influence significantly. Grid resolutions of 2.5 m or smaller (Fig. 12c-f) do not alter this sphere anymore, which is in agreement with the more quantitative analysis in Sect. 3.2 (see Fig. 8c).
Temporally averaged vertical cross-sections of the normalized temperature difference (θ − θ min )/(θ max − θ min ) in Fig. 13 indicate the highest temperatures at ground level in the vortex core. The quantities θ max and θ min depict the maximum and minimum temperature occurring in the analyzed region of 140 × 100 m 2 . Δθ is defined as θ max − θ min . Temperature gradients tend to increase with lower grid spacings, which is a consequence of the better resolution of the super-adiabatic layer close to the surface. Throughout the whole displayed vertical range, warmer temperatures are apparent above the central region compared to temperatures outside the vortex core at the same height level. This is caused by the predominant vertical flow, which lifts the super-adiabatic surface layer of warm air. In agreement with the results Fig. 12 Time-averaged horizontal cross-sections of the normalized perturbation pressure at detection height derived from different dust devil-like vortices that occur in R10, R10N5, R10N2.5, R10N2, R5N1.25, R5N1, and R10N2.5N0.625. Arrows indicate horizontal wind vectors. Their distance to each other is always 10 m from Sect. 3.1, the vertical extent of this layer outside the dust devil core is simulated to be very small for resolutions of around 2 m and below, but expands to several 10 m in height after being sucked into the vortex. A downward flow, which has often been reported (Kaimal and Businger 1970;Sinclair 1973;Balme and Greeley 2006), does not occur in the mean fields. However, reduced positive vertical velocities are visible in the vortex core compared to the values next to the center (not shown), which is in agreement with Raasch and Franke (2011). This is caused by instantaneous central downdraughts (see Sect. 3.3.1). Additionally, strong radial inflow is limited to the lower vortex regions. Air parcels that move near the surface heat up due to the prescribed positive surface heat flux and reach their highest temperatures when they approach the updraft region. Again, the size of the vortices is overestimated for resolutions larger than 2.5 m, which is also supported by vertical cross-sections of ζ /|ζ | max and p * / p * min (not shown). Isolines of the perturbation pressure in Fig. 13 reveal a minimum of the pressure at about 10 m height but only for grid spacings of 2 m and less. Such a pronounced and lifted minimum does not occur for coarser resolutions. Because DNS simulations of dust devils in Rayleigh-Bénard convection have also shown this feature (Giersch and Raasch 2021), we interpret it as an inherent characteristic of dust devil-like vortices. According to our knowledge, neither LES studies nor case studies of real dust devils have shown a lifted pressure minimum yet, probably because of a poor resolution or poorly available three-dimensional measurement data. A simple assumption of a steady circular vortex in cyclostrophic balance can explain this finding. Due to surface friction, maximum rotational velocities occur above the surface and not directly adjacent to it, which results in maximum radial pressure gradients at the same height level and, thus, stronger pressure drops. This situation is first captured for a resolution of 2 m. Also the mathematical model of whirlwinds by Pandey and Maurya (2017) includes a negative pressure gradient along z, which is especially important for the whirlwind to grow vertically.
To conclude, the mean three-dimensional structure of single dust devils reaches a nearly converged state at resolutions of approximately 2 m. For grid spacings above 2 m, especially the high spatial gradients and the vortex size are not captured appropriately. The overall convergent and spiraling upward flow is realistically simulated for each selected resolution, at least from a purely qualitative perspective.

Summary and Conclusions
In this study, we numerically investigated atmospheric dust devil-like vortices and their statistical properties. We focused on the resolution dependent convergence of various dust devil parameters. Simulations were performed with the large-eddy simulation (LES) model PALM. By using the nesting feature of PALM, we explored grid spacings between 10 and 0.625 m within a domain of approximately 4 × 4 × 2 km −3 . This domain and resolution captured the large-scale cellular pattern of the convective planetary boundary layer (PBL), in which dust devils naturally form.
As a first step, we developed a revised and resolution-independent version of the dust devil detection and analysis algorithm of Raasch and Franke (2011) and Giersch et al. (2019). It improveed the comparability of dust devil statistics and properties between simulations of different resolutions and facilitateed a direct comparison with field measurements and observations. We showed how careful this algorithm must be designed to adequately capture the natural spatial and temporal vortex scales and highlighted the challenges to create compa- Fig. 13 Time-averaged vertical cross-sections of the normalized temperature difference through the center (y = 0) derived from different dust devil-like vortices that occur in R10, R10N5, R10N2.5, R10N2, R5N1.25, R5N1, and R10N2.5N0.625. Arrows indicate vertical wind vectors derived from u and w. Their distance to each other is always 10 m. Isolines with a distance of 8 Pa show the perturbation pressure rable results between model runs with different resolutions. Also, the analyzed dust devil population for which the derived statistics are valid, must be clearly defined. Such a definition should include the detection algorithm itself, the current knowledge about real dust devils and changes in the dust devil physics that follow from a change in model resolution.
Convergence of the general flow in the convective PBL agrees well with previous studies (e.g., Sullivan and Patton 2011;Bopape et al. 2020;Wurps et al. 2020). For the boundary layer interior, grid spacings on the order of 10 m were sufficient to capture the mean flow and low-order moment statistics of the convective boundary layer appropriately. However, properties of the near surface layer, where dust devils mainly occur, still vary with resolution. Additionally, previous studies never analyzed convergence of peak values, like maximum wind velocities. Although not thoroughly investigated, our results indicated that extrema do not converge. Cursory analysis showed that most of the maxima are located within dust devils at heights below 30 m. Maxima at these heights are likely related to buoyancy caused by the super-adiabatic layer that is drawn into the cores of dust devils. Due to the non-convergence of the peak values, we concentrated our statistical analysis to mean dust devil properties. These mean properties were averaged over both the whole vortex lifetime and sample at a height of 10 m. We chose this height to avoid surface interference.
We defined the dust devil population to be investigated as all vortices with lifetimes of at least 120 s and lifetime-averaged radii of at least 5 m. The properties of this population were consistent with the results of field studies if resolutions of 1.25, 1, and 0.625 m were chosen. Therefore, we generally recommend a grid spacing below 2 m if quantitative results are desired from numerical simulations. However, the quantitative comparison between different measurement campaigns and model results remains challenging (e.g., because of poor threedimensional measurement data or different analysis heights).
With a grid spacing of 1.25 or below, our results showed a converged dust devil occurrence rate of 19 km −2 h −1 or 77 km −2 d −1 . Typical lifetimes were 4-5 min. The mean converged radius was between 6 and 8 m. Typical temperature excursions were 3-4 K with maxima of approximately 8 K. The mean strength of vortices, rated by the pressure drop at convergence, was approximately 14 Pa. Maximum instantaneous pressure drops of several hundreds of pascals existed in some vortices. The vertical vorticity averaged over a horizontal plane of 20 × 20 m 2 around the dust devil center typically ranged between 0.1 and 1 s −1 . Finally, converged mean tangential and radial velocities were 2.1 m s −1 and 1.35 m s −1 , respectively, while maxima were 13 and 4 m s −1 , respectively. The vertical velocity never converged, even for the mean values. They constantly increased as grid spacing decreased from 1 to 2 m s −1 . At the finest resolution of 0.625 m, peak vertical velocities of nearly 20 m s −1 were simulated. A comparison with observations suggested that convergence is expected at resolutions of just less than 0.5 m, which might be affordable within the next years.
Finally, the analysis of the instantaneous and time-averaged three-dimensional flow structure indicates realistic results for grid spacings of 2.5 m or smaller. With such a resolution, all vortex features are qualitatively captured. This especially includes central downdrafts, sub-vortices and the near-ground convergent flow that becomes a dominant spiraling upward motion near the center. Also, the vertically-thin, super-adiabatic layer with high vertical temperature gradients in the dust devil surroundings, which is drawn into and stretched within the core, is captured appropriately. In instantaneous pressure fields, we frequently observe tail-like structures that separate strong updraughts from downdraughts. These structures have not been reported in the current literature.
Follow-up grid convergence studies should focus especially on convergence of extrema. Namely, do extrema converge at resolutions of just less than 0.5 m? These very high-resolution studies could then investigate the small-scale structures of dust devils in more detail. The results presented here are useful for deriving better and more convincing quantitative estimates of dust devil phenomena (e.g., the amount of dust that is typically released by dust devil-like vortices). It it now clear that results derived from simulations at 10 m grid spacing or larger are not meaningful because pressure drops and wind velocities, which are directly related to turbulent dust emission, are too weakly simulated. In ongoing work, we are currently estimating the sediment fluxes and particle concentrations within dust devils. Generally speaking, this study can be taken as the basis for the resolution requirements that are needed to capture coherent vortex structures like dust devils in future LES studies.