LES study on turbulent dust deposition and its dependence on atmospheric boundary-layer stability

It is increasingly recognized that atmospheric boundary-layer stability (ABLS) plays an important role in aeolian processes. While the effects of ABLS on dust emission have been documented in several studies, those on dust deposition are less well studied. By means of large-eddy simulation, we investigate 10 how ABLS influences the probability distribution of surface shear stress and hence dust deposition. Statistical analysis of the model results reveals that the shear stress can be well approximated by using a Weibull distribution and the ABLS influences on dust deposition can be estimated by considering the shear stress fluctuations. The model-simulated dust depositions are compared with the predictions of a dust-deposition scheme and measurements, and the findings are then used to improve the dust-deposition scheme. This 15 research represents a further step towards developing dust schemes that account for the stochastic nature of dust processes.

observed that when the background wind speeds are similar, deposition velocity under convective conditions is in general larger than under neutral and stable conditions. Pellerin et al. (2017) suggested that cospectral 35 similarities exist between heat and particle-deposition fluxes and that atmospheric turbulence plays a role in dust deposition. It is therefore necessary to find a link between instantaneous wind and particle deposition and to correctly represent this link in particle-deposition schemes, i.e., to introduce and account for the effect of turbulence on particle deposition. Models for turbulent dust emission ( Klose andShao, 2012, 2013) and sand saltation (Liu et al., 2018;Li et 40 al., 2020;Rana et al., 2020) have been developed, but to our knowledge, although turbulent dust deposition is now perceived to be important, a scheme is yet to be constructed for its quantitative estimate.
The turbulent wind flow in a particle-deposition scheme is reflected in the turbulent shear stress (or vertical momentum flux). It is well-known that apart from gravitational settling, particle deposition is driven by turbulent diffusion which is intimately related to the vertical momentum transfer in the atmospheric boundary 45 layer (ABL) (Wyngaard, 2010). Based on the Prandtl mixing-length theory, the shear stress can be parameterized in neutral conditions. However, it is known that for a given mean wind speed (at a reference height) in the ABL, both the mean value and the perturbations of shear stress depend on the atmospheric boundary-layer stability (ABLS), for instance, shear stress shows generally larger fluctuations in convective ABLS. Klose and Shao (2012) pointed out that, under convective conditions, large eddies have coherent 50 structures of dimensions comparable to boundary-layer depth, which are efficient entities in generating localized momentum fluxes to the surface. Although the eddies only occupy fractions of time and space, the momentum fluxes to these fractions can be many times the average. Hicks et al. (2016) mentioned that ABLS is of immediate concern in the micrometeorological community, because of its influences on the intermittency, gustiness and diurnal cycle of particle deposition. Similar to dust emission and sand saltation, 55 intermittent dust deposition also occurs as a result of fluctuating surface shear stress. The current particledeposition schemes only consider the mean behavior of wind, and how this mean behavior varies with ABLS via the Monin-Obukhov similarity theory, (Monin et al., 2007;Monin and Obukhov, 1954), but not the fluctuations of the associated shear stress and how they vary with ABLS.
We argue that focusing only on the effects of ABLS on mean wind is insufficient to model particle deposition 60 accurately. In this study, we explore the influences of ABLS on the turbulent behavior of particle deposition and attempt to improve an existing particle deposition scheme. A large-eddy simulation (LES) model is used to simulate turbulence and particle deposition under various ABLS conditions. The dust particle depositions simulated using the LES model and predicted using the particle-deposition scheme of Zhang and Shao (2014, ZS14 hereafter) are compared with each other and with measurements. Here, we address the following three 65 issues: (1) How ABLS affects the probability distribution of surface shear stress; (2) How ABLS impacts on particle deposition; and (3) How the ZS14 scheme can be improved to account for the ABLS effect. On this basis, an improvement to the ZS14 scheme (also applicable to other schemes) is proposed. The remaining part of the paper is organized as follows: Sect. 2 gives a brief description of the Weather Research and https://doi.org/10.5194/acp-2021-809 Preprint. Discussion started: 29 November 2021 c Author(s) 2021. CC BY 4.0 License.
Forecast -Large-Eddy Simulation Model with Dust module (WRF-LES/D), the ZS14 scheme, and the design 70 of the numerical experiments. Sect. 3 discusses the findings of the numerical simulations and the improvement to the ZS14 scheme. The concluding remarks are given in Sect. 4.

