We investigate the role of slippery boundaries, quantified by the Navier boundary friction coefficient $\beta$, in regulating heat transport and flow structures in rotating Rayleigh–Bénard convection. Owing to the Ekman pumping effect arising from viscous boundary layers that is intensified with increasing boundary friction, it is found that the properties of global heat transport exhibit two distinct parameter regimes separated by a transitional Rayleigh number ($ \textit{Ra}_t$). In the rotation-dominated regime ($ \textit{Ra} \lt \textit{Ra}_t$), enhanced viscous friction increases the efficiency of Ekman pumping, significantly elevating the Nusselt number and lowering the convection onset threshold. Conversely, in the buoyancy-dominated regime ($ \textit{Ra} \gt \textit{Ra}_t$), boundary-induced viscous dissipation suppresses convective motions, thereby reducing heat transport. Large-scale vortices (LSVs), prevalent under free-slip conditions, progressively dissipate as $\beta$ increases, revealing that viscous friction disrupts the inverse energy cascade from baroclinic to barotropic modes. Through kinetic energy partitioning analysis, the transition between quasi-two-dimensional and three-dimensional turbulent states is identified, with the parameter $\beta _{\textit{cr}}$ following a generic scaling relation on the Prandtl (Pr) and Ekman (Ek) numbers $\beta _{\textit{cr}}\sim \textit{Pr}^{-0.67}\textit{Ek}^{-1.18}$. This relation enables us to predict LSV emergence across different parameter spaces. Furthermore, it is reported that the heat-transport scaling exponent, the convection onset and the partitioning of kinetic energy between barotropic and baroclinic components undergo a smooth flow transition at $\beta _{\textit{cr}}$. These results also indicate a direct correlation between Ekman pumping efficacy and the friction coefficient $\beta$, demonstrating that controlling boundary friction can modulate global transport properties and reshape flow structures.

We present three-dimensional velocity gradient statistics from turbulent Rayleigh–Bénard convection experiments in a horizontally extended cell of aspect ratio 25, a paradigm for mesoscale convection with its organisation into large-scale patterns. The Rayleigh number ${\textit{Ra}}$ ranges from $3.7 \times 10^5$ to $4.8 \times 10^6$, the Prandtl number ${\textit{Pr}}$ from 5 to 7.1. Spatio-temporally resolved volumetric data are reconstructed from moderately dense Lagrangian particle tracking measurements. All nine components of the velocity gradient tensor from the experiments show good agreement with those from direct numerical simulations, both conducted at ${\textit{Ra}} = 1 \times 10^6$ and ${\textit{Pr}} = 6.6$. As expected, with increasing ${\textit{Ra}}$, the flow in the bulk approaches isotropic conditions in the horizontal plane. The focus of our analysis is on non-Gaussian velocity gradient statistics. We demonstrate that statistical convergence of derivative moments up to the sixth order is achieved. Specifically, we examine the probability density functions (PDFs) of components of the velocity gradient tensor, vorticity components, kinetic energy dissipation and local enstrophy at different heights in the bottom half of the cell. The probability of high-amplitude derivatives increases from the bulk to the bottom plate. A similar trend is observed with increasing ${\textit{Ra}}$ at fixed height. Both indicate enhanced small-scale intermittency of the velocity field. We also determine derivative skewness and flatness. The PDFs of the derivatives with respect to the horizontal coordinates are found to be more symmetric than the ones with respect to the vertical coordinate. The conditional statistical analysis of the velocity derivatives with respect to up-/down-welling regions and the rest did not display significant difference, most probably due to the moderate Rayleigh numbers. Furthermore, doubly logarithmic plots of the PDFs of normalised energy dissipation and local enstrophy at all heights show that the left tails follow slopes of 3 / 2 and 1 / 2, respectively, in agreement with numerical results. In general, the left tails of the dissipation and local enstrophy distributions show higher probability values with increasing proximity towards the plate, in comparison with those in the bulk.

