2.1. Carrier phase

In this work, we use the NCAR Turbulence with Lagrangian Particles Model (Richter & Chamecki Reference Richter and Chamecki2018) to model one-way coupled inertial particles settling through a TBL. This code has been validated and used in multiple studies focused on inertial particle settling and transport in TBLs (Richter & Chamecki Reference Richter and Chamecki2018; Wang et al. Reference Wang, Fong, Coletti, Capecelatro and Richter2019; Bragg et al. Reference Bragg, Richter and Wang2021a ; Gao, Samtaney & Richter Reference Gao, Samtaney and Richter2023; Grace, Richter & Bragg Reference Grace, Richter and Bragg2024). For the carrier phase, we use DNS to solve the three-dimensional, incompressible Navier–Stokes equations in a turbulent open channel flow set-up,

(2.1) $\frac {D\boldsymbol{u}}{Dt} = -\frac {1}{\rho _a}\nabla p + \nu \nabla ^2\boldsymbol{u} - \frac {1}{\rho _a}\frac {{\textrm d} P}{{\textrm d}x}\hat {\boldsymbol{x}}, $

(2.2) $ \nabla \cdot \boldsymbol{u} = 0. $

A schematic of the set-up is presented in figure 1. In the above equations, ${D}/{Dt}$ represents the material derivative, $\boldsymbol{u}$ represents the three-dimensional flow velocity, $p$ represents the turbulent pressure field, $\rho _a$ is the carrier phase density and $\nu$ is the kinematic viscosity. Unconditional fluid velocities will be referred to by their components ( $u,v,w$ ) with no subscript, and fluctuating quantities will be denoted by a prime. For example, as we are primarily focused on wall normal motion in this work, we will refer to fluctuating vertical fluid velocities as $w^\prime$ . The code uses a pseudospectral method in the horizontal directions and a second-order finite difference method in the vertical. We have clustered points near the solid lower boundary using an algebraic grid stretching procedure, so that $\Delta z_1 =0.8\nu /u_\tau$ ( $u_\tau$ is the friction velocity, defined below), where $\Delta z_1$ is the first grid point and $(\Delta x_\eta ,\Delta y_\eta ,\Delta z_\eta ) = (2.8,1.4,0.9)$ at the midplane of the flow, where the subscript $\eta$ represents the grid spacing in terms of the local Kolmogorov scale, defined below.

Figure 1. A schematic of the numerical set-up. The domain is a rectangular channel of height $H$ , streamwise length $2\pi\!{H}$ and spanwise width $\pi\!{H}$ . The flow is periodic in the horizontal and is driven by a constant pressure gradient in the streamwise direction, while the no-stress and no-slip boundary conditions are enforced at the top and bottom boundaries, respectively. Particles are injected at the upper boundary at a random horizontal location with an initial velocity equal to the fluid velocity at their location and removed when they contact the bottom boundary. They are allowed to rebound elastically off the upper boundary.

At the lower boundary, a no-slip boundary condition is enforced, while at the upper boundary, a no-stress boundary condition ( $\partial u/\partial z = \partial v/\partial z = 0$ at $z = H$ ) is enforced. The domain is periodic in the $x$ and $y$ directions. The background state of the carrier phase is established by accelerating the flow with an imposed pressure gradient, $-{\textrm d}P/{\textrm d}x\gt 0$ (note that $\hat {\boldsymbol{x}}$ is the unit vector in the streamwise direction) and allowing the flow to become turbulent. The magnitude of the pressure gradient allows us to define a friction velocity $u_\tau = \sqrt {\tau _w/\rho _a}$ , where $\tau _w$ is the stress at the lower boundary. Using the friction velocity, the height of the domain, $H$ , and viscosity of the carrier phase, we can define a friction Reynolds number of $\textit{Re}_\tau = {u_\tau H}/{\nu }$ . Friction Reynolds numbers for each simulation presented in this work can be found in table 1. This set-up is identical to that used in Grace et al. (Reference Grace, Richter and Bragg2024).

Table 1. Table of cases discussed throughout the work. Parameter definitions can be found in the main text. The case with $\textit{Re}_\tau = 1260$ was run on a $512^3$ grid, all cases with $\textit{Re}_\tau =630$ were run on a $256^3$ grid, while cases with $\textit{Re}_\tau =315$ were run on a $128^3$ grid.

We can define the local Kolmogorov time scale, velocity scale and acceleration scale,