WRF-LES/D
The WRF-LES/D used here is initially developed by Shao et al. (2013) and Klose and Shao (2013) by 75 coupling the WRF-LES (Moeng et al., 2007;Skamarock et al., 2008) with a land-surface module and dust module. As demonstrated in the earlier studies, WRF-LES/D is a reasonably well-established system for applications to simulating turbulence, turbulent dust emission and transport for various ABLS conditions. WRF-LES is a three-dimensional and non-hydrostatic model for fully compressible flow. The model separates the turbulent flow into a grid-resolved component and a subgrid component. The k-l subgrid closure 80 (Deardorff, 1980) together with the TKE equation (Skamarock et al., 2008) based on nonlinear backscatter and anisotropic (Kosović, 1997;Mirocha et al., 2010) where ( , , ) is grid-resolved flow velocity along xi (x, y, z) refer to the streamwise, spanwise, and 90 vertical directions, respectively; g is the acceleration due to gravity; is the air density; f is the Coriolis parameter; p is pressure; is subgrid stress tensor modeled using an eddy viscosity approach where the eddy viscosity is represented as the product of a length scale and a velocity scale characterizing the subgridscale (SGS) turbulent eddies, with the velocity scale being derived from the SGS TKE and the length scale from the grid spacing; ν is the kinematic viscosity; is the Kronecker operator and is the alternating 95 operator; cp is the specific heat of air at constant pressure; T is air temperature; Hj is the jth component of subgrid heat flux; c is dust concentration; wt is dust particle terminal velocity; Fj is the jth component of subgrid dust flux; sT and sr are the source or sink terms for heat and particles, respectively. The eddy https://doi.org/10.5194/acp-2021-809 Preprint. Discussion started: 29 November 2021 c Author(s) 2021. CC BY 4.0 License. diffusivity is obtained using eddy viscosity dividing the Prandtl number. For the surface layer, an important parameterization to solve the governing equations for high-Reynolds-number turbulence is embedded in the 100 surface boundary condition, which computes the instantaneous local surface shear stress using the bulk transfer method (Kalitzin et al., 2008;Kawai and Larsson, 2012;Piomelli et al., 2002;Zheng et al., 2020) as follows, where being eddy viscosity and being the MOST stability function, = √ + . Even though Shao et al. (2013) questioned the application of the MOST in LES, it is still used here, as our emphasis is on the variance of shear stress in the simulation domain. Several land-surface models (LSMs) can be selected (e.g., Chen and Dudhia, 2000;Pleim and Xiu, 2003) in WRF-LES/D, and the 5-layer thermal diffusion 110 (Dudhia, 1996) is used in this study. Surface heat flux in this study is artificially given. In addition, we denote surface heat flux as H0 and dust dry deposition flux on grand in each grid as Fd.

Particle-deposition scheme of ZS14
The dust particle deposition on the surface is more complicated than momentum flux as the dust concentration changing close to the surface is unclear. To solve the dust conservation equation, Eq. (5), the emission and 115 deposition fluxes at the surface need to be specified. The problem of dust emission has been dealt with elsewhere (Shao, 2004;Klose and Shao, 2013) and is not considered here. For our purpose, dust emission is assumed to be zero. This section gives the parameterization scheme of surface settlement proposed by ZS14.
The detail of the scheme is as described in ZS14, only the main results are given here for completeness. In general, we can express dust deposition flux Fd as 120 where Kp and kp are eddy diffusivity and molecular diffusivity, respectively. By analogy with the bulk-transfer formulation of scalar fluxes in ABL, Fd can be parameterized as where c(z) is the dust concentration at height z (the center height of the lowest model level in this study), 125 Vd(z) is the corresponding dry deposition velocity.
The surface layer is divided into an inertial layer and a roughness layer. Integrating Eq. (8)

 
with ra being the aerodynamic resistance for the inertial layer. Using the MOST, we have 130 where zd is the displacement height, h is the height of roughness element ѱm is the integral of stability function in the inertial layer, = ⁄ (Csanady, 1963), and κ is the von Karman constant. The gravitational resistance rg is defined as the reciprocal of the gravitational settling and depends mainly on particle size and density. A free-falling particle is subject to gravitational and aerodynamic drag forces. When these forces are 135 in equilibrium, the gravitational settling velocity reaches the terminal velocity given by the Stokes formula where Dp is the particle diameter, ρp is the particle density, μa is the air dynamic viscosity, Cu is the Cunningham correction factor that accounts for the slipping effect affecting the fine particles.
In the roughness layer, the collection process is reflected in collection resistance, defined by = − ( ) with 140 an assumption of dust concentration is zero on roughness elements or ground. In addition to the meteorological factors and land-use category, Zhang and Shao (2014) established a relationship between aerodynamic and surface-collection processes by using an analogy between drag partition and deposition flux partition, which can describe surface heterogeneity.
Here, R is the reduction of collection caused by particle rebound, E is the collection coefficient of the roughness elements, which includes the contributions of Brownian motion, impaction and interception, Sc is the Schmidt number which is the ratio of air viscosity to molecular diffusion, uh is the wind speed at the top of roughness layer, Cd is the drag coefficient for isolated roughness element, 10 represents the turbulent impaction efficiency with being the dimensionless particle relaxation time. The ratio τc/τ describes the drag 150 partition with τc being the pressure drag (the force exerted on roughness elements) and can be calculated according to Yang and Shao (2006) with β (= 200) is the ratio of the pressure drag coefficient for isolated roughness element to that of bare surface, λ is the frontal area index of the roughness elements, η is the basal area index or the fraction of cover.
From Eqs. (10)-(15), it can be seen that Vd and τ are nonlinearly related. As example, for the particle of diameter 1 μm, analysis shows that when τ is small, Vd is dominated by wt. As τ increases, wt and τ are both important in Vd. As τ increases further, the effect of τ becomes much larger than gravity settling, thus the Vd 160 is mainly determined by τ.

