Multi-scale interactions in turbulent mixed convection drive efficient transport of Lagrangian particles

Abstract

When turbulent convection interacts with a turbulent shear flow, the cores of convective cells become aligned with the mean current, and these cells (which span the height of the domain) may interact with motions closer to the solid boundary. In this work, we use coupled Eulerian–Lagrangian direct numerical simulations of a turbulent channel flow to demonstrate that, under conditions of turbulent mixed convection, interactions between motions associated with ejections and low-speed streaks near the solid boundary and coherent superstructures in the interior of the flow interact and lead to significant vertical transport of strongly settling Lagrangian particles. We show that the primary suspension mechanism is associated with strong ejection events (canonical low-speed streaks and hairpin vortices characterised by $u'\lt 0$ and $w'\gt 0$, where $u'$ and $w'$ are the streamwise and vertical turbulent velocity fluctuations), whereas secondary suspension is strongly associated with large-scale plume structures aligned with the mean shear (characterised by $w'\gt 0$ and $\theta '\gt 0$, where $\theta$ represents temperature fluctuations). This coupling, which is absent in the limiting cases (pure channel flow or free convection) is shown to lead to a sudden increase in the interior concentration profiles as ${Ri}_\tau$, the friction Richardson number, increases, resulting in concentrations that are larger by roughly an order of magnitude at the channel midplane.

1. Introduction

It is well known through numerical simulation, experiment and theory that fluid turbulence significantly influences the transport of heavy Lagrangian particles with appreciable inertia. For example, turbulence generated near the Earth’s surface can lead to long-range transport ( ${\gt } 2000$ km) of giant dust grains ( ${\gt}100\,\unicode{x03BC}$ m) (Adebiyi et al. 2023; Van Der Does et al. 2018), while simple scaling theories predict a much shorter displacement due to their strong settling behaviour. Understanding and predicting the ultimate fate of heavy particles is particularly interesting to a wide range of scientific disciplines due to the climatic and health impacts of atmospheric aerosols and particles, and accurate modelling of the primary mechanisms responsible for this transport are of key interest. In this ‘mysterious long-range transport’ problem, Van Der Does et al. (2018) hypothesised several mechanisms which may influence particle transport including strong convection and high wind speeds, on which we place our primary focus. Interestingly, under convective conditions and moderate to strong mean shear, coherent, roll-like structures arise, and are well documented in the Earth’s atmosphere (Kuo 1963; Salesky et al. 2017 and references therein). In idealised studies (discussed more below), this phenomenon is often referred to as ‘mixed convection’. In this work, our goal is to investigate the consequences of mixed convection on the suspension and transport of strongly settling inertial particles.

When turbulence production by both convection and shear are comparable, the flow is termed ‘mixed convection’ (the limiting cases being a pure channel flow and free convection), and it is this dynamic regime on which we focus our study. Although our knowledge of free and forced convection is extensive (see Grossmann & Lohse (2000); Lohse & Shishkina (2024), for example), the mixed convection literature is comparatively small. For example, the first study focused on the role of shear in Rayleigh–Bénard (RB) convection was undertaken by Domaradzki & Metcalfe (1988), and since then, there have been numerous studies focused on the impacts of mixed convection of the vertical transport of heat in channel flows (Scagliarini et al. 2014, 2015; Pirozzoli et al. 2017), and atmospheric boundary layers (Moeng & Sullivan 1994; Salesky et al. 2017). A primary hallmark of mixed convection is that convective plumes generated near a solid boundary become aligned with the background flow within the interior of the flow, leading to large, coherent streamwise rollers, often referred to as superstructures.

Studies focused on the role of convection on the transport of inertial particles are relatively few. For example, see the recent work by Denzel et al. (2023), which focused on developing a stochastic model for the lifetime of small particles in an experimental chamber designed to study clouds in RB flow (Chang et al. 2016). To the authors’ knowledge, there may only be one study focused on the role of mixed convection in particle transport (Zaza & Iovieno 2024). In that work, the authors focused on the role of the particles in heat transfer throughout a turbulent boundary layer, and while the authors traversed a relatively large range of Richardson number (the key parameter quantifying the strength of convective turbulence to shear generated turbulence), they were restricted to relatively low Reynolds number ( $Re_\tau=180$ ), and ignored particle settling. In this work, we aim to describe a new mechanism by which strongly settling inertial particles are efficiently mixed through the boundary layer, in a way that is not present for either pure channel flow or free convection. Specifically, we are interested in the interactions between the near-boundary structures (associated with low-speed streaks and ejection events, see e.g. Wallace (2016)), and the interior superstructures generated by convective plumes, and how they couple to lead to vertical transport of isothermal inertial particles settling under the action of gravity. To investigate the dynamics, we use a series of coupled Eulerian-Lagrangian direct numerical simulations to simulate settling inertial particles channel flow ranging from free convection to pure shear.

2. Technical background

2.1. Carrier phase