(2.3) \begin{equation} \tau _\eta = \left (\frac {\nu }{\epsilon }\right )^{1/2},\quad v_\eta = \left (\nu \epsilon \right )^{1/4},\quad a_\eta = \left (\frac {\epsilon ^3}{\nu }\right )^{1/4}, \end{equation}

respectively, which represent the smallest relevant length scales of the turbulence. These parameters will be used to characterize the turbulent flow below. In statistically stationary HIT, the above scales are constants, but for wall-bounded turbulence they depend on height. Since the turbulence intensity decreases with height outside of the very thin viscous sublayer adjacent to the wall, so too does the kinetic energy dissipation rate resulting in a height variation of the Kolmogorov microscales. When the TBL is horizontally homogeneous, the mean dissipation rate $\epsilon$ is a function of the distance from the solid boundary. Within the logarithmic layer, it scales as

(2.4) \begin{equation} \epsilon \sim O\left (\frac {u_\tau ^3}{\kappa z}\right ), \end{equation}

where $u_\tau$ is the friction velocity, and $\kappa$ is the von Kármán constant. However, when calculating the Kolmogorov time scale, velocity scale and length scale, we use the dissipation computed in the DNS.

2.2. Dispersed phase

Our main focus is towards on-coarse dust transport in the atmospheric surface layer. Dust particles can range in size, but even coarse grains (roughly 30–100 μm) are significantly smaller than the local Kolmogorov scale, which can be in the range of several millimetres. Indeed, these particles are also significantly denser than the carrier phase, though their volume fractions can be quite low once they are above the emission layer. With these assumptions in mind, for each particle (the dispersed phase), we apply the point-particle approximation and apply the conservation of momentum for a rigid spherical particle subjected to linear hydrodynamic drag and gravity. Note that we have purposefully omitted particle lift as it has been shown to play negligible role in the dynamics discussed in the present work (Costa, Brandt & Picano Reference Costa, Brandt and Picano2020). Furthermore, as we are concerned with the dilute limit, two-way coupling and particle–particle interactions may be ignored. Since these particles are assumed to be much denser than the carrier phase, we may also ignore added mass and Basset history forces. The one-way coupled point particle approximation also has the added benefit that each particle is independent from each other particle, effectively removing the volume fraction as a governing parameter as the number of particles is increased. This allows us to increase the number of particles to ensure convergence of the statistics of interest without affecting the flow. As such, our results are assumed to be valid in the dilute limit, where particle densities are much greater than that of the fluid phase.

The equations of motion are

(2.5) \begin{equation} \frac {{\textrm d}\boldsymbol{u}_p}{{\textrm d}t} = \frac {\Psi }{\tau _p}(\boldsymbol{u}_f - \boldsymbol{u}_p ) - \boldsymbol{g}, \end{equation}

(2.6) \begin{equation} \frac {{\textrm d}\boldsymbol{x}_p}{{\textrm d}t} = \boldsymbol{u}_p. \end{equation}

Here, $\boldsymbol{u}_p = (u_p,v_p,w_p)$ is the three-dimensional velocity vector for each particle, $\boldsymbol{x}_p = (x,y,z)$ is the location of each particle in space, $\boldsymbol{g}$ is the gravitational acceleration (which only affects accelerations in the $z$ direction), $\boldsymbol{u}_f = (u_f,v_f,w_f)$ is the three-dimensional instantaneous flow velocity evaluated at the location of the particle. Much of our focus will be placed on the vertical component which we denote as $w_f$ (not bold).

Particles are injected at the upper boundary at a random horizontal location with an initial velocity equal to the fluid velocity at their location. We have performed a sensitivity test to investigate the importance of the particle initial condition by comparing results when injected at their laminar Stokes settling velocity. We have concluded that the choice of initial particle velocities at injection makes a negligible quantitative change to our results, and does not affect the results qualitatively.

Particles are removed when they contact the bottom boundary, and once removed, are reinjected at a random horizontal location at the upper boundary. This maintains a constant particle number in the domain, and a constant vertical particle flux. Particles are allowed to come to statistical equilibrium, and statistics are recorded for no less than 100 eddy turnover times after this time. Finally, particles are allowed to rebound elastically off the upper boundary, and the particle boundary conditions are periodic in the horizontal directions.

The relaxation time scale of the particles, $\tau _p$ , is defined as

(2.7) \begin{equation} \tau _p = \frac {\rho _pd_p^2}{18\rho _a\nu }, \end{equation}