Simulation Set-up
Numerical experiments are carried out with WRF-LES/D for various atmospheric stability and backgroundwind conditions for two different roughness lengths (Table 1) For model initialization, the wind and dust concentration (Chamberlain, 1967;Monin, 1970;Kind, 1992) are 170 assumed to be logarithmic in the vertical and uniform in the horizontal direction. For each experiment, a constant surface heat flux is specified. A 300 m deep Rayleigh damping layer is used at the upper boundary with a damping coefficient of 0.01. The wind speed at the top boundary, U, is given in Table 1    and z0 = 0.76 mm for Exp (21-35) in field observation (Bergametti et al., 2018)

Turbulent shear stress 185
In the first set of the analysis, we examine the impact of atmospheric stability on shear stress fluctuations.
Early dust deposition studies considered only the time average of surface shear stress, , with the assumption that shear stress is horizontally homogeneous. In WRF-LES/D, the corresponding mean resultant shear stress can be obtained as below: The shorthand notation 1 ( , , )  addition, the insert plots in Fig. 1 show that the autocorrelation function, ACF, is oscillated during decrease.
The oscillation periodicity is longer under weak wind conditions (Fig. 1a-c) than strong wind (Fig. 1d-f). The 200 ACF in neutral conditions decreases rapidly than in convective conditions. Recall the definition of coherent motion given by Robinson (1991) -the correlation of variables over a range of long time larger than the smallest scales of flow is an evidence of coherent oscillating motion. Thus, the regular oscillation and a long time correlation of τ are closely related to the evolvement of the coherent structure. This indicates that in a convective ABL, stronger large-scale coherent structures exist even under weak wind conditions. 205 To gain insight into the behavior of the unsteady shear stress field, we introduce the turbulence intensity of surface shear stress (TI-S) defined as the ratio of the standard deviation of fluctuating surface shear stress, , to the mean resultant stress , i.e., / . Analysis shows that / increases as atmospheric conditions become more unstable and decreases with wind speed (e.g., Fig. 1). High wind speeds tend to force the ratio to be more similar to that in neutral ABLs, as the mean-wind induced shear stress becomes dominant over 210 the large-eddy induced shear-stress fluctuations. For a weak TI-S, is dominated by and the stress fluctuations are small compared to . As TI-S increases, the contribution of momentum transport by large eddies becomes significant because in unstable ALBS, buoyance generated large eddies penetrate to high levels and intermittently enhance the momentum transfer to the surface. The intermittent surface shear stress can directly cause localized dust deposition. Therefore, dust deposition is also intermittent in space and time. However, to our knowledge, in existing dust-deposition schemes (e.g., ZS14 used here), the dust-deposition velocity is calculated using only the mean resultant shear stress 220 instead of the instantaneous shear stress. We denote this deposition velocity as , . The mean deposition As shown, the fluctuation behavior of Vd is consistent with that of . Moreover, Fig. 2a shows a substantial difference between , and , , while Fig. 2b shows , is similar with , . This suggests that the ZS14 scheme can more accurately estimate the deposition velocity for weak TI-S but underestimates the 230 deposition velocity for strong TI-S. The reason for this is that in the case of strong TI-S, dust deposition caused by the gusty wind plays an important role as Vd and are non-linearly related, which is not reflected in , . Since fluctuates and sometimes strongly, a bias always exists in conventional dust-deposition schemes and the magnitude of the bias depends on turbulence intensity. Therefore, in order to estimate dust deposition accurately, we need to first describe and parameterize the shear stress. 235