Turbulence is an out-of-equilibrium flow state that is characterised by non-zero net fluxes of kinetic energy between different scales of the flow. These fluxes play a crucial role in the formation of characteristic flow structures in many turbulent flows encountered in nature. However, measuring these energy fluxes in practical settings can present a challenge in systems other than the case of unrestricted turbulence in an idealised periodic box. Here, we focus on rotating Rayleigh–Bénard convection, being the canonical model system to study geophysical and astrophysical flows. Owing to the effect of rotation, this flow can yield a split cascade, where part of the energy is transported to smaller scales (direct cascade), while another fraction is transported to larger scales (inverse cascade). We compare two different techniques for measuring these energy fluxes throughout the domain: one based on a spatial filtering approach and an adapted Fourier-based method. We show how one can use these methods to measure the energy flux adequately in the anisotropic, aperiodic domains encountered in rotating convection, even in domains with spatial confinement. Our measurements reveal that in the studied regime, the bulk flow is dominated by the direct cascade, while significant inverse cascading action is observed most strongly near the top and bottom plates, due to the vortex merging of Ekman plumes into larger flow structures.

Turbulent Rayleigh–Bénard convection in an extended layer of square cross-section with moderate aspect ratio $L/H=8.6$ ($L$ is the length of the cell, $H$ is its height) is studied numerically for Rayleigh numbers in the range ${\textit{Ra}}= 10^6{-}10^8$. We focus on the influence of different types of boundary conditions, including asymmetrical ones, on the characteristics of Rayleigh–Bénard convection with and without an immersed freely floating body. Convection without a floating body is characterised by the formation of stable thermal superstructures with preferred location. The crucial role of the symmetry of the boundary conditions is revealed. In the case of thermal boundary conditions of different types at the upper and lower boundaries, the flow pattern in Rayleigh–Bénard convection has a regular shape. The immersed body makes random wanderings and actively mixes the fluid, preventing the formation of superstructures. The mean flow structure with an immersed body is similar for all combinations of boundary conditions except for the case of a fixed heat flux at both boundaries. The floating disk does not change the tendency of turbulent convection to form a circulation of the maximal available scale under symmetric Neumann-type conditions. The type of boundary conditions has a weak influence on the Nusselt and Reynolds number values, significantly changing the ratio of the mean and fluctuating components of the heat flux. As the Rayleigh number increases, the motions of the body become more intensive and intermittent. The increase of $Ra$ also changes the structure of the mean flow without the body but the additional mixing provided by the floating body preserves the flow structure.

Using direct numerical simulations, we systematically investigate the inner-layer turbulence of a turbulent vertical buoyancy layer (a model for a vertical natural convection boundary layer) at a constant Prandtl number of $0.71$. Near-wall streaky structures of streamwise velocity fluctuations, synonymous with the buffer layer streaks of canonical wall turbulence, are not evident at low and moderate Reynolds numbers (${\textit{Re}}$) but manifest at high ${\textit{Re}}$. At low ${\textit{Re}}$, the turbulent production in the near-wall region is negligible; however, this increases with increasing ${\textit{Re}}$. By using domains truncated in the streamwise, spanwise and wall-normal directions, we demonstrate that the turbulence production in the near-wall region at moderate and high ${\textit{Re}}$ is largely independent of large-scale motions and outer-layer turbulence. On a fundamental level, the near-wall turbulence production is autonomous and self-sustaining, and a well-developed bulk is not needed to drive the inner-layer turbulence. Near-wall streaks are also not essential for this autonomous process. The type of thermal boundary condition only marginally influences the velocity fluctuations, revealing that the turbulence dynamics are primarily governed by the mean-shear induced by the buoyancy field and not by the thermal fluctuations, despite the current flow being solely driven by buoyancy. In the inner layer, the spanwise wavelength of the eddies responsible for positive shear production is remarkably similar to that of canonical wall turbulence at moderate and high ${\textit{Re}}$ (irrespective of near-wall streaks). Based on these findings, we propose a mechanistic model that unifies the near-wall shear production of vertical buoyancy layers and canonical wall turbulence.