where $\rho _p$ is the particle density, $d_p$ is the particle diameter and the Stokes settling velocity is defined as $v_g = \tau _p g$ . Here $\Psi = 1 + 0.15\textit{Re}_p^{0.687}$ is the Schiller–Neumann correction to the drag force, and $\textit{Re}_p$ is the particle Reynolds number. The particle Reynolds number remains small enough such that $\Psi \approx 1$ . For the theory discussed below in § 2.3, we follow Bragg et al. (Reference Bragg, Richter and Wang2021a ) and make the assumption that $\Psi = 1$ for analytical tractability.

We can now define a set of non-dimensional parameters characterizing the system,

(2.8) \begin{equation} {S{\kern-0.5pt}t}^+ = \frac {\tau _pu_\tau ^2}{\nu },\quad {S{\kern-0.5pt}v}^+ = \frac {v_g}{u_\tau },\quad {St}_\eta = \frac {\tau _p}{\tau _\eta },\quad {S{\kern-0.5pt}v}_\eta = \frac {v_g}{v_\eta }. \end{equation}

These are the friction Stokes number and the settling velocity parameter based on the viscous scales, and the Stokes number and the settling velocity parameter based on the local Kolmogorov scales. Here ${St}_\eta$ and ${S{\kern-0.5pt}v}_\eta$ are functions of the distance from the solid boundary, whereas ${S{\kern-0.5pt}t}^+$ and ${S{\kern-0.5pt}v}^+$ are constant parameters.

Other studies focused on particle settling define other parameters such as a Froude number, $\textit{Fr} = a_\eta /g$ , (Bec et al. Reference Bec, Homann and Ray2014; Berk & Coletti Reference Berk and Coletti2021), or a scaled gravity

(2.9) \begin{equation} g^+ = \frac {g\nu }{u_\tau ^3}, \end{equation}

which is simply the ratio of ${S{\kern-0.5pt}v}^+$ to ${S{\kern-0.5pt}t}^+$ . These parameters are useful as they describe the relative role of gravitational accelerations and turbulent accelerations without referring to $\tau _p$ . We will refer to $g^+$ at various points throughout the discussion when necessary. The values for all parameters considered in this work, including the associated ranges of ${St}_\eta$ and ${S{\kern-0.5pt}v}_\eta$ , can be found in table 1.

2.3. Statistics of inertial particles in a turbulent boundary layer

We adopt the particle phase space approach used in Bragg et al. (Reference Bragg, Richter and Wang2021a ), focusing only on the vertical component of the particle equations of motion. First they define the particle probability density function (PDF) in position-velocity space as

(2.10) \begin{equation} {\mathcal{P}} = \langle \delta (z_p-z)\delta (w_p-w)\rangle , \end{equation}

which describes the distribution of the vertical components of the particle position and velocity, $z_p(t)$ and $w_p(t)$ , in the phase space with coordinates $z$ , $w$ , and $\langle \cdot \rangle$ represents an ensemble average over all realizations of the system (note that $\delta$ represents the Dirac delta function). Note we make frequent use of conditional averages throughout this work, denoted by $\langle \cdot \rangle _{z,w}$ which is short hand for $\langle \cdot | z_p = z,w_p = w\rangle$ . We can form an evolution equation for the PDF,

(2.11) \begin{equation} \frac {\partial {\mathcal{P}}}{\partial t} + \frac {\partial }{\partial z}\left (w{\mathcal{P}}\right ) + \frac {\partial }{\partial w}\left (\langle a_p\rangle _{z,w}{\mathcal{P}}\right ) = 0 \end{equation}

where we have defined $\langle a_p\rangle _{z,w} = \tau _p^{-1} (\langle w_f\rangle _{z,w}-w ) - g$ as the vertical particle acceleration conditioned on $z_p = z$ and $w_p = w$ based on the vertical component of (2.5). The utility of this equation comes from the fact that we can derive evolution equations for each moment. Recall that the $n$ th moment is defined as

(2.12) \begin{equation} \langle w_p^n\rangle _z = \frac {1}{\varrho }\int _{-\infty }^{\infty }w^n{\mathcal{P}}\,{\textrm d}w, \end{equation}

where the notation $\langle \cdot \rangle _z$ represents an ensemble average conditioned on $z_p = z$ and $\varrho$ is zeroth moment, which is related to the particle number concentration. Of importance to the present work will be the first and second moments (i.e. $n=1$ and $n=2$ ). The evolution equation for the first moment is

(2.13) \begin{equation} \langle w_p\rangle _z = \langle w_f\rangle _z - v_g -\frac {\tau _p}{\varrho }\frac {{\textrm d}}{{\textrm d}z}\varrho \langle w_p^2\rangle _z. \end{equation}