240
As a main predisposing factor for aeolian processes, turbulent shear stress has attracted much attention (e.g., Klose et al., 2014;Li et al., 2020a;Liu et al., 2018;Rana et al., 2020;Zheng et al., 2020). Similar to previous studies, we use the probability density function p(τ) to characterize the stochastic variable τ. Figure 3 shows that the variability of τ increases as atmospheric instability increases. The statistic moments of τ, including its mean resultant value , standard deviation , skewness γ1 of Exp (1-20) are listed in Table 2. and 245 increases with increased instability, and the distribution is positively skewed. Positive skewness is characterized by the distribution having a longer positive tail as compared with the negative tail and the https://doi.org/10.5194/acp-2021-809 Preprint. Discussion started: 29 November 2021 c Author(s) 2021. CC BY 4.0 License. distribution appears as a left-leaning (i.e., tends toward low values) curve. This indicates that large negative fluctuations are not as frequent as large positive fluctuations. The data also shows γ1 generally shows a downward trend as TI-S decreases, which is consistent with (Monahan, 2006), i.e., as TI-S decreases, p(τ) 250 becomes increasingly Gaussian.  The parameterization of surface shear stress has attracted intense interests, as example, Klose et al. (2014) 260 reported that τ in unstable conditions is Weibull distributed based on large-eddy simulations. Shao et al. (2020) found that p(τ) is skewed to small τ values (i.e., positively skewed) based on field observations. Li et al. (2020) suggested that τ in neutral conditions is Gauss distributed based on a wind-tunnel experiment. Colella and Keith (2003) explained that in turbulent shear flows, the non-linear interaction between the eddies gives rise to a departure from Gaussian behavior. Our results show that the Gaussian approximation is inadequate in 265 representing the skewed p(τ), especially for the conditions of strong turbulence intensity (e.g., unstable cases in Fig. 3a). Therefore, p(τ) here is approximated using a Weibull distribution, i.e., where α and β are the shape and scale parameters, respectively. The values of α and β for the numerical experiments Exp (1-20) are listed in Table 2. It can be seen that both α and β depend on wind speed and 270 atmospheric stability. However, β is mainly determined by wind conditions when wind is strong, while it is affected by ABL stability when wind is weak. The behavior of α and β are shown in Fig. 4. In both stable and unstable atmospheric conditions, analysis shows that the scale parameter α is related to ABL stability as the power of |1 ⁄ | where is the Monin-Obukhov length. Fig. 4a shows that α decreases with the |1 ⁄ |, satisfying approximately Eq. (19). For neutral conditions, Lo goes to infinity, Eq. (19) no longer applies. 275 Therefore, the shape parameter obtained by the fitting was directly used for pdf reproduction for the neutral https://doi.org/10.5194/acp-2021-809 Preprint. Discussion started: 29 November 2021 c Author(s) 2021. CC BY 4.0 License. cases instead of the approximated α used for stable and unstable conditions. As Fig. 4b shows, the β parameter increases almost linearly with * + 0.001 • * but can be best approximated using Eq. (20) with * = ⁄ .
Using Eqs. (18)-(20), we can approximately describe the turbulent surface shear stress in non-neutral cases. 285 3.2 Improvement to dust deposition scheme Figure 5a shows the performances of WRF-LES/D by comparing the simulated deposition velocity, , , with wind tunnel experiments  and field observation (Bergametti et al., 2018). The observed data are measured under neutral conditions and similar wind flow. As shown, the simulation results agree well with the observed data. On this basis, we further evaluate the performance of the ZS14 scheme, 290 and show that the accuracy of ZS14 scheme decreases as instability increases. As examples, Fig. 5b compared , , of Exp (5,9,17) and Exp (24,27,33) with the ZS14 scheme result , which is calculated using .
with p(τ) is as given by Eqs. (18)-(20). As Fig. 5c shows, the improved scheme results , and the simulation 300 value , are shown a remarkable congruence.
Analysis shows that the value of relative error, RE, depends on surface conditions, wind conditions, atmospheric stabilities, and particle sizes. It increases obviously with increased atmospheric instability under 315 https://doi.org/10.5194/acp-2021-809 Preprint. Discussion started: 29 November 2021 c Author(s) 2021. CC BY 4.0 License. weak wind conditions, while becomes less sensitive to stability when wind is strong. Through the analysis, we find that the RE of ZS14 scheme generally increases with the shear stress turbulence intensity, TI-S, and the value depends on particle size, as shown in Fig. 5d (left). Thus, we compared the RE of some different sized particles to investigate that the particle in which size range is strongly affected (Fig. A2). The result shows that RE first increases and then decreases with increasing particle size, and the particles with size 320 normally in the range of 0.01 to 10 are strongly affected by turbulent shear stress and p(τ) needs to be considered. After modification, the errors are limited to less or about 10%. For example, the relative error of Exp (17, i.e., U = 4 m s -1 and = 600 W m -2 ) for particles of 1.46 μm is reduced from ~ 25% to ~ 3%. The relative error of Exp (33, i.e., U = 5.44 m s -1 and = 600 W m -2 ) for particles of 0.5 μm is reduced from ~ 50% to ~ 12%. 325 To further analyze if the RE of ZS14 in unstable conditions is dominated by kinetic instability or dynamic instability, the Richardson number is calculated. Analysis shows that TI-S is positively correlated to gradient Richardson number Ri, and RE of ZS14 is increasing with the magnitude of Richardson number Ri under convection predominant unstable condition associating weak winds and strong vertical motion (Fig. A3). The relationship between Ri and TI-S needs further study. Consequently, the results illustrate that the modified 330 scheme , tends to be more accurate than the unmodified scheme , .