Turbulent convection under strong rotation can develop an inverse cascade of kinetic energy from smaller to larger scales. In the absence of an effective dissipation mechanism at the large scales, this leads to the pile up of kinetic energy at the largest available scale, yielding a system-wide large-scale vortex (LSV). Earlier works have shown that the transition into this state is abrupt and discontinuous. Here, we study the transition to the inverse cascade at Ekman number ${Ek}=10^{-4}$ and using stress-free boundary conditions, in the case where the inverse energy flux is dissipated before it reaches the system scale, suppressing the LSV formation. We demonstrate how this can be achieved in direct numerical simulations by using an adapted form of hypoviscosity on the horizontal manifold. We find that, in the absence of the LSV, the transition to the inverse cascade becomes continuous. This shows that it is the interaction between the LSV and the background turbulence that is responsible for the earlier observed discontinuity. We furthermore show that the inverse cascade in absence of the LSV has a more local signature compared with the case with LSV.

We report our finding from direct numerical simulations that polygonal cell structures are formed by inertial particles in turbulent Rayleigh–Bénard convection in a large aspect ratio channel at Rayleigh numbers of $10^6, 10^7$ and $10^8$, and Prandtl number of 0.7. The settling of small particles modified the flow structures only through momentum interactions. From the simulations performed for various sizes and mass loadings of particles, we discovered that for small- and intermediate-sized particles, cell structures such as square, pentagonal or hexagonal cells were observed, whereas a roll structure was formed by large particles. As the mass loading increased, the sizes of the cells or rolls decreased for all particle sizes. The Nusselt number increased with the mass loading of intermediate and large particles, whereas it decreased with the mass loading of small particles compared with the value for particle-free convection. A detailed investigation of the effective feedback forces of the settling particles inside the hot and cold plumes near the walls revealed that the feedback forces break the up–down symmetry between the hot and cold plumes near the surfaces. This enhances the hot plume ascent while not affecting the cold plume, which leads to the preferred formation of cellular structures. The energy budget analysis provides a detailed interaction between particles and fluid, revealing that the net energy is transferred from the fluid to particles when the particles are small, while settling intermediate and large particles drag the fluid so strongly that energy is transferred from particles to fluid.

Motivated by the need for a better understanding of the melting and stability of floating ice bodies, we experimentally investigated the melting of floating ice cylinders. Experiments were carried out in a tank, with ice cylinders with radii between 5 and 12 cm, floating horizontally with their axis perpendicular to gravity. The water in the tank was at room temperature, with salinities ranging from 0 to 35 g l−1. These conditions correspond to Rayleigh numbers in the range 10$^5\lesssim$ Ra $\lesssim$ 10$^9$. The relative density and thus the floating behaviour was varied by employing ice made of H$_2$O–D$_2$O mixtures. In addition, we explored a two-layer stable stratification. We studied the morphological evolution of the cross-section of the cylinders and interpreted our observations in the context of their interaction with the convective flow. The cylinders only capsize in fresh water but not when the ambient is saline. This behaviour can be explained by the balance between the torques exerted by buoyancy and drag, which change as the cylinder melts and rotates. We modelled the oscillatory motion of the cylinders after a capsize as a damped nonlinear oscillator. The downward plume of the ice cylinders follows the expected scalings for a line-source plume. The plume’s Reynolds number scales with Rayleigh number in two regimes, namely Re $\propto$ Ra$^{1/2}$ for Ra $\lt \mathcal{O}(10^7)$ and Re $\propto$ Ra$^{1/3}$ for Ra $\gt \mathcal{O}(10^7)$, and the heat transfer (non-dimensional as Nusselt number) scales as Nu $\propto$ Ra$^{1/3}$. Although the addition of salt substantially alters the solutal, thermal and momentum boundary layers, these scaling relations hold irrespectively of the initial size or the water salinity. While important differences exist between our experiments and real icebergs, our results can qualitatively be connected to natural phenomena occurring in fjords and around isolated icebergs, especially with regard to the melting and capsizing behaviour in stratified waters.

We study convection in a volumetrically heated fluid which is cooled from both plates and is under rotation through the use of direct numerical simulations. The onset of convection matches similar systems and predictions from asymptotic analysis. At low rotation rates, the fluid becomes more organised, enhancing heat transport and increasing boundary layer asymmetry, whereas high rotation rates suppress convection. Velocity and temperature statistics reveal that the top unstably stratified boundary layer exhibits behaviour consistent with other rotating convective systems, while the bottom boundary shows a unique interaction between unstable stratification and Ekman boundary layers. Additional flow statistics such as energy dissipation are analysed to rationalise the flow behaviour.