The details of the derivation of this equation can be found in Bragg et al. (Reference Bragg, Richter and Wang2021a ), and the general form of this equation (for the case where $v_g=0$ ) for arbitrary moments can be found in Johnson et al. (Reference Johnson, Bassenne and Moin2020). Equation (2.13) says that the average settling velocity comes from the average fluid velocity sampled by the particles, the laminar Stokes settling velocity, and a term that depends on the derivative of the mean square particle velocity which may be important near a solid boundary. For compactness, we can rearrange (2.13) into a relationship describing how the mean slip velocity varies with height,

(2.14) \begin{equation} \langle w_s\rangle _z = v_g + \frac {\tau _p}{\varrho }\frac {{\textrm d}}{{\textrm d}z}\varrho \langle w_p^2\rangle _z, \end{equation}

where $\langle w_s\rangle _z=\langle w_f\rangle _z - \langle w_p\rangle _z$ .

Our goal is to derive a continuum equation for the slip velocity variance. To do this, we multiply (2.11) by $w^2$ and integrate over all $w$ . The full details of this operation can be found in Johnson et al. (Reference Johnson, Bassenne and Moin2020). After integrating, we are left with the following relationship:

(2.15) \begin{equation} \frac {{\textrm d}}{{\textrm d}z}\varrho \langle w_p^3\rangle _z - 2\varrho \langle a_pw_p\rangle _z = 0. \end{equation}

To expand the acceleration–velocity covariance, we expand to get $\langle a_pw_p\rangle _z = \langle w_fw_p\rangle _z - \langle w_p^2\rangle _z$ and use the fact that $\langle w_s^2\rangle _z =\langle w_f^2\rangle _z - 2\langle w_fw_p\rangle _z + \langle w_p^2\rangle _z$ to arrive at

(2.16) \begin{equation} 2\langle a_pw_p\rangle _z = \frac {1}{\tau _p}(\langle w_f^2\rangle _z - \langle w_s^2\rangle _z - \langle w_p^2\rangle _z) - 2\langle w_p\rangle _zg. \end{equation}

Putting (2.15) and (2.16) together, we get an equation for the mean squared slip velocity,

(2.17) \begin{equation} \langle w_s^2\rangle _z = \langle w_f^2\rangle _z - \langle w_p^2\rangle _z -2\langle w_p\rangle _z v_g-\frac {\tau _p}{\varrho }\frac {{\textrm d}}{{\textrm d}z}\varrho \langle w_p^3\rangle _z. \end{equation}

Equation (2.17) indicates that the mean squared slip velocity variance is a function of the mean squared sampled velocity $\langle w_f^2\rangle _z$ , the mean squared particle velocity $\langle w_p^2\rangle _z$ , a drift due to the non-zero average vertical velocity $\langle w_p\rangle _z$ and the vertical derivative of the mean cubed particle velocity $\langle w_p^3\rangle _z$ , all of which are implicit functions of particle inertia and gravity.

At this point, an important distinction must be made. Though the above equation is valid for inertial particles settling through a TBL, it should be noted that when particles settle under gravity, the mean squared quantities are not equal to their variances in general. This arises from the non-zero average settling velocity. Thus, to derive a relationship between the variances, we must do a Reynolds decomposition of each term. The details of this process are omitted here, but can be found in Appendix A . The final relationship is

(2.18) \begin{equation} \langle {w^{\prime}_s}^2\rangle _z = \underbrace {\overbrace {\underbrace {\langle {w^{\prime}_f}^2\rangle _z - \langle {w^{\prime}_p}^2\rangle _z}_{(1)} + R_t}^{(2)} + R_g}_{(3)} \end{equation}

where $R_t$ and $R_g$ are defined as

(2.19) \begin{eqnarray} R_t &=& -\frac {\tau _p}{\varrho }\frac {{\textrm d}}{{\textrm d}z}\varrho \langle {w^{\prime}_p}^3\rangle _z ,\end{eqnarray}

(2.20) \begin{eqnarray} R_g &=& -\frac {\tau _p}{\varrho }\frac {{\textrm d}}{{\textrm d}z}\big(\varrho \langle w_p\rangle ^3_z + 3\varrho \langle w_p\rangle _z\langle {w^{\prime}_p}^2\rangle _z\big) + 2\langle w_p\rangle _z( \langle w_s\rangle _z -v_g), \end{eqnarray}

respectively. This model contains many terms and is quite complex, but it can be broken down into three parts, arranged in order of increasing problem complexity, and the level of complexity highlighted by the over and underbraces.