Conclusion
The present study was designed to determine the effect of ABL stability on dust particle deposition. For this purpose, the WRF-LES/D was used to model atmospheric boundary-layer turbulence under the presence of atmospheric stability effects to recover statistics of shear stress variability. We then presented an improved 335 dust-deposition scheme with the consideration of turbulent shear stress. While ABLS can broadly represent levels of atmospheric turbulence, its effect on dust deposition is wind speed dependent. Through a series of numerical experiments, we have shown the turbulent characteristics of dust deposition velocity caused by the turbulent wind flow and pointed out existing dust-deposition schemes have deficiencies in representing dust deposition under convective conditions. The relative error RE increases as the ABL instability increases for 340 low wind conditions, i.e., RE increases with shear stress turbulence intensity, especially for a certain size range of particles.
Since the dependency of dust deposition on micrometeorology is imbedded in the application of the surface shear stress, we believe that the dependency of dust deposition on ABL stability is ultimately attributed to the statistical behavior of shear stress τ. Therefore, in this study, a model including the effects of surface shear The project is the first comprehensive investigation of the turbulent characters of dust deposition and the findings will be of interest to improve the accuracy of dust-deposition predictions in regional or global scales.
One source of weakness in this study is the variation of τ may change with surface roughness and needs further study. In spite of this limitation, the study adds to our understanding of the influence caused by ABLS on particle deposition. 355 Appendix Figure A1 shows the probability density distribution of surface shear stress for experiments (21-35); Figure   A2 shows the changing of relative error with particle size; Figure A3 shows the variation of relative error (RE) of the ZS14 scheme (Eq. (10)) and improved scheme (Eq. (21)) with gradient Richardson number Ri.   Figure A3. Comparison of relative error as a function of Ri, estimated by ZS14 scheme (circles) and the improved scheme (crosses) for Exp (1-20) (left) and Exp (24, 27, 30, 33) (right).

Code and data availability 370
The source code used in this study is the WRF-chem version 3.7 in the LES mode coupled with a new deposition scheme. WRF-LES model can be download at https://www2.mmm.ucar.edu/wrf/users/download/get_sources.html. The code of the coupled deposition scheme and data set obtained by the simulation are available online at https://github.com/YinXin2021/WRF-

LES-DustDepositionScheme. 375
Author contributions XY, YPS and JZ were responsible for the formal analysis, Methodology. XY and CJ were responsible for the data curation, software, validation and visualization. YPS, JZ and NH were responsible for the supervision, project administration and funding acquisition. XY was responsible for investigation and Writing -original draft preparation. XY, YPS and JZ were responsible for the Writing -review & editing. 380

Competing interests
The authors declare that they have no conflict of interest.

Disclaimer
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.