Plumes generated from a point buoyant source are relevant to hydrothermal vents in lakes and oceans on and beyond Earth. They play a crucial role in determining heat and material transport and thereby local biospheres. In this study, we investigate the development of rotating point plumes in an unstratified environment using both theory and numerical simulations. We find that in a sufficiently large domain, point plumes cease to rise beyond a penetration height $h_{{f}}$, at which buoyancy flux from the heat source is leaked laterally to the ambient fluid. The height $h_{{f}}$ is found to scale with the rotational length scale $h_{ \!{ f}}\sim L_{ \!\textit{ rot}}^p\equiv ({F_0}/{f^3})^{{1}/{4}},$ where $F_0$ is the source buoyancy flux, and $f=2\varOmega$ is the Coriolis parameter ($\varOmega$ is the rotation rate). In a limited domain, the plume may reach the top boundary or merge with neighbouring plumes. Whether rotational effects dominate depends on how $L_{\textit{rot}}^{p}$ compares to the height of the domain $H$ and the distance between the plumes $L$. Four parameter regimes can therefore be identified, and are explored here through numerical simulation. Our study advances the understanding of hydrothermal plumes and heat/material transport, with applications ranging from subsurface lakes to oceans in icy worlds such as Snowball Earth, Europa and Enceladus.

The momentum dispersion model for flows in isotropic porous media has been validated and successfully applied by Rao & Jin (2022, J. Fluid Mech., vol. 937, A17). However, the anisotropic coupled models concerning heat–fluid–solid interactions in turbulent forced convection requires further development. This research proposes various anisotropic physical coefficient tensors to model the total drag ${R}_{i}$, interphase energy resistance $H$, momentum dispersion and thermal dispersion accounting for both anisotropic and isotropic scenarios. The effective physical coefficients of the Darcy–Forchheimer equation regarding ${R}_{i}$ are adapted to accommodate anisotropy. The heat transfer coefficient $h$ between the solid and fluid, despite being a scalar, is also required to depend on the local flow direction in anisotropic cases. Two scaling laws of $h$ with respect to a local Reynolds number ${\textit{Re}_{K}}$ are found: $h\sim \textit{Re}_K^2$ for the Darcy regime, and $h\sim \textit{Re}_{K}^{1/2}$ for the Forchheimer regime, with a transition at ${\textit{Re}_{K}}\sim 1$. The influence of momentum and thermal dispersions, along with the modelling errors of ${R}_{i}$ and $H$ originating from heterogeneity, are approximated using a second-order pseudo-stress tensor and a pseudo-flux vector, respectively. The effective viscosity and thermal diffusivity tensors are simplified into longitudinal and transverse components using tensor symmetries, and are assumed to rely mainly on another local Reynolds number ${\textit{Re}_{d}}$. Both components of the effective viscosity are positive in isotropic cases, whereas the longitudinal component may be negative in anisotropic cases, mainly serving as a compensation of overestimated drag. The coupled models are applied to simulate turbulent forced thermal convection in porous media with one or two length scales across a wide range of Reynolds numbers. The comparisons with direct numerical simulations results imply that the coupled macroscopic models can accurately predict not only statistically stationary distributions but also real-time changes in velocity and temperature.

A model for obtaining scaling laws for Rayleigh–Bénard convection (RBC) at high Rayleigh numbers in tall, slender cells (cells with low aspect ratio, $\varGamma = d/H \ll 1$) is presented. Traditional RBC ($\varGamma \gtrsim 1$) is characterised by large-eddy circulation scaling with the height of the cell, a near-isothermal core and almost all the thermal resistance provided at the horizontal walls. In slender RBC cells, on the other hand, away from the horizontal walls, tube-like convection with eddies scaling with the tube diameter and a linear temperature gradient driving the convective flow is present. The crux of our approach is to split the cell into two components: (i) ‘wall convection’ near the top and bottom horizontal walls and (ii) ‘tube convection (TC)’ in the central part away from the walls. By applying the scaling relations for both wall convection and TC, and treating the total thermal resistance as a sum of their contributions, unified scaling relations for Nusselt number, Reynolds number and mean vertical temperature gradient in slender RBC cells are developed. Our model is applicable for high enough Rayleigh numbers, such that convection both at the wall and in the tube are turbulent. Our model predictions compare well with the data from various studies in slender RBC cells where these conditions are satisfied. In particular, the effects of $\varGamma$ and Prandtl number are well captured. We propose a scaled aspect ratio using which we obtain ‘universal’ correlations for the heat flux and for the fractional temperature drop in the tube that include the effects of Rayleigh and Prandtl numbers. The profiles of suitably scaled horizontal and vertical velocity fluctuations, along with estimates for boundary layer thickness near the horizontal walls, and the radial distribution of the velocity fluctuations in the tube part are also presented.