First, grouped under $(1)$ , are the terms that would appear for particles settling through homogeneous turbulence. In this limit, the slip velocity variance is determined exclusively by the difference between $\langle {w^{\prime}_f}^2\rangle _z$ and $\langle {w^{\prime}_p}^2\rangle _z$ . This is the case for particles both with and without gravity, since gravity and inertia implicitly modify these terms. Next, grouped under $(2)$ , are the terms that appear for particles dispersing vertically through a TBL in the absence of gravity. Note that all terms encompassed by $(1)$ are included in $(2)$ , but when considering those terms covered under $(1)$ in the context of a TBL, they gain implicit height dependence since their magnitudes vary with the distance from the boundary. Furthermore, at this level, a new term appears, denoted by $R_t$ . This term is proportional to the derivative of the product of the concentration and the particle velocity triple moment and increases with particle inertia. As ${Sv^+}\rightarrow 0$ , $\langle {w^{\prime}_p}^2\rangle _z$ approaches $\langle w_p^2\rangle _z$ , but $R_t$ remains, regardless. Finally, by incorporating gravity, the mean particle velocity is no longer zero, leading to a new term grouped under $(3)$ , denoted by $R_g$ . The quantities composing $R_g$ are explicitly dependent on both the inhomogeneity of the flow through the vertical derivative, and the non-zero particle settling velocity. For clarity of interpretation, the second term on the right-hand side of (2.20) is written in terms of $\langle w_s\rangle _z - v_g$ , which we can see from (2.14) is identical to ${\tau _p}/{\varrho }{\rm d}/{{\textrm d}z}\varrho \langle w_p^2\rangle _z$ . In summary, by considering the continuum equations for the first and second moment of the particle velocity, we have been able to derive an equation for the particle slip velocity. We have identified a hierarchy of terms that appear in homogeneous turbulence and TBLs with and without settling (grouped under $(1)$ ), those that appear in a TBL without settling (grouped under $(2)$ ) and those that appear in a TBL with settling (grouped under $(3)$ ). In the following section, we will identify the importance of these terms throughout the TBL.

4.1. Summary

Motivated by coarse particle transport in the atmospheric surface layer, we used coupled Eulerian–Lagrangian simulations to simulate the dynamics of ensembles of inertial particles in boundary layer turbulence. We examined the impact of particle inertia and settling on the mean and fluctuating particle slip velocity. We adapted a mathematical model discussed in Bragg et al. (Reference Bragg, Richter and Wang2021a ) and Johnson et al. (Reference Johnson, Bassenne and Moin2020) for the slip velocity variance for settling inertial particles in a TBL and highlighted the controlling factors throughout the domain. We showed that to leading order, the slip variance above of the viscous sublayer was determined by the difference between the seen variance, $\langle {w_f^\prime }^2\rangle _z$ , and the particle velocity variance $\langle {w_p^\prime }^2\rangle _z$ , except for the largest value of ${S{\kern-0.5pt}v}^+$ , where all terms in (2.18) became comparable. Consequently, as changes in the seen variance were relatively small, changes in the slip variance within the logarithmic layer were primarily governed by a decrease of the particle velocity variance, which was implicitly a function of particle inertia and the particle settling velocity (more on this below). Within the viscous sublayer, the balance became more complicated. In all cases, $\langle {w_f^\prime }^2\rangle _z$ tended towards zero to adhere to the no-slip condition enforced at the bottom boundary. However, the slip variance remained finite as the particles tended towards $z^+ = 0$ as $\langle {w_p^\prime }^2\rangle _z$ , $R_t$ and $R_g$ remained finite. The higher-order terms tended to peak within this layer, and the magnitude of the peak tended to increase with ${S{\kern-0.5pt}v}^+$ . However, by using domain averages, we demonstrated that the relative magnitude of the higher-moment terms tended to decrease as the Reynolds number increased, reflecting the fact that at higher Reynolds number, the viscous sublayer becomes much thinner, leading to a smaller contribution when averaged across the domain. As discussed above, these terms may still be important within the viscous sublayer depending on particle parameters, though.