In this paper, we use the NCAR Turbulence with Lagrangian Particles Model (Richter & Chamecki 2018) to simulate one-way coupled inertial particles emitted from the lower solid boundary in a turbulent closed-channel flow. This code has been validated and used in multiple studies focused on inertial particle settling and transport in turbulent boundary layers (Wang et al. 2019; Bragg et al. 2021; Gao et al. 2023; Grace et al. 2024). For the carrier phase, we use direct numerical simulations (DNS) to solve the three-dimensional, incompressible Navier–Stokes equations under the Boussinesq approximation in a turbulent channel flow set-up of streamwise length $L_x$ , spanwise extent $L_y$ and total height $2h$ . At the upper and lower boundaries, a no-slip boundary condition is enforced, while 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, $-{\rm d}P/{\rm d}x\gt 0$ Change to $-\partial P/\partial x > 0$ for consistency with figure 1 (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 and $\rho _a$ is the fluid density. Using the friction velocity, the channel half height $h$ and the applied temperature difference at the walls $\,\Delta T$ , we non-dimensionalise the velocity, spatial coordinates and the temperature field, resulting in the following set of non-dimensional equations:

(2.1) \begin{align} &\frac {{\rm D}\boldsymbol {u}}{{\rm D}t} = -\nabla p - {Ri}_\tau \theta \hat {\boldsymbol {z}} +\frac {1}{{Re}_\tau }\nabla ^2\boldsymbol {u} - \hat {\boldsymbol {x}}, \quad \frac {{\rm D}\theta }{{\rm D}t} = \frac {1}{{Pr\!Re_\tau }}\nabla ^2\theta ,\quad \nabla \cdot \boldsymbol {u} = 0, \end{align}

where the governing parameters of the carrier phase are

(2.2) \begin{equation} {Ri}_\tau = \frac {{Ra}}{{Pr\!Re}_\tau ^2},\quad {Ra} = \frac {g\alpha _\theta \,\Delta T(2h)^3}{\nu \kappa }, \quad Re_\tau = \frac {u_\tau h}{\nu },\quad {Pr} = \frac {\nu }{\kappa }. \end{equation}

Respectively, these are the friction Richardson number, the Rayleigh number, the friction Reynolds number and the Prandtl number. In these parameters, $\alpha _\theta$ is the isobaric thermal expansion coefficient. Throughout this work, we assume $Pr=0.715$ for all cases, where $\nu$ is the kinematic viscosity for dry air, and $\kappa$ is the thermal diffusivity of dry air. As mentioned previously, $Ri_\tau $ is an important parameter as it characterises the relative importance of turbulence generated by convection (through $Ra$ ) to that generated by friction along the solid boundary (through $Re_\tau )$ . Practically speaking, Pirozzoli et al. (2017) noted the appearance of coherent superstructures in both the temperature field and vertical fluctuating velocity field within the regime $1\lt Ri_\tau \lt 1000$ . In this work, we aim to investigate the regime $Ri_\tau\sim \mathcal {O}(10-100)$ as it allows us to attain computationally feasible values for $Re_\tau $ while also allowing for the formation of the coherent roll structures. Values for $Ra$ , $Re_\tau$ and $Ri_\tau$ used throughout this work can be found in table 1. In the table, the prefix of the case names correspond to the dynamic regime of the flow, i.e. CF (channel flow), FC (free convection) and MC (mixed convection). The suffix of the cases indicates particle specific information for that particular case, and this will be discussed in § 2.2. Finally, using an algebraic grid stretching procedure, we cluster grid points near the top and bottom boundaries so that $\,\Delta z_1 =0.8\nu /u_\tau$ , 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 mid-plane of the flow, where the subscript $\eta$ represents the grid spacing in terms of the local Kolmogorov scale. Note that, in Shishkina et al. (2010), for $Pr\sim \mathcal {O}(1)$ , the local Kolmogorov scale is the limiting factor in determining the resolution, and we have ensured that in all simulations, that $\,\Delta z_1\lt \eta$ as well.

Table 1. Cases discussed in this manuscript; CF stands for ‘channel flow’, MC stands for ‘mixed convection’ and FC stands for ‘free convection’. The suffix of each case follows the naming convention $-{St^+}-{Sv}^+$ . Parameter definitions can be found in the main text. All cases were run on a $512^3$ grid, with $Pr = 0.715$ and a reservoir concentration of $\mathcal{C}$ , except cases MC-0.5-0.05, MC-5-0.05, which had $\mathcal{C}$ due to computationally infeasible particle numbers in the low $Sv^+$ limit.

2.2. Dispersed phase

The applications of this work are towards coarse dust transport in the atmospheric surface layer. These dust particles range in size, but even coarse and giant grains (roughly 30–100 $\,\mathrm {\unicode{x03BC} m}$ ) are significantly smaller than the local Kolmogorov scale, which can be in the range of several millimetres. These particles are also significantly denser than the carrier phase, so we may ignore added mass and Basset-history forces, as well as two-way coupling and particle-particle interactions. Given these assumptions, we apply the point-particle approximation and apply the conservation of momentum for a rigid spherical particle subjected to linear hydrodynamic drag and gravity. The one-way coupled point-particle approach also has the added benefit that each particle is independent from each other particle, which removes the volume fraction as a governing parameter. This allows us to increase the number of particles to ensure convergence of the statistics of interest without affecting the flow. As a consequence, the results here would only be valid for physical systems with sufficiently low volume fraction.

The particles dynamics is characterised by two parameters; the Stokes number and the settling parameter, which quantify the importance of the particle inertia and the particle settling, respectively. The dynamical regime we are concerned with in this work has two fundamental processes working in concert (turbulence driven by mean shear and convection), which means we can consider different versions of Stokes numbers and settling parameters. First, we can define $Sv^+$ and $St^+$ , which are the settling number and Stokes number based on the friction velocity, $u_\tau$ , and the viscous time scale, $\nu /u_\tau ^2$ , respectively. These parameters characterise the particle inertia and settling with respect to the shear-driven turbulent scales. We can also define $Sv_*$ and ${St_*}$ , which are the same parameters instead based on the turbulent convective velocity scale $w_* = (g\alpha_\theta h \langle w'\theta '\rangle _s)^{1/3}$ (also known as the Deardorff convective velocity scale, where $\langle w'\theta '\rangle _s$ is the surface heat flux) and the large-scale turnover time, $h/w_*$ . These parameters are associated with the bulk flow motion, and characterise the large-scale convective motions of the flow. These parameters are defined as

(2.3) \begin{equation} {St}^+ =\frac {\tau _pu_\tau ^2}{\nu }, \; {Sv}^+ =\frac {v_g}{u_\tau },\; {St}_* =\frac {\tau _pw_*}{h}, \; {Sv}_* =\frac {v_g}{w_*}, \end{equation}

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

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

where $\rho _p$ is the particle density and $d_p$ is the particle diameter. The Stokes settling velocity is defined as $v_g = \tau _p g_p$ where $g_p$ is the gravitational acceleration applied to the particle. In simulations, $g_p$ need not be equivalent to $g$ (the gravitational acceleration applied to the fluid), thus allowing for turbulence, settling and inertial properties to be specified independently. The values for each parameter used throughout this work can be found in table 1. The suffix of each case name follows the naming convention - $St^+$ - $Sv^+$ , except for case MC-Sc, which is discussed below.

Using the friction scales to non-dimensionalise the particle equations of motion, we arrive at

(2.5) \begin{equation} {St}^+\frac {{\rm d}\boldsymbol {v}_p}{{\rm d}t} =\Psi \left (\boldsymbol {u}_f(\boldsymbol {x}_p(t),t)- \boldsymbol {v}_p\right ) - {Sv}^+\hat {\boldsymbol {z}} .\end{equation}

Here, $\boldsymbol {v}_p = (v_1,v_2,v_3)$ is the three-dimensional velocity vector for each particle, $\boldsymbol {x}_p$ is the location of each particle in space and $\boldsymbol {u}_f(\boldsymbol {x}_p(t),t) = (u,v,w)$ is the three-dimensional instantaneous flow velocity evaluated at the location of the particle. In the expressions above, $\Psi = 1 + 0.15{Re}_p^{0.687}$ is the Schiller-Neumann correction to the drag force, and $Re_p$ is the particle Reynolds number. As we keep the particle diameter fixed for all cases in this work, the particle Reynolds number remains small, meaning that $\Psi \approx 1$ .

Next, we must consider the boundary conditions for the particles, specifically the scheme by which we are to lift the particles into the domain. Here, we take an approach focused on simplicity; following Richter & Chamecki (2018) and Bragg et al. (2021), we add a Brownian-like jump term to the particle equations of motion. Using this term, particles take a discontinuous jump with zero mean and variance $\sqrt {2\kappa _p {\rm d}t}$ where, for practical purposes, ${\rm d}t$ is the timestep of the model and $\kappa _p$ is a parameter modelling the particle diffusivity. Mathematically, the location of the particle centroid is advanced according to the following equation:

(2.6) \begin{equation} {\rm d}\boldsymbol {x}_p = \boldsymbol {v}_p {\rm d}t + \left (\frac {Sc}{2}\right )^{-1/2} {\rm d}t^{1/2}{\rm d}\boldsymbol {\xi }, \end{equation}

where ${\rm d}\boldsymbol {\xi }$ is a Wiener process. Particles are initialised in a thin reservoir of thickness $D$ held at a fixed concentration beneath the lower solid boundary and are emitted through the lower surface. This approach effectively creates a Dirichlet condition for the particle concentration which is maintained throughout the simulation. For the upper boundary, particles reflect creating a Neumann no-flux condition. This model introduces several new parameters, specifically the particle Schmidt number, $Sc = \nu /\kappa _p$ and the fixed reservoir concentration, $\mathcal {C}$ , that are defined as

(2.7) \begin{equation} Sc=\frac {\nu }{\kappa _p},\quad \mathcal {C} = \frac {N_ph^2\nu }{u_\tau L_xL_yD}, \end{equation}

respectively. Here, $N_p$ is the number of particles contained in the reservoir. We note that the Brownian term is purely artificial and acts as a mechanism to diffuse particle into the domain where their dynamics is then governed by the fluid motion. We have investigated the role of Sc in the Appendix (see case MC-Sc), and we see that it only plays a minor role below the quasi-homogeneous layer, where the dynamics of the carrier phase takes over. We wish to re-iterate that some sort of emission mechanism is required to lift particles into the domain, and this is meant to emulate the role of more realistic (but significantly more complex) emission mechanisms, such as those discussed in Dupont et al. (2013). In this work, we hold $\mathcal {C}$ fixed for all cases except MC-0.5-0.05 and MC-5-0.05, where $\mathcal {C}$ is reduced by an order of magnitude for computational feasibility. Note that, under the one-way coupled approximation, the choice of $\mathcal {C}$ does not affect the results and is chosen to ensure statistical convergence. A schematic highlighting the salient features of the model is shown in figure 1.

Figure 1. A rectangular channel of dimensions $4\pi h\times 2\pi h\times 2h$ . The flow is periodic in the horizontal and is driven by a constant streamwise pressure gradient. No-slip boundary conditions are enforced at the top and bottom boundaries. The solid boundaries are held at a temperature difference $\Delta T=T_h-T_c$ . Particles are emitted randomly from a reservoir (fixed $\mathcal {C}\approx 4000$ ) into the domain through the bottom boundary. Particles reflect elastically off the upper boundary.

For case FC, a friction velocity cannot be defined, as there is no mean surface stress. Therefore, for that case specifically, we cannot define $St^+$ and $Sv^+$ . For the following results, any plots containing data plotted against $z^+$ (i.e. the height scaled in wall unts), data from FC will be plotted against $z^+$ from the CF case. Likewise, for CF, there is no associated heat flux, so a convective velocity cannot be defined meaning that $St_*$ and $Sv_*$ are undefined for that case.

3. Results

3.1. Primary and secondary suspension

We first examine how increasing $Ri_\tau$ induces the transition from the pure channel flow state (CF) to the free convection state (FC). Here we consider the effect of increasing $Ri_\tau$ on the vertical component of the turbulent kinetic energy, the mean temperature profiles and the vertical momentum flux in figure 2(a–c), respectively. As $Ri_\tau$ is increased, the intensity of the near-boundary fluctuations decreases while there is a corresponding increase in the intensity of fluctuations within the bulk, shown in figure 2(a). The increase in interior fluctuations occurs because of the appearance of large-scale convective plumes. These convective plumes result in a mechanism that efficiently re-distributes heat vertically, which results in the homogenisation of the mean temperature profiles, highlighted in figure 2(b). It is clear from these plots that there are two separate regions of the flow: the shear dominated region where the local maxima in the turbulent kinetic energy occur, and the bulk, which becomes progressively more important as $Ri_\tau$ increases. Interestingly, we see that there are no significant differences in the vertical momentum flux as $Ri_\tau$ increases, shown in figure 2(c). Throughout this work, we will refer to the quasi-homogeneous region of the flow, where the Reynolds stresses are thought to be roughly constant (Kunkel & Marusic 2006). Due to the relatively low Reynolds numbers of our DNS, the Reynolds stresses are not constant over a substantial region of the flow, so when we refer to this region, we are instead referring to where the Reynolds stresses are instead slowly varying. This will be discussed later. Finally, note we have omitted the momentum flux for the FC cases ( $Ri_\tau =\infty$ ), as there is no mean momentum flux in this case.

Figure 2. Horizontally averaged profiles of the vertical fluid velocity variance (a), the temperature profile (b) and the vertical momentum flux (c) for increasing values of $Ri_\tau$ . Note that the case FC ( $Ri_\tau =\infty$ ) is instead scaled by $w_*$ in (a), and is omitted from (c) for clarity.

These profiles serve to show that, as $Ri_\tau$ crosses a threshold of $\mathcal {O}(10)$ , the turbulent kinetic energy begins to shift from a pure channel flow like structure, to one that shares similarities to both a channel flow (local peaks in turbulent kinetic energy near the boundary) and free convection (local maximum in turbulent kinetic energy in the bulk). A similar transition was investigated in set-ups of the planetary boundary layer (Salesky et al. 2017), and DNS studies (Pirozzoli et al. 2017; Blass et al. 2020). Below, we highlight these structures in these regions of the flow in more detail.

Figure 3. Slices of the fluctuating vertical velocity at the mid-plane for pure channel flow (a,b), free convection (c,d) and mixed convection (e,f). Fluctuating velocities are normalised by $u_\tau$ in panels (a,b) and (e,f) and by $w_*$ (the convective velocity scale) in panels (c,d). The left column shows $x\!-\!z$ slices and the right column shows $x\!-\!y$ slices. Shaded contours in panels (c,f) show regions of $ w'\theta ' \gt 0.12\kappa \,\Delta T(2h)^{-1}{Ra}^{1/3}$ .

To examine the consequences of mixed convection on the suspension of settling particles, we first provide a comparison of the coherent structures present in the turbulent flow by examining the fluctuating vertical fluid velocity for a pure turbulent channel flow (CF; figure 3 a,b), free convection (FC; figure 3 c,d) and mixed convection (MC; figure 3 e,f). For each of these cases, the left column shows $y\!-\!z$ slices along the centre of the domain, while the right column shows $x\!-\!y$ slices at the mid-plane, indicated by small schematics at the top of the figure columns. Finally, for the mixed and free convection cases, we have indicated regions where $w'\theta ' \gt 0.12\kappa \,\Delta T(2h)^{-1}{Ra}^{1/3}$ with contour shading. This criterion is derived from the scaling relationship discussed in Pirozzoli et al. (2017), who studied turbulent channel flows under unstable stratification for similar $Re_\tau$ and $Ra$ . This criterion serves as a quantitative indicator of large positive turbulent heat fluxes in the domain. By isolating regions where the turbulent heat flux is larger than this value (the shaded regions), we can identify structures in the flow that exhibit very strong updrafts and downdrafts. Other criteria have been used for defining plume structures in the literature. For example, Van Der Poel (2015) and Huang et al. (2013) used similar but more restrictive criteria that included information of the relative temperature of the plume structures within high heat flux regions. We have compared our results with their criteria and found no substantial differences.

The CF case, figure 3(a,b), exhibits characteristic flow structures that vary in vertical scale across the domain, including small-scale instabilities associated with low-speed streaks near the solid boundaries, and larger features within the interior that scale with $h$ . Conversely, in the FC case, shown in figure 3(c,d), we can see evidence of organised convective cells in both the horizontal and vertical slices, which tend to scale with the full domain height $2h$ and dominate the interior motion. Importantly, there is much less activity near the solid boundary, as the boundary stresses induced by the convective plumes are not significant at this $Ra$ (Pirozzoli et al. 2017; Blass et al. 2020). The important insight is that MC, figure 3(e,f), shares aspects of both the CF and FC cases. For example, we can see large-scale structures within the interior of the MC case (strong updrafts and downdrafts), figure 3(e), reminiscent of the domain size convective cells from the FC case. However, in the horizontal, figure 3(f), these interior plumes becomes strongly aligned in the streamwise direction, with weaker fluctuations between the coherent roll structures. Furthermore, we can see that the near-boundary structures (also associated with strong heat fluxes), are qualitatively similar to the CF case. To summarise, MC exhibits large-scale convective interior plumes that scale with $2h$ , characteristic of FC, but also exhibits significant activity near the solid boundaries, characteristic of CF. However, due to the mean shear, the convective plumes align in the streamwise direction, creating large-scale superstructures.

Figure 4. As in figure 3. Contours are regions where $w'\theta ' \gt 0.12\kappa \,\Delta T(2h)^{-1}{Ra}^{1/3}$ , and are coloured based on the sign of the vertical fluid velocity (red is positive and blue is negative). Particles (not to scale) are overlaid highlighting their clustering behaviour.

We now highlight the role that these structures have in the transport of particles into the interior of the channel in MC, as compared with FC and CF. Figure 4 shows the snapshots of the high heat flux contours discussed in § 3, except here we colour the contours based on the direction of the vertical velocity fluctuation (red indicates positive fluid velocities while blue indicates negative fluid velocities). Overlaid are particles in a slab of non-dimensional thickness 0.02. It is clear from slices in the CF case, figure 4(a), that particles are ejected by the structures near the boundary layer, resulting in some clustering in the mid-plane, figure 4(b) (Lee & Lee 2019). When compared with the FC case, figure 4(c,d), particles are suspended much higher in the domain when they coincide with a strong updraft, shown in figure 4(c). However, these updrafts are much less efficient at generating ejections near the boundary, resulting in far fewer particles within the bulk, shown in figure 4(d). However, since both the near-boundary instabilities and interior convective motions are present in MC, figure 4(e,f), these structures can interact cooperatively, and we see a marked increase in the number of particles within the bulk, as well as strong spatial clustering within the longitudinally aligned updrafts. This observation suggests a coupling between the canonical near-wall ejections and the large-scale convective plumes in the bulk, which act cooperatively to lift particles far from the solid wall.

Our hypothesis is that in CF and MC, particles are initially suspended by the near-wall fluid ejection events, which we refer to as primary ejection. Next, if a particle undergoing primary ejection aligns with a large-scale plume in the bulk, it is entrained into the plume and experiences significant vertical transport, termed secondary suspension. Conversely, at these values of $St^+$ and $Sv^+$ (see table 1), particles that do not align with a convective plume after primary suspension instead settle back towards the lower boundary. Importantly, since the near-boundary instabilities in the FC case are much weaker for this $Ra$ (i.e. there is no mean shear to generate fluid ejections) the mechanism responsible for primary suspension is much weaker (as the interior convective instabilities are much less efficient at entraining them from this layer), resulting in far fewer particles suspended in the domain interior.

Figure 5. Profiles of mean concentration (a) and slabwise correlation coefficients conditioned on ejection events focusing on the constant flux region ( $z^+\lt 300$ ) (b) for increasing $Ri_\tau$ and fixed $St^+=5$ and $Sv^+=0.5$ . The dashed curve in panel (b) represents a profile of the magnitude of the fluid phase vertical momentum flux, normalised by its maximum. Note that panel (b) is plotted in terms of wall units, $z^+$ .

In figure 5(a), we examine the response of the global particle concentration to increasing $Ri_\tau$ . Figure 5(a) shows the horizontally and temporally averaged concentration profiles for fixed $St^+\approx 5$ and $Sv^+\approx 0.5$ as $Ri_\tau$ is increased from zero (CF) to infinity (FC). Here, we choose to define the Stokes number and the settling parameter based on the viscous scales, as it is expected that the particle dynamics of the near-wall region of the flow depend primarily on these scales. Importantly, we can see that as $Ri_\tau$ is increased, there is a significant increase in the bulk concentration above a height of roughly $z/h\approx 0.2$ . Recall that as $Ri_\tau$ increases, bulk convective plumes align with the direction of mean streamwise current, as demonstrated in visually figure 4(f). This effect leads to progressively larger differences in the bulk concentration, and even up to an order of magnitude when MC and CF are compared at the mid-plane (a height of $z/h=1$ ). Upon further increasing $Ri_\tau$ , here achieved with FC, the concentration within the interior region decreases as $Ri_\tau$ relative to MC, while the mid-channel slope of the MC and FC profiles remains similar in character. These results imply that that the efficacy of secondary suspension is linked to the coherence of the bulk superstructures, specifically their alignment with the mean shear.

In spite of the differences between concentration profiles in the bulk as $Ri_\tau$ is increased, the concentration profiles for all cases (aside from FC) are coincident below $z/h\approx 0.2$ . This is a reflection of the hypothesis that fluid ejections near-wall region are responsible for primary suspension of particles, and these fluid ejections maintain a similar character until $Ri_\tau$ is large enough that free convection takes over everywhere in the domain. This suggests that, in the shear dominated region of the flow (i.e. where the flow share strong similarities to a channel flow), the particle dynamics may be independent of $Ri_\tau$ . To illustrate this point, we use a conditional averaging technique to show slabwise correlation coefficients for each case discussed in table 1. We define the conditional correlation coefficient as

(3.1) \begin{equation} R(\alpha ,\beta |\gamma ) = \frac {\langle \alpha \beta |\gamma \rangle }{\sqrt {\langle \alpha ^2|\gamma \rangle \langle \beta ^2|\gamma \rangle }}, \end{equation}

where $\alpha$ and $\beta$ are dummy variables used for demonstration in (3.1), and $\langle \cdot \rangle$ indicates a slabwise ensemble average. Figure 5(b) show profiles of the correlation coefficient between the vertical particle velocity and the Reynolds stresses conditioned on ejection events (i.e. $u'\lt 0$ and $w'\gt 0$ , where a prime indicates a fluctuating quantity) to a height of 300 wall units. This region highlights the quasi-homogeneous region (Bragg et al. 2021; Pope 2000), where the vertical momentum flux achieves its maximum magnitude. Following the notation of Salesky et al. (2017), we use $(u'w')^{II}$ to represent this conditioning. Additionally, we plot the magnitude of the fluid-phase momentum flux normalised by its maximum value.

Profiles in figure 5(b) show a rapid increase in the correlation between vertical particle velocities and ejections as the vertical momentum flux increases, suggesting that below this region, the artificial Brownian diffusion controls the particle dynamics, but is quickly replaced by that of fluid ejections. Above roughly $z^+=10$ , the correlations for all cases except FC are nearly coincident within the region of strong vertical momentum fluxes, suggesting an independence to variation of $Ri_\tau$ . The important insight here is that when the shear is strong enough to produce the ejection events, then these ejections govern the transport of the particles through the quasi-homogeneous region of the flow, and this process is independent of $Ri_\tau$ (and is not influenced by the artificial diffusion).

As particles are lifted into the bulk, secondary suspension takes over. As demonstrated in figure 5(a), the efficacy of the secondary suspension relies on the degree to which the superstructures are aligned, which is a non-monotonic function of $Ri_\tau$ (Blass et al. 2020). Importantly, the interior entrainment by secondary suspension remains weak until $Ri_\tau$ becomes large enough that the interior convective plumes align – a process that happens very rapidly as $Ri_\tau$ surpasses roughly 10. Then, as $Ri_\tau$ becomes increasingly large, the shear becomes relatively weak, and ejection events near the solid boundary are no longer present, resulting in an increasingly weak primary suspension mechanism. As there are far fewer particles made available to the bulk by primary suspension, there is an associated decrease in the number concentration within the interior.

The main conclusion is that primary suspension is a function of the magnitude of the mean shear and occurs independently of $Ri_\tau$ , at least up to roughly 50. Primary suspension makes particles available to the bulk, and it is not until the bulk plumes are aligned by the mean shear that secondary suspension can occur. As such, secondary suspension is a strong function of $Ri_\tau$ .

3.2. Sensitivity to particle inertial and settling

We now consider how the efficacy of primary suspension is affected by particle inertia and settling. First, we consider the horizontally and temporally averaged concentration profiles at fixed $Ri_\tau \approx 50$ as $St^+$ and $Sv^+$ are varied independently. We have chosen to investigate the dynamics in terms of the frictional parameters, as we assume that the particle dynamics is governed by those scales in this region of the flow at this value of $Ri_\tau$ . We consider particles covering three orders of magnitude of $St^+$ , highlighted in table 1. For reference, particles in the range ${St^+}=0.5{-}50$ have Stokes numbers based on the local Kolmogorov scale in the range, $St_\eta = 0.1{-}10$ . We also choose values of $Sv^+$ which we qualitatively categorise into weakly settling ( $Sv^+\approx 0.05$ ) and moderate settling ( $Sv^+\approx 0.5$ ), which have settling parameters based on the local Kolmogorov scale in the range $Sv_\eta \approx 0.1{-}1$ . These values correspond to the range where it is known that particle inertia and settling become dynamically relevant (Good et al. 2014). We also include a case that we can categorise into strongly settling ( $Sv^+\approx 1$ ), where we expect that fluid ejections are simply not strong enough to overcome the mean gravitational drift.

Figure 6. Profiles of mean concentration for varied $St^+$ at fixed $Sv^+$ (a) and vice versa (b) for fixed $Ri_\tau$ .

Figure 6(a) shows the concentration profiles at fixed $Sv^+$ traversing the entire range of $St^+$ considered in this work. The results in figure 6(a) show that within the bulk, where the dynamics is strongly impacted by the coherent roll structures, there is a weak dependence of the magnitude of the concentration as $St^+$ increases. However, in figure 6(b), which shows the concentration profiles at two values of $St^+$ for each value of $Sv^+$ considered, there is a much stronger dependence on $Sv^+$ at fixed $St^+$ . Specifically, by increasing $Sv^+$ by roughly an order of magnitude from 0.05 to 0.5, there is over an order of magnitude drop in the particle concentration across the entire domain. Importantly, we see that as $Sv^+$ becomes $\mathcal {O}(1)$ , we see a massive decrease in the bulk concentration, as fluid ejections do not supply enough upward momentum to suspend such strongly settling particles. From these plots, we can surmise that particle inertia plays a minor role in the dynamics, while settling has a more influential role in the dynamics. We hypothesise that the overall decrease in the particle concentration in the bulk is not wholly a consequence of secondary suspension, but is instead influenced by the efficacy of primary suspension.

Figure 7. Profiles of the vertical particle flux conditioned on fluid phase ejections for varied $St^+$ at fixed $Sv^+$ (a) and vice versa (b) at fixed $Ri_\tau$ . Note that MC- 5-1 is omitted from this plot.

To see why this is, we consider the average vertical particle flux conditioned on ejection events, written as $\langle \Phi | (u'w' )^{II}\rangle$ . To compute this quantity, we count number of particles in a slab of thickness $dv$ that are sampling ejection events, and multiply by the vertical velocities of those particles. From this, we subtract a diffusive flux due to the artificial Brownian term in the particle equations of motion (see Bragg et al. 2021), although this does not make a substantial difference.

Figure 7(a) shows the vertical flux when $St^+$ is varied at fixed $Sv^+$ , where we have chosen to focus on the region of the flow within 300 wall units as this is the region of the flow undergoing primary suspension. Importantly, the average vertical flux is positive for particles sampling fluid ejections, which agrees with our primary suspension hypothesis in § 3.1. Here, we note that, until $St^+= 50$ , there is almost no difference in the vertical flux of the particles. It is known that as $St^+$ increases, particles become increasingly insensitive to the high-frequency flow fluctuations, often referred to as inertial filtering (Bec et al. 2006). This is also evident in figure 7(b), as we can see that there is a significant sensitivity of $\langle \Phi | (u'w' )^{II}\rangle$ to $Sv^+$ at fixed $St^+$ . When $Sv^+$ is small, there is a significant vertical flux, as particles can be easily suspended by the fluid ejections, but by increasing $Sv^+$ by an order of magnitude, there is a large decrease in the vertical flux of particles within this region of the flow. These results have significant implications for the rate at which particles may be supplied to the bulk motions. They suggest that, for primary suspension, the particle Stokes number has a very small role to play, acting to decrease the suspension rate due to inertial filtering, but the dynamics is more strongly influenced by $Sv^+$ , with the efficacy of primary suspension being very high for low $Sv^+$ , but decreasing quickly as $Sv^+$ increases. Note we have omitted the conditional flux for $Sv^+=1$ , as it is significantly lower than the other cases.

In summary, both particle inertia and particle settling tend to reduce the efficacy of the primary suspension. Increasing $St^+$ leads to inertial filtering, meaning that particles become insensitive to the fluid ejections responsible for primary suspension, however, this occurs slowly as $Sv^+$ is increased. Alternatively, by increasing $St^+$ the mean downward drift increases, relative to the upward ejection velocities, leading to a decrease in the conditional flux. These effects working in concert lead to decreases in concentration of the particles within the bulk. However, once particles are suspended into the bulk by the primary suspension mechanism, the particle inertia and settling appear to play a secondary role in the dynamics, as their relaxation time and settling velocities are both much smaller relative to the dominant motions in the bulk, quantified by $St_*$ and $Sv_*$ .

4. Conclusions and discussion

In this work, we have provided the first demonstration of the role of turbulent mixed convection in the suspension of strongly settling Lagrangian particles. Through the use of coupled Eulerian–Lagrangian DNS of turbulent mixed convection, achieved here with $Ra\approx 10^7$ and $Re_\tau\approx 500$ , we provided evidence for a distinct, multi-scale mechanism that can serve as an efficient means of inertial particle suspension. Our observations (at $Ri_\tau \approx 50$ ) suggest that the action of the near-boundary ejections work in concert with streamwise aligned interior plumes to lead to an efficient global suspension mechanism for Lagrangian particles.

By considering correlation coefficients of vertical particle velocities conditioned on fluid ejections (typically associated with canonical low-speed streaks and hairpin vortices), we showed that particles are ejected near the solid boundary, which we termed primary suspension. We showed that primary suspension mechanism was independent of $Ri_\tau$ up to at least 50, and relied on the presence of strong fluid ejections near the boundary. We hypothesised that this process acts as a mechanism to supply the bulk region of the flow with particles, and was relatively insensitive to particle inertia (at least up to $St^+=50$ ), but was much more sensitive to particle settling. We showed that strong settling leads to a lower upward particle flux by ejections, which results in few particles supplied to the bulk from below. This led to a lower concentration of particles that could be entrained the coherent superstructures in the bulk, which we termed secondary suspension. Importantly, secondary suspension was only efficient when there was strong alignment of the convective plumes and is effectively absent in the limiting cases (pure channel flow and free convection), and when convection is present, but weak enough that convective plumes do not become aligned (observed at $Ri_\tau \approx 5.9$ here). Therefore, while primary suspension was found to be independent of $Ri_\tau$ , secondary suspension was shown to be strongly dependent on this parameter.

These results have important implications for both fundamental and applied studies on the influence of mixed convection on particles, including a potentially strong influence on clustering and dispersion, as well as determining particle residence times. Moreover, it is known that fluid structures responsible for secondary suspension appear ubiquitously at the field scale in the planetary boundary layer (Moeng & Sullivan 1994; Salesky et al. 2017), so more work should take place to understand these cooperative dispersed phase transport mechanisms in the natural environment. For example, as we have demonstrated, secondary suspension occurs when $Ri_\tau \sim 10$ − $100$ , meaning that, for field-scale Reynolds numbers (which are much larger than what a DNS can feasibly simulate), the required Rayleigh number must also be very large. As these superstructures have been documented in numerical studies of the planetary boundary layer (see Salesky et al. (2017), for example), it is plausible that the mechanism described in this work will scale to the field. It has also been hypothesised that, at large Reynolds number, the quasi-homogeneous region (where $\langle u'w'\rangle \approx u_\tau ^2$ , which is the dynamical region responsible for primary suspension) scales with the boundary layer height (see Kunkel & Marusic (2006), for example), it is plausible that primary ejection may occur at the field scale. Interestingly, this work also suggests that, as long as there is any primary suspension mechanism (such as jet and spume droplets from surface wave breaking in the marine atmospheric surface layer (Veron 2015), for example), secondary suspension can still lead to enhanced mixing at altitude.

Moreover, at higher $Ra$ , characteristic of the atmospheric surface layer, it is plausible that there is an associated increase in the convective velocity scale within the bulk, $w_*$ , since $\langle w'\theta '\rangle _{s}$ is an increasing function of $Ra$ in this regime. So, in fact, if particles can be ejected via any primary suspension mechanism, their secondary suspension may be even more significant in the field, as the importance of inertia and settling may further decrease.