We simulate thermal convection in a two-dimensional square box using the no-slip condition on all boundaries, and isothermal bottom and top walls, and adiabatic sidewalls. We choose 0.1 and 1 for the Prandtl number $Pr$ and vary the Rayleigh number $Ra$ between $10^6$ and $10^{12}$. We particularly study the temporal evolution of integral transport quantities towards their steady states. Perhaps not surprisingly, the velocity field evolves more slowly than the thermal field, and its steady state – which is nominal in the sense that large-amplitude low-frequency oscillations persist around plausible averages – is reached exponentially. We study these oscillation characteristics. The transient time for the velocity field to achieve its nominal steady state increases almost linearly with the Reynolds number. For large $Ra$, the Reynolds number itself scales almost as $Ra^{2/3}\, Pr^{-1}$, and the Nusselt number as $Ra^{2/7}$.

Studying rotating convection under geo- and astrophysically relevant conditions has proven to be extremely difficult. For the rotating Rayleigh–Bénard system, van Kan et al. (J. Fluid Mech., vol. 1010, 2025,A42)have now been able to massively extend the parameter space accessible by direct numerical simulations. Their progress relies on a rescaling of the governing Boussinesq equations, which vastly improves numerical conditioning (Julien et al., arXiv:2410.02702). This opens the door for investigating previously inaccessible dynamical regimes and bridges the gap to the asymptotic branch of rapidly rotating convection.

We perform direct numerical simulations of centrifugal convection with an oscillating rotational velocity of small amplitude to study the effects of oscillatory boundary motion. The oscillation period is the main control parameter, with its range corresponding to a Womersley number in the range $1\lt Wo\lt 300$. Oscillating boundaries generate a circumferential shear flow, which significantly inhibits heat transfer, with maximum suppression $87\,\%$ observed in the present parameter space. Through analysis of the background flow, we find that as the oscillation period increases, the increasing penetration depth of the oscillation and weakening local shear strength result in non-monotonic changes in heat transfer. Under high-frequency oscillation, the characteristic length scale of the viscous layer induced by the oscillation is smaller than the convective length scale, and shear manifests primarily as a continuous suppression of the boundary layer. In contrast, under low-frequency oscillation, the shear flow covers the entire region but with weak strength. The suppression effect of such shear flow exhibits periodicity, leading to alternating phases of convection inhibition and convection generation. The present findings explore the physical mechanisms behind the suppression of convective heat transfer by oscillation, and offer a new strategy for controlling convection systems, with potential implications for both fundamental research and industrial applications.

Direct numerical simulations are performed for turbulent forced convection in a half-channel flow with a wall oscillating either as a spanwise plane oscillation or to generate a streamwise travelling wave. The friction Reynolds number is fixed at $Re_{\tau _0} = 590$, but the Prandtl number $Pr$ is varied from 0.71 to 20. For $Pr\gt 1$, the heat transfer is reduced by more than the drag, 40 % compared with 30 % at $Pr=7.5$. This outcome is related to the different responses of the velocity and thermal fields to the Stokes layer. It is shown that the Stokes layer near the wall attenuates the large-scale energy of the turbulent heat flux and the turbulent shear stress, but amplifies their small-scale energy. At higher Prandtl numbers, the thinning of the conductive sublayer means that the energetic scales of the turbulent heat flux move closer to the wall, where they are exposed to a stronger Stokes layer production, increasing the contribution of the small-scale energy amplification. A predictive model is derived for the Reynolds and Prandtl number dependence of the heat-transfer reduction based on the scaling of the thermal statistics. The model agrees well with the computations for Prandtl numbers up to 20.