We also showed that the fluid velocity variance along the particle trajectories exhibited only small changes with ${S{\kern-0.5pt}t}^+$ , ${S{\kern-0.5pt}v}^+$ and $\textit{Re}_\tau$ , and was approximately 80 % of the unconditional fluid velocity variance when averaged across the domain. The differences between the seen and unconditional variances are largest at the smallest ${S{\kern-0.5pt}v}^+$ considered in this work, though the differences are still relatively small. These differences are likely a result of the relatively large spectrum of turbulent motions at high Reynolds number, and the impact of crossing trajectories, which works to implicitly affect the seen variance (see the discussion in Csanady (Reference Csanady1963), for example). This conclusion is quantitatively consistent to the conclusions presented in Berk & Coletti (Reference Berk and Coletti2021) for their laboratory experiments in homogeneous turbulence. This correspondence is significant as the seen variance is not known a priori. There exist approaches to modelling this quantity (Pozorski & Minier Reference Pozorski and Minier1998; Minier & Peirano Reference Minier and Peirano2001), but the results of our work suggest that even in a TBL where there is spatial dependence in the variance of the turbulent quantities, we can approximate the fluid variance along the particle trajectories with the unconditional fluid variance without introducing significant errors. While this reduces the overall accuracy of the slip variance estimate, it increases the predictive power at the field scale, as models for the unconditional variance, such as those discussed in Kunkel & Marusic (Reference Kunkel and Marusic2006), can be employed.

To examine the relative importance of the fluctuating and mean slip, we considered the ratio of the slip variance to the total mean squared slip velocity, $\langle {w_s}^2\rangle _z$ (i.e. the square of the mean plus the fluctuation). We found that relative to the mean, the vertical slip variance decreased much faster than the horizontal as ${S{\kern-0.5pt}v}^+$ was varied. However, for both components, we showed that the overall magnitudes in the average sense did not change significantly, indicating that at relatively large ${S{\kern-0.5pt}v}^+,$ the mean slip was the determining factor in the vertical, while the fluctuating slip was the determining factor in the horizontal. This effect was also accentuated for the smallest particles considered, due to their relatively small slip variance. To further complicate the behaviour, the relative size of the slip variance tended to decrease as the Reynolds number was increased, though due to computational restrictions, we can only provide limited guidance on this issue.

We also compared the slip variance computed by the DNS with a model derived for HIT by Berk & Coletti (Reference Berk and Coletti2021) (as no such model currently exists for a TBL and is the focus of future work). The main conclusion is that the globally averaged differences (in the absolute sense) were relatively small, but the higher moment terms act as a confounding factor to reduce differences between the model and the DNS data. Thus, care should be taken when extrapolating results from low Reynolds number DNS in TBLs to higher Reynolds number experiments in homogeneous turbulence. However, as we know the size of higher moment terms tends to decrease as Reynolds number increase when integrated across the domain, we expect that DNS at higher Reynolds number should tend towards the results derived in Berk & Coletti (Reference Berk and Coletti2021) outside of the thin viscous sublayer.

Additionally, due to the inhomogeneous nature of the turbulence, non-local effects implicit to $\langle {w_f^\prime }^2\rangle _z$ and $\langle {w_p^\prime }^2\rangle _z$ may occur at large ${S{\kern-0.5pt}v}^+$ and ${S{\kern-0.5pt}t}^+$ . Isolating the importance of non-local effects in a TBL is the focus on ongoing research (for a recent example for a model of $\langle {w_p^\prime }^2\rangle _z$ in a TBL, see Zhang, Bragg & Wang (Reference Zhang, Bragg and Wang2023) and references therein), but incorporating them in a model is beyond the scope of the current article. These effects may arise due to the fact that the statistics of the turbulence may change significantly as the particle travels vertically. For example, consider the distance a settling particle travels over one relaxation time: $\delta \sim |\tau _p\langle w_p\rangle _z|$ , where $\langle w_p\rangle _z$ is the average particle settling velocity conditioned on a height $z$ given by (2.13). In order for the particle trajectory to be altered, turbulent fluctuations must be correlated over this distance. However, if this distance is comparable to the distance over which the characteristics of the turbulence change, then we expect that the particle feels the inhomogeneous nature of the flow. To formalize this quantitatively, consider the local turbulent kinetic energy at a height $z$ , $k$ . Taylor expanding about this point and truncating after the second term, we have

(4.1) \begin{equation} k(z-\delta ) = k(z)- \left .\delta \frac {{\textrm d}k}{{\textrm d}z}\right |_z. \end{equation}

To make a locally homogeneous approximation about the turbulence, we must have that

(4.2) \begin{equation} k(z)\gg \left .\delta \frac {{\textrm d}k}{{\textrm d}z}\right |_z, \end{equation}

i.e. the kinetic energy in a small neighbourhood about $z$ (defined by the distance $\delta$ ) is primarily defined by the kinetic energy measured at a height $z$ . Using (4.1), we can estimate under what conditions a locally homogeneous approximation would be appropriate by looking for cases where the second term is small compared with the first. By assuming that the gradient of the turbulent kinetic energy scales as $u_\tau ^2z^{-1}$ (also implying $k$ can be scaled by $u_\tau ^2$ ) (Smits et al. Reference Smits, McKeon and Marusic2011) we have that $z\gg \delta$ . By normalizing both sides of this inequality by the r.m.s. turbulent velocity, $u^\prime = k^{1/2}$ , we can write $\delta$ in terms of the sum of the settling enhancement, $E = ({\langle w_p\rangle _z + v_g})/{u^\prime }$ and $v_g$ , (Good et al. Reference Good, Ireland, Bewley, Bodenschatz, Collins and Warhaft2014; Loth Reference Loth2023) as

(4.3) \begin{equation} \tau _p\left |E + \frac {v_g}{w^\prime }\right | \ll \frac {z}{u^\prime }. \end{equation}

Now normalizing by $\tau _{\eta }$ , and observing that $u^\prime \sim u_\tau$ , we can simplify both side of this inequality to reveal that

(4.4) \begin{eqnarray} {St}_\eta \left |E + {S{\kern-0.5pt}v}_\ell \right | \ll \left (\frac {zu_\tau }{\nu }\right )^{1/2}, \end{eqnarray}

where we have used the dissipation scaling in (2.4) to relate $\tau _\eta$ to the vertical coordinate, $z$ . This relationship indicates that both particle inertia and gravity have an explicit role, and an implicit role (through $E$ ) to play in potential non-local effects.

We know from Good et al. (Reference Good, Ireland, Bewley, Bodenschatz, Collins and Warhaft2014) and Loth (Reference Loth2023) that in small scale laboratory experiments, $E\sim 0.2$ as ${St}_\eta$ and ${S{\kern-0.5pt}v}_\ell$ approach unity, but as both of these parameters increase, $E$ tends back towards zero, and may even become negative (Ferran et al. Reference Ferran, Machicoane, Aliseda and Obligado2023). Therefore, for the coarse particles we are concerned with in this work, we can make a locally homogeneous approximation when ${St_\eta Sv_\ell }\ll ({zu_\tau }/{\nu } )^{1/2}$ . This may not particularly restrictive for the atmospheric surface layer as the Reynolds numbers are ${O}(10^6)$ , but for laboratory experiments, the integral scales tend to scale with the size of the experimental domain (i.e. $z\sim h$ where $h$ could be the half-height of a channel). This could present a problem making a locally homogeneous approximation for coarse particles in a wind tunnel set-up.

For the DNS presented in this work, the above relationship shows that non-local effects are likely only important for cases when both ${St^+}$ and ${S{\kern-0.5pt}v}^+$ are large, as ${S{\kern-0.5pt}v}_\ell$ tends to scale with ${S{\kern-0.5pt}v}^+$ since $u^\prime \sim u_\tau$ . For example, if we consider cases with ${St^+}=100$ and ${S{\kern-0.5pt}v}^+ = 2.5$ (and assume $E\approx 0$ ), we can see immediately that (4.4) is not satisfied. This may explain the differences between the DNS and the model in Berk & Coletti (Reference Berk and Coletti2021) in figures 6(a) and 6(b) for these particles. One of the main conclusions of this work is that for the governing continuum equation for the particle slip velocity in a TBL, there are tendencies that arise due to the inhomogeneities in the turbulence associated with the presence of the wall and the fact that the particle settling velocity is non-zero. However, as we have shown, for moderately sized particles (characterized by ${S{\kern-0.5pt}t}^+$ or ${St}_\eta$ ), and ${S{\kern-0.5pt}v}^+\lt 1$ , these terms are subleading outside of the viscous sublayer. Moreover, the magnitudes of these terms in the logarithmic layer tend to diminish as $\textit{Re}_\tau$ increases. Thus, outside of the viscous sublayer, and at moderate local ${St}_\eta$ and ${S{\kern-0.5pt}v}^+$ , we can extend models designed for homogeneous turbulence (like that described in Berk & Coletti (Reference Berk and Coletti2021)) to a TBL, where we must interpret the model as local to a height $z$ . However, outside of this regime (i.e. for very large and strongly settling particles), there are implicit non-local effects that appear as particles tend to settle through the flow due to the vertical variation of the turbulent statistics along the particle trajectories.

4.2. Implications for modelling coarse particle transport in the atmospheric surface layer