In rotating convection, analysis of heat transfer reveals a distinct shift in behaviour as the system transitions from a steep scaling regime near the onset of convection to a shallower scaling at higher Rayleigh numbers ($Ra$), irrespective of whether the top and bottom plates have stress-free, no-slip or no boundaries (homogeneous convection). However, while most research on this transition focuses on no-slip boundary conditions, geophysical and astrophysical flows commonly involve stress-free and homogeneous convection models as well. This study delves into the transition from the rapidly rotating regime to the non-rotating one with both stress-free and homogeneous models, leveraging direct numerical simulations (DNS) and existing literature data. Our findings unveil that for stress-free boundary conditions, the transitional Rayleigh number ($Ra_T$) exhibits a relationship $Ra_T\sim Ek^{-12/7}$, whereas for homogeneous rotating convection, $Ra_T\sim Ek^{-2}\, Pr$, where $Ek$ denotes the Ekman number, and $Pr$ denotes the Prandtl number. Both of these relationships align with the data obtained through DNS.

Geophysical and astrophysical fluid flows are typically driven by buoyancy and strongly constrained at large scales by planetary rotation. Rapidly rotating Rayleigh–Bénard convection (RRRBC) provides a paradigm for experiments and direct numerical simulations (DNS) of such flows, but the accessible parameter space remains restricted to moderately fast rotation rates (Ekman numbers ${ {Ek}} \gtrsim 10^{-8}$), while realistic ${Ek}$ for geo- and astrophysical applications are orders of magnitude smaller. On the other hand, previously derived reduced equations of motion describing the leading-order behaviour in the limit of very rapid rotation ($ {Ek}\to 0$) cannot capture finite rotation effects, and the physically most relevant part of parameter space with small but finite ${Ek}$ has remained elusive. Here, we employ the rescaled rapidly rotating incompressible Navier–Stokes equations (RRRiNSE) – a reformulation of the Navier–Stokes–Boussinesq equations informed by the scalings valid for ${Ek}\to 0$, recently introduced by Julien et al. (2024) – to provide full DNS of RRRBC at unprecedented rotation strengths down to $ {Ek}=10^{-15}$ and below, revealing the disappearance of cyclone–anticyclone asymmetry at previously unattainable Ekman numbers (${Ek}\approx 10^{-9}$). We also identify an overshoot in the heat transport as ${Ek}$ is varied at fixed $\widetilde { {Ra}} \equiv {Ra}{Ek}^{4/3}$, where $Ra$ is the Rayleigh number, associated with dissipation due to ageostrophic motions in the boundary layers. The simulations validate theoretical predictions based on thermal boundary layer theory for RRRBC and show that the solutions of RRRiNSE agree with the reduced equations at very small ${Ek}$. These results represent a first foray into the vast, largely unexplored parameter space of very rapidly rotating convection rendered accessible by RRRiNSE.

We conducted a series of pore-scale numerical simulations on convective flow in porous media, with a fixed Schmidt number of 400 and a wide range of Rayleigh numbers. The porous media are modeled using regularly arranged square obstacles in a Rayleigh–Bénard (RB) system. As the Rayleigh number increases, the flow transitions from a Darcy-type regime to an RB-type regime, with the corresponding $Sh$–$Ra_D$ relationship shifting from sublinear scaling to the classical 0.3 scaling of RB convection. Here, $Sh$ and $Ra_D$ represent the Sherwood number and the Rayleigh–Darcy number, respectively. For different porosities, the transition begins at approximately $Ra_D = 4000$, at which point the characteristic horizontal scale of the flow field is comparable to the size of a single obstacle unit. When the thickness of the concentration boundary layer is less than approximately one-sixth of the pore spacing, the flow finally enters the RB regime. In the Darcy regime, the scaling exponent of $Sh$ and $Ra_D$ decreases as porosity increases. Based on the Grossman–Lohse theory (J. Fluid Mech. vol. 407, 2000, pp. 27–56; Phys. Rev. Lett. vol. 86, 2001, p. 3316), we provide an explanation for the scaling laws in each regime and highlight the significant impact of mechanical dispersion effects during the development of the plumes. Our findings provide some new insights into the validity range of the Darcy model.

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.