Interpreting DNS results in terms of the laboratory or field scales must be done with care, as the Reynolds numbers in DNS numbers are much smaller than those found at these scales. However, by scaling up the results in this work, we can gain valuable qualitative insights into the drag on inertial settling dust particles. For example, using estimates of turbulent dissipation for an atmospheric surface layer of roughly $10^{-3}$ $\mathrm{m^2\, s^{-3}}$ , we can define a rough Kolmogorov time scale as $10^{-1}$ s. Thus, for quartz dust particles ( $\rho _p = 2650$ $\mathrm{kg\,m^{-3}}$ ) that range between 30–100 μm, we should expect a value of ${St}_\eta$ to range between 0.1–10. We can see that the ranges in our DNS are in the correct neighbourhood to model these same coarse dust particles. Moreover, we can use the values of $g^+$ from table 1 (recall $g^+ = {Sv^+/St^+}$ ) to estimate an equivalent friction velocity, $u_\tau ^*$ (note that since we are rescaling, $u_\tau ^*$ is necessarily different than the value of $u_\tau$ used in this work), which effectively gives us a qualitative estimate of the intensity of the turbulence in an atmospheric surface layer. The effective friction velocity is given by

(4.5) \begin{equation} u_\tau ^* = \left (\frac {g\nu }{g^+}\right )^{1/3}. \end{equation}

Assuming $g = 9.81\mathrm{\,m\,s^{-2}}$ and $\nu = 1.57 \times 10^{-5}\mathrm{\,m^2\, s^{-1}}$ , and some relationship between 10 m wind velocity and the friction velocity (see Kantha & Clayson (Reference Kantha and Clayson2000) for example), this gives us a proxy for wind speed at the field scale. For the values of $g^+$ in this manuscript (see table 1), the effective friction velocities vary between 0.09 $\mathrm{m\,s^{-1}}$ and 0.84 $\mathrm{m\,s^{-1}}$ , which covers a wide range of friction velocities on the Earth (Vickers, Mahrt & Andreas Reference Vickers, Mahrt and Andreas2015). Since $g^+$ is proportional to ${S{\kern-0.5pt}v}^+$ , we can see the effective wind speed increases as ${S{\kern-0.5pt}v}^+$ decreases.

Therefore, the insight we can gain is that the slip velocity in high wind conditions (small ${S{\kern-0.5pt}v}^+)$ should be primarily governed by the fluctuations associated with the turbulence, as opposed to the mean induced by gravitational settling and the presence of the solid boundary. Conversely, at lower wind speeds, the drag induced by turbulent fluctuations is much smaller relative to the mean slip. Thus, the magnitude of the slip velocity should instead be controlled by the average, which itself is controlled by the Stokes settling velocity and the turbophoretic term. Likewise, we expect the higher moment terms governing the slip variance (i.e. $R_t$ and $R_g$ ) to be more important to the dynamics farther away from the surface in this limit, relatively speaking.

This is significant when applying models like BC2021 to particle transport in field scale systems. For example, as we have described previously, BC2021 can be used (in conjunction with a model for the unconditional fluid velocity variance) to predict the slip velocity variance for inertial settling particles in homogeneous turbulence. Under low ${S{\kern-0.5pt}v}^+$ conditions, our results show that the magnitude of the slip velocity is primarily governed by its fluctuating component, which is in turn associated with the interactions with the turbulence. Moreover, as we have discussed, a locally homogeneous approximation may be used when ${S{\kern-0.5pt}v}^+$ is small enough (see (4.4)), our work suggests that BC2021 can also be applied to inhomogeneous turbulence, like that of the atmospheric surface layer, provided we are not concerned with dynamics too close to the ground and the wind conditions are strong enough. Another interesting related application is towards modelling the particle Reynolds number, which is known to affect the associated drag on the particles (Balachandar Reference Balachandar2009; Berk & Coletti Reference Berk and Coletti2024). For example, it is known that loitering effects are typically associated with large particle Reynolds number (Rosa et al. Reference Rosa, Parishani, Ayala and Wang2016), and these loitering effects work to reduce the average particle settling velocity (Good et al. Reference Good, Ireland, Bewley, Bodenschatz, Collins and Warhaft2014). Accurate modelling of loitering effects could explain discrepancies between numerical simulations and laboratory experiments with respect to the measurement of settling velocities (Ferran et al. Reference Ferran, Machicoane, Aliseda and Obligado2023). Moreover, our results may gain some insights into further than expected horizontal transport of giant dust particles off of the West African coast (Van Der Does et al. Reference van der Does, Knippertz, Zschenderlein, Giles Harrison and Stuut2018), which could be linked to loitering effects.