2. Simulation set-up and numerical methods

We perform numerical simulations of the flow arising in a fluid confined between two parallel no-slip vertical walls. The walls are separated by a horizontal distance $H$ and are held at fixed temperatures, with a temperature difference $\varDelta$ between the plates. As in the schematic in figure 1(b), we take the $x$-coordinate to be horizontal and parallel to the plates, the $y$-coordinate to be normal to the boundaries, and $z$ to be in the vertical direction. We consider a fluid satisfying the Oberbeck–Boussinesq approximation, such that changes in density are significant only in the buoyancy term, and are linear with respect to changes in temperature. We therefore treat the velocity field $\boldsymbol {u}=(u,v,w)$ as divergence-free ($\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {u} = 0$) and satisfying the Navier–Stokes momentum equation

(2.1)\begin{equation} \partial_t \boldsymbol{u} + (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla}) \boldsymbol{u} ={-}\rho^{{-}1}\,\boldsymbol{\nabla} p + \nu\,\nabla^2 \boldsymbol{u} + g \alpha \theta \boldsymbol{\hat{z}} + G \boldsymbol{\hat{x}}, \end{equation}

where $\rho$ is the mean fluid density (assumed constant), $p$ is the pressure, $\nu$ is the kinematic viscosity, $g$ is the acceleration due to gravity, and $\alpha$ is the thermal expansion coefficient. A time-dependent, spatially uniform forcing $G(t)$ is applied in the streamwise ($x$) direction to maintain a constant mean flow $\langle u \rangle = U$. The magnitude of this forcing is computed at every time step to exactly cancel out any variation in the mean flow. Previous work has shown that such a forcing produces results near-identical to those of a constant pressure gradient (Quadrio, Frohnapfel & Hasegawa Reference Quadrio, Frohnapfel and Hasegawa2016), but allows us to use the mean flow strength as an input parameter. The temperature field $\theta$ satisfies the advection–diffusion equation

(2.2)\begin{equation} \partial_t \theta + (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla}) \theta = \kappa\,\nabla^2 \theta, \end{equation}

where $\kappa$ is the thermal diffusivity. Periodic boundary conditions are applied to both $\boldsymbol {u}$ and $\theta$ in the $x$ and $z$ directions. Unless otherwise stated, the aspect ratio of the domain in these periodic directions is taken as $\varGamma =L/H=8$, such that the length of the domain $L$ in $x$ and $z$ is equal, and eight times the plate separation distance $H$.

Figure 1. Schematics of canonical mixed convection systems. (a) Mixed RB convection with Poiseuille-type forcing as studied by e.g. Domaradzki & Metcalfe (Reference Domaradzki and Metcalfe1988), Pirozzoli et al. (Reference Pirozzoli, Bernardini, Verzicco and Orlandi2017) and Yerragolam et al. (Reference Yerragolam, Howland, Stevens, Verzicco, Shishkina and Lohse2024). (b) Mixed VC with Poiseuille-type forcing, as used for the new simulations presented in this paper. Large arrows highlight the imposed horizontal pressure gradient, and smaller arrows highlight the large-scale flows driven by buoyancy in each configuration. The acceleration due to gravity is highlighted by a double arrow.

We perform direct numerical simulations of (2.1) and (2.2) using a highly parallelised flow solver that computes spatial derivatives using second-order central finite differences on a staggered grid configuration. The wall-normal diffusive terms are evolved in time using a Crank–Nicolson scheme, and all other terms are treated explicitly using a three-stage Runge–Kutta method. An adaptive time step is chosen using the constraint of a maximum Courant number of 1. The velocity is kept divergence-free to machine precision using a pressure correction step at each time step that is implemented with fast Fourier transforms in the periodic directions and a tridiagonal matrix solver in the wall-normal direction. A multiple-resolution technique is applied to evolve the velocity and temperature fields on independent grids, with cubic Hermite interpolation used for the buoyancy forcing and the advection of temperature. Detailed overviews of the numerical discretisation, the domain decomposition strategy and the multiple-resolution technique can be found in Verzicco & Orlandi (Reference Verzicco and Orlandi1996), van der Poel et al. (Reference van der Poel, Ostilla-Mónico, Donners and Verzicco2015) and Ostilla-Monico et al. (Reference Ostilla-Monico, Yang, van der Poel, Lohse and Verzicco2015), as well as in our software documentation.

The physical input parameters of the system are the Rayleigh number, the Prandtl number and the Reynolds number:

(2.3a–c)\begin{equation} \textit{Ra}= \frac{g\alpha\,\Delta H^3}{\nu \kappa},\quad \textit{Pr}= \frac{\nu}{\kappa},\quad \textit{Re}=\frac{U H}{\nu}. \end{equation}

When considering the strength of the flow driven by convection, it is often useful to consider the Grashof number $Gr$ as the relevant input parameter, and when comparing the relative strengths of buoyancy to pressure driving, we can construct a Richardson number as

(2.4a,b)\begin{equation} \textit{Gr} = \frac{\textit{Ra}}{\textit{Pr}} = \frac{g \alpha\,\Delta H^3}{\nu^2} = \left(\frac{U_f H}{\nu}\right)^2,\quad \textit{Ri} = \frac{\textit{Gr}}{\textit{Re}^2} = \frac{g \alpha\,\Delta H}{U^2} = \left(\frac{U_f}{U}\right)^2 . \end{equation}

These can both be considered as input parameters. Above, we have written $U_f = \sqrt {g\alpha \, \Delta H}$ as the free-fall velocity scale to give insight into the interpretation of these parameters.

In this study, we perform two sets of simulations to compare the relative impacts of the various input parameters. First, we fix $Gr=10^6$ and vary $1\leq Pr \leq 10$ along with $10^{2.5} \leq Re \leq 10^4$, which correspond to Richardson numbers $10^{-2} \leq Ri \leq 10$. For the second set, we fix $Pr=1$ and increase $Gr$ up to $10^8$ while again varying $Re$ up to $10^4$. A detailed overview of the parameters used for each simulation is given in table 1 of Appendix B. Each simulation is run to a statistically steady state, in which the flow statistics are computed and averaged for a minimum of 200 advective time units. For high $Ri\geq 1$, the relevant time unit is $H/U_f$, whereas for low $Ri\leq 1$, the relevant time unit is $H/U$. The wall-normal grid spacing is stretched following Pirozzoli & Orlandi (Reference Pirozzoli and Orlandi2021) to ensure sufficient resolution close to the wall, such that $\Delta y^+ < 0.1$ at the wall for the base velocity grid. In the periodic directions, the uniform grid spacing satisfies $\Delta x^+ \leq 5.4$ in every simulation. At the centre of the domain, the spacing of the refined grid satisfies $\Delta y_r \leq 1.05 l_B$, $\Delta x_r < 1.4 l_B$, where $l_B$ is the Batchelor scale computed using the plane-averaged turbulent kinetic energy (TKE) dissipation rate over the midplane. Full details of the grid sizes are provided in Appendix B.

3. Flow visualisation

We begin with a qualitative comparison of the simulations through visualisations of the temperature field and the local heat flux. Figure 2 displays the instantaneous local wall-normal heat flux at the boundary plate $y=0$ for cases with fixed $Pr=1$, and ranges $10^6\leq Gr\leq 10^8$, $10^3\leq Re\leq 10^4$. The dimensionless heat flux plotted here is defined as

(3.1)\begin{equation} q_\theta(x,z,t) ={-}\frac{H}{\varDelta} \left.\frac{\partial \theta}{\partial y}\right|_{y=0}, \end{equation}

such that its time- and plane-averaged value is equivalent to the Nusselt number, $Nu=\overline {q_\theta }$.

Figure 2. Vertical planes of the instantaneous local wall-normal heat flux at the boundary $y=0$. All simulations shown are with fixed $Pr=1$. The Richardson number $Ri=Gr/Re^2$ is comparable along diagonals from the upper left to the lower right, with the largest value ($Ri=100$) in the top right and the lowest value ($Ri=0.01$) in the bottom left. This figure is also available as an interactive JFM notebook: https://www.cambridge.org/S0022112024005986/JFM-Notebooks/files/Figure-2.ipynb.

As mentioned above, the relative strength of convection to the horizontal flow can be characterised by the Richardson number $Ri=Gr/Re^2$, which is constant along diagonals in figure 2. At high $Ri$, as in figure 2(c), the horizontal flow has little impact on the local distribution of the wall heat flux. The near-wall temperature structures are the same as in the case of no crossflow, with regions of large local heat flux (dark spots) separated by longer, streaky structures aligned in the vertical (Howland et al. Reference Howland, Ng, Verzicco and Lohse2022). As $Ri$ decreases, such as in figures 2(e,i) where $Ri=1$, the prominence of the large heat flux regions diminishes, and the streaks become aligned in the diagonal. This visualises the local direction of the flow along the wall, which at $Ri=1$ is due to a combination of the VC and the horizontal pressure gradient. At lower $Ri$, these structures become more aligned with the horizontal, eventually spanning the domain as in figure 2(g), which is reminiscent of classical low-speed streaks in turbulent channel flow (Kline et al. Reference Kline, Reynolds, Schraub and Runstadler1967; Antonia, Abe & Kawamura Reference Antonia, Abe and Kawamura2009). A more quantitative analysis of the change in the near-wall heat flux distribution will be provided in § 7.

In figure 3, we present visualisations of the same simulations but now at the midplane of the simulation domain $y=H/2$. As would be expected, the fields at higher $Gr$ and $Re$ exhibit structures with a wider range of spatial scales. Aside from this dynamical range, the effect of increasing $Re$ is less noticeable at the midplane than at the wall in figure 2. At the centre of a turbulent channel, the mean profile of horizontal velocity is relatively flat, with zero mean shear and a local minimum in TKE. By contrast in VC, a mean shear in the vertical velocity drives the generation of turbulence in the bulk, argued by Li et al. (Reference Li, Jia, Liu, Jiao and Zhang2023) to follow a mixing-layer-like behaviour. In figure 3, greater mixing at higher $Gr$ leads to smaller values of the temperature perturbations in the midplane, although the same effect is not evident as $Re$ increases. Compared to the mixed RB system, where gravity is orthogonal to the plates and the temperature field organises into large, streamwise-aligned coherent rolls, the fields in mixed VC appear rather featureless. A hint of such large-scale rolls is noticeable only for the cases dominated by strong pressure driving with high $Re$ and low $Ri$ (figures 3g,h), where the temperature structures appear more aligned along the streamwise ($x$) axis. This contrast in the organisation of mixed convection systems will be investigated more quantitatively in § 6 through spectral analysis, but we first turn to the global responses of the two mixed convection systems in the next section.

Figure 3. Vertical planes of the temperature field at the channel centre $y=H/2$. See the caption of figure 2 for details of the simulations presented. Here, $\theta _0$ is the arbitrary reference temperature that is the midpoint of the two boundary values. This figure is also available as an interactive JFM notebook: https://www.cambridge.org/S0022112024005986/JFM-Notebooks/files/Figure-3.ipynb.

4. Global response quantities in mixed convection systems

In this section, we compare the responses of the mixed VC set-up with those of mixed RB using the simulation data of Yerragolam et al. (Reference Yerragolam, Howland, Stevens, Verzicco, Shishkina and Lohse2024). Those simulations cover comparable ranges of parameters $10^6\leq Gr \leq 10^8$, $0.5\leq Pr \leq 5$ and $Re\leq 10^4$ in a domain of streamwise aspect ratio $\varGamma _x=8$ and spanwise aspect ratio $\varGamma _y=4$. The flow solver shares an identical code base except for the multiple-resolution technique, which is used only in the newly reported mixed VC simulations.

4.1. Friction coefficients

A key global response parameter in shear-driven flows is the friction coefficient

(4.1)\begin{equation} C_f = \frac{2 \tau}{\rho U^2}, \end{equation}

where $\tau$ is the wall shear stress and $U$ is the velocity magnitude in the bulk. Since the friction coefficient is determined solely by the velocity field, we expect the dependence on the Prandtl number to be relatively weak, and consider a relationship $C_f(Re)$. Power-law scalings for $C_f(Re)$ can be derived for laminar boundary layer flows, with e.g. $C_f\sim Re^{-1}$ applicable to Couette or Poiseuille flows, and $C_f \sim Re^{-1/2}$ arising from the classical Blasius boundary layer solution (Schlichting & Gersten Reference Schlichting and Gersten2016). For a turbulent boundary layer in the sense of Prandtl and von Kármán, the friction coefficient satisfies the relation

(4.2)\begin{equation} \sqrt{\frac{2}{C_f}} = \frac{1}{\kappa_u} \log{\left( \textit{Re}\,\sqrt{\frac{C_f}{8}} \right) } + B, \end{equation}

known as the Prandtl friction law after Prandtl (Reference Prandtl1932). The von Kármán constant $\kappa _u$ typically takes a value of approximately 0.4, and the intercept $B$ is close to 4, but the exact values, their universality and the way in which they are fit to data remain an active topic of research (Monkewitz & Nagib Reference Monkewitz and Nagib2023). Due to our similar set-up and numerical methods, we take the values suggested by Pirozzoli, Bernardini & Orlandi (Reference Pirozzoli, Bernardini and Orlandi2014), namely $\kappa _u = 0.41$ and $B=5$.

Since the mixed VC flow is driven in orthogonal directions by the pressure gradient and by buoyancy, we can construct separate friction coefficients for each component of the wall shear stress. Understanding the response of the friction coefficients in this context relies on choosing an appropriate Reynolds number for each component of the flow. For the streamwise ($x$) direction in which the mean flow is imposed by a pressure gradient, this Reynolds number is simply the input parameter defined in (2.3a–c).

We consider the response of the streamwise friction coefficient in figure 4, where only the streamwise component of the shear stress ${\tau =\rho \nu \,\partial _y \bar {u}}$ is applied to the definition (4.1). The global response of mixed VC is near identical to that of mixed RB, with a transition from a laminar power-law scaling to the Prandtl friction law of (4.2). For comparable parameter values, the largest difference in $C_f$ between the mixed RB and mixed VC cases is 16 %. In the laminar scaling regime, stronger buoyancy driving (characterised by higher $Gr$) leads to an increase in the streamwise skin friction for a given $Re$. As suggested earlier, the dependence of $C_f$ on $Pr$ is very weak compared to the other control parameters. Unlike in standard Poiseuille or Couette flow, where a subcritical transition arises due to instability of the laminar base flow and leads to a jump in $C_f$, the streamwise boundary layer transition in mixed convection systems appears smooth. We anticipate that the laminar scaling regime remains relevant until its intersection with the friction law (4.2). The increase of $C_f$ with $Gr$ in the laminar regime would therefore delay the transition to this ‘fully turbulent’ Prandtl friction law (4.2) to higher $Re$. At low $Re$, although the relationship exhibits a laminar-like scaling, one should recall that the convection flow in the interior remains turbulent. In figure 4(b), we focus on clarifying this low $Re$ regime and the increase in $C_f$ with stronger convection. Across both mixed convection systems and for a range of $Pr$, we find a collapse of the data upon rescaling by $Gr^{1/4}$. The $Re^{-1}$ scaling that arises from laminar profiles in Couette/Poiseuille flow appears somewhat too steep to describe the data accurately. Blass et al. (Reference Blass, Zhu, Verzicco, Lohse and Stevens2020) reported a scaling $C_f\sim Re^{-0.90}$ in C-RB, but at this time there is no theoretical basis for such a scaling. Note that one could equivalently express the simplified $Re^{-1}$ collapse as $C_f \sim Ri^{1/4}\,Re^{-1/2}$ using the definitions of (2.3a–c). In the case of mixed convection, the buoyancy-driven flow generates non-zero Reynolds stresses in the equation for the mean profile, which can be expected to modify the mean momentum budget and lead to such an increase in $C_f$. This will be analysed in further detail for mixed VC in § 5.

Figure 4. (a) Friction coefficients calculated for the streamwise component of velocity as a function of Reynolds number. (b) Friction coefficients normalised by $Gr^{1/4}$ to collapse the data at low $Re$. Large markers are used for the mixed VC simulations, whereas small dots show the data from the mixed RB cases of Yerragolam et al. (Reference Yerragolam, Howland, Stevens, Verzicco, Shishkina and Lohse2024). Grashof numbers are highlighted by the colour of the markers, and for the mixed VC cases, different markers signify different $Pr$. Black dashed lines show a scaling law of $Re^{-1}$, and the black solid line marks the Prandtl friction law of (4.2). This figure is also available as an interactive JFM notebook: https://www.cambridge.org/S0022112024005986/JFM-Notebooks/files/Figure-4.ipynb.

We now turn to the friction coefficient associated with the buoyancy-driven component of the flow. For the mixed VC system, we can simply take the peak velocity $W_{max}$ of the mean vertical velocity $\bar {w}(y)$ as the relevant velocity scale, and directly measure the mean vertical shear stress at the wall, $\tau = \rho \nu \,\partial _y \bar {w}$. In the RB configuration, the convection has no preferential direction along the walls, resulting in zero mean shear stress. However, we can still construct a friction coefficient associated with the persistent large-scale circulation by using the root mean square (r.m.s.) horizontal velocity profile $u_H(y) = (\overline {u^2 + v^2})^{1/2}$. When a horizontal crossflow is added to the RB system, the large-scale circulation aligns itself perpendicular to the imposed flow (Pirozzoli et al. Reference Pirozzoli, Bernardini, Verzicco and Orlandi2017; Yerragolam et al. Reference Yerragolam, Verzicco, Lohse and Stevens2022), so for the mixed RB cases, we construct a friction coefficient using only the spanwise r.m.s. velocity. Defining the friction coefficient in this way will be appropriate only for cases where the convectively driven flow is stronger than the spanwise velocity fluctuations induced by the turbulent shear flow, i.e. for $Ri\geq O(1)$.

In terms of the Reynolds number, the plate separation $H$ is no longer the appropriate length scale for describing the boundary layer dynamics of the convectively driven flow. As shown in the inset of figure 5(a), the mean profile $\bar {w}(y)$ of the vertical velocity in (mixed) VC reaches its peak value at a certain wall-normal distance $\delta$. From this, we can define a boundary layer Reynolds number

(4.3)\begin{equation} \textit{Re}_\delta = \frac{W_{max} \delta}{\nu} \end{equation}

that drives the behaviour at the wall. We construct an analogous Reynolds number for the mixed RB system using the spanwise r.m.s. velocity profile. The friction coefficients of the convective flow component are plotted against $Re_\delta$ in figure 5. Note that $Re_\delta$ is not known a priori, but is itself a response parameter of the system that varies with $Gr$, $Pr$ and $Re$. Similar to the low-$Re$ regime for the streamwise friction, we observe a power-law scaling close to $C_f\sim Re_\delta ^{-1}$.

Figure 5. (a) Friction coefficients associated with the convectively driven flow as a function of the boundary layer Reynolds number $Re_\delta =W_{max} \delta /\nu$, where $W_{max}$ and $\delta$ are computed as in the mean velocity profile shown in the inset. (b) The same data normalised against the values from the corresponding natural convection system. For the mixed VC cases, different symbols denote different $Pr$ values as outlined in the legend of figure 4. The colour of the markers represents $Gr$, and simulations with $Ri<1/4$ are plotted as semi-transparent. This figure is also available as an interactive JFM notebook: https://www.cambridge.org/S0022112024005986/JFM-Notebooks/files/Figure-5.ipynb.

This is made clearer in figure 5(b), where we collapse the data using the friction coefficient $C_0$ and Reynolds number $Re_{\delta,0}$ obtained from the corresponding natural convection flows, matching $Ra$ and $Pr$. Deviations from the scaling relation are observed for cases where $Ri<1/4$, highlighted by translucent symbols in the figure. As mentioned above, for the mixed RB system this is likely an artefact of using the spanwise r.m.s. velocity to construct the friction coefficient. However, we also observe the discrepancy for low $Ri$ in mixed VC, suggesting that at low Richardson numbers, the turbulence generated by the imposed horizontal flow disrupts the near-wall vertical velocity. Within the range of parameters explored here, the pressure-driven horizontal flow does not modify the vertical $Re_\delta$ by more than a factor 4, suggesting that even in the case of mixed convection, the primary control parameters determining $Re_\delta$ are $Gr$ and $Pr$. As discussed in the Appendix of Howland et al. (Reference Howland, Ng, Verzicco and Lohse2022), a ‘fully turbulent’ transition of this boundary layer may be possible, analogous to the so-called ‘ultimate regime’ in RB convection (Lohse & Shishkina Reference Lohse and Shishkina2023), but only at very high $Gr$. A more in-depth analysis of the mean vertical momentum budget in mixed VC will be presented in § 5.

4.2. Nusselt number

The dimensionless heat flux through the system is characterised by the Nusselt number, defined as

(4.4a,b)\begin{equation} \textit{Nu} = \frac{F_\theta}{\kappa \varDelta / H},\quad F_\theta ={-}\kappa\,\frac{\partial \bar{\theta}}{\partial y} + \overline{v'\theta'} . \end{equation}

Here, $F_\theta$ is the horizontal heat flux through the system (normalised by the specific heat $\rho c_p$), and the overbar denotes averaging in time and in directions parallel to the plates. Integration of the mean temperature equation (3.1) shows that $F_\theta$ is constant across the domain in a statistically steady state.

In figure 6(a), we plot the Nusselt number compensated by $Pr^{1/2}$. This prefactor seems appropriate for the $Pr$ dependence of passive scalar transport in turbulent boundary layers when $Pr\lesssim O(1)$ (Kays et al. Reference Kays, Crawford, Weigand and Kays2005), although at higher $Pr$, one expects a transition towards a $Pr^{1/3}$ dependence (Kader & Yaglom Reference Kader and Yaglom1972; Alcántara-Ávila & Hoyas Reference Alcántara-Ávila and Hoyas2021). The data of both systems converge towards the expression

(4.5)\begin{equation} \textit{Nu}\approx \frac{C_f(\textit{Re})}{4}\,\textit{Re}\,\textit{Pr}^{1/2}, \end{equation}

where $C_f(\textit {Re})$ follows the Prandtl law of (4.2). This expression draws a parallel between the transport of heat and momentum at the wall, known as the Reynolds analogy, and describes the heat transport in ‘forced convection’ when buoyancy no longer affects the flow. From these data, we anticipate that the forced convection expression (4.5) applies for Reynolds numbers greater than that where the friction coefficient begins to follow the turbulent friction law shown in figure 4. At low $Re$ (or more precisely high $Ri$), the Nusselt number responses of the two systems (VC and RB) do not match as precisely as the friction coefficients. Indeed, in the absence of an external flow, VC and RB do not exhibit the same $Nu(Ra,Pr)$ response due to the lack of coupling between the kinetic energy budget and the heat flux in VC (Ng et al. Reference Ng, Ooi, Lohse and Chung2015).

Figure 6. Nusselt numbers plotted as functions of Reynolds number. (a) The Nusselt number is normalised by $Pr^{1/2}$ to focus on the response in the shear-dominated regime at high $Re$. The black dashed line marks the Reynolds analogy $Nu\approx \frac {1}{4}C_f\,Re\,Pr^{1/2}$, where $C_f$ satisfies the Prandtl friction law (4.2). (b) The data are normalised by the values associated with natural convection. Here, the black dashed line marks the recently proposed scaling relation $Nu/Nu_0 \sim (\sqrt {1 + (Re/Re_0)^2})^{-1/5}$ from Yerragolam et al. (Reference Yerragolam, Howland, Stevens, Verzicco, Shishkina and Lohse2024). This figure is also available as an interactive JFM notebook: https://www.cambridge.org/S0022112024005986/JFM-Notebooks/files/Figure-6.ipynb.

However, we observe a more universal behaviour when the Nusselt numbers in the mixed convection systems are normalised by the values for the equivalent natural convection systems $Nu_0(Gr, Pr, Re) = Nu(Gr, Pr)|_{Re=0}$. In figure 6(b), the normalised Nusselt numbers are plotted as a function of the input Reynolds number normalised by the Reynolds number of the natural convection case $Re_0=W_0H/\nu$. For the mixed VC cases, we define $W_0(Gr, Pr, Re) = W_{max}(Gr, Pr)|_{Re=0}$ as the peak velocity of the natural VC flow, as highlighted in the inset of figure 5(a). For mixed RB, we follow Yerragolam et al. (Reference Yerragolam, Howland, Stevens, Verzicco, Shishkina and Lohse2024) in using the volume-averaged r.m.s. velocity for $W_0$, and note that the results are insensitive to this choice of velocity scale in describing the ‘wind’ of the large-scale circulation. For $Re/Re_0=O(1)$, all the data from both configurations collapse onto a single curve, showing a drop in the heat flux of up to 25 %. Given this collapse, the $Re_0$ of the natural convection appears to be a critical $Re$ above which the Nusselt number is significantly affected by the horizontal crossflow. Yerragolam et al. (Reference Yerragolam, Howland, Stevens, Verzicco, Shishkina and Lohse2024) provide an estimate for the drop in $Nu$, derived from the kinetic energy balance in mixed RB, but this balance cannot be related to the horizontal heat flux that is relevant for mixed VC.

In summary, in this section we have demonstrated the universality in the global response parameters of mixed RB and mixed VC, namely in the friction coefficient $C_f$ and the Nusselt number, and the limitations of this universality. In the following sections, we will compare more local quantities, starting with the wall-normal profiles.

5. Wall-normal profiles in mixed VC

We now turn to the first- and second-order statistics, averaged parallel to the plates, to further investigate the dynamics behind the observed global responses. For clarity, we focus solely on the new simulations of mixed VC, and study the variation across the three-parameter space of $Gr$, $Pr$ and $Re$.

We begin with the response of the mean streamwise velocity $\bar {u}(y)$ in figure 7. For a fixed $Gr=10^6$, as in figures 7(a,b), the effect of increasing $Re$ can be seen most clearly when the mean velocity profiles are scaled by viscous wall units in figure 7(b). As $Re$ increases, the velocity profile tends towards the classical log-law profile, with the case $Re=10^4$ closely matching the profile of turbulent Poiseuille flow, as in e.g. Lee & Moser (Reference Lee and Moser2015). The effect of stronger thermal convection on the mean profile is also similar to that observed in other mixed convection systems in the literature (Scagliarini et al. Reference Scagliarini, Einarsson, Gylfason and Toschi2015; Blass et al. Reference Blass, Zhu, Verzicco, Lohse and Stevens2020). In figure 7(c), where $Re=10^3$ is fixed and $Gr$ varies between $10^6$ and $10^8$, higher $Gr$ leads to a flatter mean profile in the bulk of the channel. This is illustrated further in wall units in figure 7(d), where a significant drop in $\bar {u}^+$ is observed for $y^+ =O(10)$. Plus symbols denote scaling in viscous wall units, with velocity $u_\tau = \sqrt {\tau _u/\rho }$ and length $\nu /u_\tau$. Such a drop is consistent with the previous findings of Scagliarini et al. (Reference Scagliarini, Einarsson, Gylfason and Toschi2015) for mixed RB, who proposed a modified log law based on mixing length arguments coupled to the temperature field.

Figure 7. Mean profiles of the streamwise velocity $\bar {u}(y)$. (a,b) Data from simulations at a fixed $Gr=10^6$, varying $Re$ and $Pr$. The $Re$ variation is denoted by the line colour, and the $Pr$ variation is denoted by the line style. (c,d) Present data from simulations at a fixed $Re=10^3$, varying $Gr$ and $Pr$. (a,c) The profiles are normalised by the imposed bulk velocity $\langle u\rangle = U$ and the plate separation $H$. (b,d) The profiles are presented in viscous wall units, where $\bar {u}^+ = \bar {u}/{u_\tau }$ and $y^+=y u_\tau /\nu$. Dashed black lines represent the linear relation $u^+=y^+$ and the logarithmic region $u^+ = \kappa _u^{-1} \log {y^+} + B$, where the von Kármán constant $\kappa _u = 0.41$ and $B=5$ are taken from Pirozzoli et al. (Reference Pirozzoli, Bernardini and Orlandi2014). This figure is also available as an interactive JFM notebook: https://www.cambridge.org/S0022112024005986/JFM-Notebooks/files/Figure-7.ipynb.

Further insight for the streamwise velocity can be gained from the appropriate component of the Reynolds stress. Considering a statistically steady state, we average the streamwise component of the momentum equation (2.1) to obtain

(5.1)\begin{equation} \partial_y \overline{uv} = \nu\,\partial_{yy}\bar{u} + \bar{G}, \end{equation}

where an overbar denotes an average in the periodic ($x, z$) directions and in time. Note that wall-normal advective fluxes in incompressible channel flows are purely turbulent, i.e. $\overline {uv}\equiv \overline {u'v'}$ since $\bar {v}\equiv 0$. From volume averaging, we can also relate the mean pressure gradient forcing to the mean wall shear stress through $\bar {G}=2\tau _u/\rho H$. The first integral of (5.1) can therefore be written as

(5.2)\begin{equation} \partial_{y^+} \bar{u}^+- \overline{uv}^+= 1 - \frac{2y}{H} . \end{equation}

From (5.1) and (5.2), the close coupling of the mean streamwise velocity and the Reynolds stress $\overline {uv}$ is evident. We therefore present the profiles of $\overline {uv}(y)$ in figure 8. As expected, the Reynolds stress dominates the viscous contribution to (5.2) away from the walls, leading to a balance of $-\overline {uv}^+\approx 1 - 2y/H$ as shown in figures 8(a,c). Relative to $H$, the boundary layer in which the viscous term is relevant becomes thinner as both $Re$ and $Gr$ increase. The near-wall behaviour of $\overline {uv}$ exhibits a remarkable collapse when scaled by $Ri^{1/4}$ as in figures 8(b,d), except for the highest $Re$ case with $Ri=0.01$.

Figure 8. Wall-normal profiles of the streamwise Reynolds stress component $\overline {uv}(y)$ scaled by the streamwise friction velocity $\overline {uv}^+ = \overline {uv}/u_\tau ^2$. As in figure 7, (a,b) are at fixed $Gr=10^6$ and (c,d) are at fixed $Re=10^3$. Colours and line styles are as detailed in the caption of figure 7. (a,c) Profiles relative to the plate separation $H$; (b,d) profiles in viscous wall units scaled by $Ri^{1/4}$. The inset of (d) presents the data of (b,d) with a logarithmic scale on both axes to highlight the near-wall collapse of the data. This figure is also available as an interactive JFM notebook: https://www.cambridge.org/S0022112024005986/JFM-Notebooks/files/Figure-8.ipynb.

The additional factor $Ri^{1/4}$ suggests that the appropriate near-wall length scale for the Reynolds stress is modified from the standard viscous wall unit as

(5.3a,b)\begin{equation} \textit{Ri}^{1/4} y^+= \frac{y}{l},\quad l =\frac{\nu}{\sqrt{\dfrac{U_f}{U}\,\dfrac{\tau_u}{\rho}}}. \end{equation}

The additional prefactor $U_f/U$ in front of the shear stress suggests that the vertical, convectively driven component of shear cannot be neglected when considering the streamwise Reynolds stress. An improved collapse to that seen in figures 8(c,d) can be found by computing a viscous length scale using the total shear stress at the wall, $\tau = \sqrt {\tau _u^2 + \tau _w^2}$. As shown explicitly in Appendix A, with this scaling, the Reynolds stress for the highest $Re$ case also matches the other curves. The presence of the convectively driven flow increases the near-wall Reynolds stress, which in turn leads to the flattened mean velocity profiles observed in figure 7(d). This result explains qualitatively the origin of the change in $C_f$ with $Gr$ seen in figure 4, where enhanced buoyancy driving led to a larger skin friction. The increase in near-wall Reynolds stress that arises due to convection thins the boundary layer of the mean horizontal velocity, which in turn produces a larger mean gradient at the wall and a larger friction coefficient $C_f$.

The effect of convection on the shear is by no means a one-way interaction. This is evident from the modification of the mean vertical velocity profile by the imposed horizontal flow. In figures 9(a,b), vertical velocity profiles are shown for fixed $Gr=10^6$ and varying $Re$, with the reference natural convection case ($Re=0$) highlighted in blue for comparison. Compared to the natural VC case, the introduction of horizontal driving at moderate $Re$ leads to an increase in the peak vertical velocity, and hence an increase in the mean shear both in the bulk and at the walls. At the highest $Re=10^4$ (the darkest red line in figure 9), a subsequent decrease is observed in the peak vertical velocity, as well as a nonlinear profile in the bulk. None of the cases studied here exhibits a log layer in the vertical velocity, with the largest $\bar {w}^+$ being approximately 6.5 for the most strongly convective case, $Gr=10^8$. For fixed $Ri=1$ (shown in figures 9c,d), all cases show a similar increase in the peak velocity, and the mean gradient in the bulk appears largely independent of $Gr$ and $Pr$. The distance of the velocity peak from the wall (compared to the channel width $H$) decreases for larger $Gr$, but the value of the peak velocity in free-fall units does not depend strongly on $Gr$. This similarity at constant $Ri$ is suggestive that the vertical velocity modification is determined primarily by the relative strength of the horizontal flow to convection.

Figure 9. Mean profiles of the vertical velocity $\bar {w}(y)$. As in figure 7, (a,b) are at fixed $Gr=10^6$, but here (c,d) are at fixed $Ri=1$. (a,c) Profiles are normalised by the free-fall velocity $U_f = \sqrt {g\alpha \,\Delta H}$ and the plate separation $H$. (b,d) Profiles are presented in viscous wall units, where $\bar {w}^+ = \bar {w}/{w_\tau }$ and $y^+=y w_\tau /\nu$. Reference data for natural VC (with $Gr=10^6$, $Re=0$) are shown in blue. This figure is also available as an interactive JFM notebook: https://www.cambridge.org/S0022112024005986/JFM-Notebooks/files/Figure-9.ipynb.

To further investigate the behaviour of the vertical velocity profiles, we now turn to the mean vertical momentum budget. Unlike the streamwise velocity in (5.1)–(5.2), the vertical velocity is not only tied to the Reynolds stress profile. Rather, the mean vertical momentum equation reads

(5.4)\begin{equation} \partial_y \overline{vw} = \nu\,\partial_{yy} \bar{w} + g\alpha \bar{\theta} . \end{equation}

Close to the wall, we expect the Reynolds stress to be negligible and a balance to arise between buoyancy and viscosity. Due to the symmetry of the boundary conditions, all three terms must be zero at the channel centre. In the bulk of the flow, the mean velocity is approximately linear, so we expect a balance between buoyancy and Reynolds stress. These features are present in each of the simulations highlighted in figure 10. As $Re$ increases, the key modification to the budget arises in the Reynolds stress term $-\partial _y \overline {vw}$. In natural VC, there is a minuscule positive peak in the Reynolds stress term close to the wall, but its amplitude is so small that it is indistinguishable in figures 10(a,f). This peak grows with $Re$, becoming visible at $Re=10^{3.5}$ in figures 10(d,i), coinciding with a flattened profile of the viscous term (shown in green). By $Re=10^4$, at which the Reynolds stresses are significantly energised by the horizontal forcing, the nonlinear term peak exceeds the contribution from the buoyancy term. This leads to a significant drop in the viscous term at $y\approx 10^{-2}H$, which determines the modified mean velocity observed in figures 9(a,b).

Figure 10. Profiles of the mean vertical momentum budget terms for $Pr=1$ and (a–e) $Gr=10^6$ or ( f–j) $Gr=10^7$. Reynolds numbers vary from (a,f) $Re=0$ to (e,j) $Re=10^4$. Colours denote the budget term being plotted as detailed in the legend, and each budget term is plotted normalised against the buoyancy scale $g\alpha \varDelta$. Within each column, semi-transparent lines are added to show the profiles for the Reynolds numbers highlighted by the other rows as comparison. The wall-normal coordinate $y$ is plotted on a logarithmic axis to highlight the variation in the near-wall region. This figure is also available as an interactive JFM notebook: https://www.cambridge.org/S0022112024005986/JFM-Notebooks/files/Figure-10.ipynb.

Changes in the mean temperature profile $\bar {\theta }$ in figure 10 are more subtle, with the most obvious feature being a slight drop in figure 10(e), coinciding with the Reynolds stress peak. A clearer picture of the mean temperature response can be found in figure 11, where profiles are plotted both on linear axes and in wall units. Since the dimensionless wall temperature gradient is equivalent to the Nusselt number, which is most strongly dependent on $Ra$, we compare the temperature profiles in figures 11(a,b) at varying $Re$ and $Pr$ but fixed $Ra=10^7$. As shown in these plots, the temperature gradient in the bulk increases with $Re$, and any possible log-law profile does not collapse to a universal slope or coefficient. For fixed $Ri=1$, however, shown in figures 11(c,d), the bulk gradient appears independent of $Gr$, suggesting that $Ri$ and $Pr$ are the key control parameters determining the flow properties away from the walls.

Figure 11. Mean temperature profiles for (a,b) fixed $Ra=10^7$ and (c,d) fixed $Ri=1$. As in figures 7–9, colours denote changes in $Re$ and $Gr$, whereas line styles show changes in $Pr$. (a,c) Temperature is normalised by the difference across the plates $\varDelta$. (b,d) Temperature is shown in terms of wall units $\theta ^+=(\theta _w - \theta ) U_\tau / F_\theta$, where $\theta _w$ is the temperature at the wall $y=0$. The $U_\tau$ value used here and in $y^+=yU_\tau /\nu$ is calculated using both components of the wall shear stress, as defined in (A1). This figure is also available as an interactive JFM notebook: https://www.cambridge.org/S0022112024005986/JFM-Notebooks/files/Figure-11.ipynb.

The mean profiles presented here do not further clarify the origin of the non-monotonic behaviour of $Nu$ with $Re$ seen in figure 6. When normalised by the wall heat flux, the mean temperature profiles in figure 11(b) show a monotonic increasing trend with $Re$ away from the wall. The non-monotonic behaviour observed in the vertical velocity profiles of figures 9(a,b) suggests an intrinsic connection between the vertical flow and the wall-normal heat transport, although the very subtle changes in the vertical momentum budget remain difficult to interpret in the context of the non-monotonic Nusselt number variation. The later sections on spectral analysis and statistics of the boundary layer aim to shed more light on this global heat transport.

To conclude our analysis of the mean profiles, we compare the new simulations of mixed VC to the universal functions of the MOST (Monin & Obukhov Reference Monin and Obukhov1954; Foken Reference Foken2006). As outlined in the Introduction, the theory considers a region of the flow where viscosity and molecular diffusion are negligible, and the appropriate length scale is the Obukhov length

(5.5)\begin{equation} L_O=\frac{u_\tau^3}{g\alpha F_\theta}. \end{equation}

On dimensional grounds, universal functions that depend only on $y/L_O$ for the first- and second-order statistics can then be constructed for the temperature and velocity fields. One key assumption used to derive the universal functions for velocity is that the turbulent momentum flux $\rho \,\overline {uv}$ is constant over the region of interest. However, as shown in figure 8, this is not true in a turbulent channel flow driven by a pressure gradient, where the Reynolds stress follows a linear relation in the bulk. We can therefore expect the universal functions of Monin & Obukhov (Reference Monin and Obukhov1954) to hold for only a small range of $y$ in our simulations.

Considering these caveats, we present a collection of mean profiles from the mixed VC simulations as a function of $y/L_O$ in figure 12. Following Pirozzoli et al. (Reference Pirozzoli, Bernardini, Verzicco and Orlandi2017), who present a similar figure (their figure 16) for mixed RB convection, we emphasise only the region of the domain where the turbulent fluxes of heat and momentum are greater than 80 % of their maximum value, with the rest of the profiles made semi-transparent. We compare our results with the updated scaling theories of Kader & Yaglom (Reference Kader and Yaglom1990), and for the mean temperature and velocity profiles, the commonly used Businger–Dyer relations outlined in Foken (Reference Foken2006). For the highlighted region where the fluxes are close to their peak, the majority of the results collapse onto a single function of $y/L_O$ close to these profiles. In figure 12(a), there are two notable discrepancies from the theoretical estimates for the mean temperature profile. First, the values at transitional $y/L_O\approx 10^{-1}$ are approximately a factor 2 larger than the suggested profiles. Second, the effective power-law scaling observed for strong buoyancy driving $y\geq L_O$ appears steeper than the $-1/3$ proposed by Kader & Yaglom (Reference Kader and Yaglom1990) (dashed black lines), and closer to the $-1/2$ scaling associated with the asymptotic limit of the Businger–Dyer relation (solid blue line). A discrepancy in the two theoretical predictions for the mean velocity profile at $y\gg L_O$, highlighted in figure 12(b), cannot be resolved from the results presented here, although we note that the data here agree well with the mixed RB data of Pirozzoli et al. (Reference Pirozzoli, Bernardini, Verzicco and Orlandi2017). The variation in $Pr$ across the simulations collapses for most quantities in figure 12, except for the r.m.s. temperature shown in figure 12(e). The origin of this discrepancy with $Pr$ is currently unclear, but understanding how the $Pr$ dependence may affect the Monin–Obukhov profiles is important in the context of their applications outside the atmosphere, for example in the surface layer of the upper ocean (Zheng, Harcourt & D'Asaro Reference Zheng, Harcourt and D'Asaro2021).

Figure 12. Comparison of the mixed VC simulation results with the ‘universal functions’ of the MOST for first- and second-order statistics. Profiles are shown as semi-transparent unless the advective wall-normal fluxes of heat $\overline {v \theta }$ and momentum $\overline {vu}$ are within 80 % of their maximum value. Consistent with the Monin–Obukhov formulation, all quantities here are normalised by the horizontal component of the friction velocity $u_\tau$ and the wall-normal heat flux $F_\theta$. This figure is also available as an interactive JFM notebook: https://www.cambridge.org/S0022112024005986/JFM-Notebooks/files/Figure-12.ipynb.

6. Spectral analysis

We now investigate the scale dependence of the thermal structures and heat flux in mixed convection through analysis of the power spectra. To ensure that we capture the full range of dynamical scales, we perform further simulations in extended aspect ratio domains, with $L_x=L_z=24 H$. The details of these simulations are outlined in table 2 of Appendix B.

In unsheared RB convection, Krug, Lohse & Stevens (Reference Krug, Lohse and Stevens2020) show that large-scale patterns, also known as superstructures, can be identified from a low-wavenumber peak in the power spectrum of the temperature field and in the co-spectrum of the wall-normal heat flux. This peak denotes the scale of the large-scale circulation or ‘wind’ of convection that is key to theories describing RB convection. Since there is no preferential horizontal flow direction in RB convection, one-dimensional spectra can be analysed, but in mixed VC, the two directions parallel to the plates must be considered separately.

We therefore compute time-averaged spectra using two-dimensional Fourier transforms in the periodic directions. For visualisation purposes, we present separate one-dimensional spectra for the streamwise ($x$) and spanwise ($z$) wavenumbers, which are computed by integrating over the other wavenumber. Precisely, the two- and one-dimensional co-spectra of any two variables $f$ and $g$ are defined as

(6.1a,b)\begin{equation} \varPhi_{fg}(k_x, y, k_z) = {\rm Re}[\,\hat{f}^* \hat{g}],\quad \varPhi_{fg}(k_x,y) = \int {\rm Re}[\,\hat{f}^* \hat{g}]\, \mathrm{d}k_z, \end{equation}

where $\hat {f}$ denotes the two-dimensional Fourier transform of $f$ in $x$ and $z$, and ${\rm Re}$ denotes the real part. The Fourier transforms are normalised such that integrating the spectrum over wavenumber space recovers the corresponding volume-averaged quantity $\langle\, f g \rangle = \iint \varPhi _{fg} \,\mathrm {d}k_x \,\mathrm {d}k_z$.

The power spectrum of temperature $\varPhi _{\theta \theta }$ is presented in figure 13 for the four extended simulations. In figure 13(a), for the standard RB configuration, we see behaviour similar to that in Krug et al. (Reference Krug, Lohse and Stevens2020), with a distinct, sharp low-wavenumber peak at $k H\approx 1$, and a more broad peak at smaller scales. As mentioned above, the lack of a preferential direction means that the two directional spectra are virtually identical. When shear is added to the RB system in figure 13(b), the horizontal isotropy in the system is destroyed. As expected from previous work on mixed RB (Pirozzoli et al. Reference Pirozzoli, Bernardini, Verzicco and Orlandi2017; Blass et al. Reference Blass, Zhu, Verzicco, Lohse and Stevens2020,  Reference Blass, Tabak, Verzicco, Stevens and Lohse2021; Yerragolam et al. Reference Yerragolam, Verzicco, Lohse and Stevens2022), coherent rolls aligned with the streamwise axis dominate the signal. Note the difference in $y$ axis between figures 13(a) and 13(b). In the streamwise spectrum $\varPhi _{\theta \theta }(k_z,y)$, a sharp peak at $kH\approx 2$ is visible at all wall-normal locations, but is particularly prominent in the near-wall region (light green) where there is a maximum in $\overline {\theta ^2}(y)$. The corresponding wavelength of these structures is $\lambda = 2{\rm \pi} H/k_z \approx 3 H$, which is consistent with that observed previously in the mixed convection literature at $Ri=1$ (e.g. Pirozzoli et al. Reference Pirozzoli, Bernardini, Verzicco and Orlandi2017). The alignment of the flow structures in the streamwise direction leads to a broad contribution to the streamwise spectrum at low $k_x$. This behaviour is purely a result of the regular alignment, and is not related to any domain size effects.

Figure 13. One-dimensional power spectra $\varPhi _{\theta \theta }(k,y)$ of temperature for the extended domain simulations described in table 2. Solid lines denote spectra as functions of $k=k_x$, and dashed lines are functions of $k=k_z$. The colour of each line is determined by its wall normal location $y$, with the specific locations highlighted on the colour bar. Each spectrum is multiplied by the wavenumber $k$ such that area under each curve is representative of the relative contribution in wavenumber space. (a,c) Simulations of natural convection ($Re=0$). (b,d) Simulations of mixed convection at $Ri=1$. (a,b) The RB cases with gravity in the wall-normal ($y$) direction. (c,d) The VC cases with gravity parallel to the wall in z. Note that P-VC here denotes the Poiseuille-VC configuration to compare with P-RB. All simulations have fixed $Gr=10^7$, $Pr=1$. This figure is also available as an interactive JFM notebook: https://www.cambridge.org/S0022112024005986/JFM-Notebooks/files/Figure-13.ipynb.

In figures 13(c,d), we present the same analysis but for natural and mixed vertical convection. The contrast to the RB cases is immediately apparent, with no sharp peaks at any wavenumber or wall-normal position. This confirms the earlier visual observation of figure 3, where the instantaneous snapshots showed no clear coherent length scale. In natural VC, shown in figure 13(c), the two one-dimensional spectra do not overlap exactly as in the RB case due to the buoyancy-driven mean flow in the vertical. The $k_x$ spectrum shows a small peak at $kH\approx 3$ that is absent from the $k_z$ spectrum and becomes more prominent towards the channel centre. Comparing figures 13(c) and 13(d), we see that this peak is suppressed by the addition of external shear. Furthermore, all of the broad mid-range peaks in the spectra exhibit a decrease in amplitude and a shift to higher wavenumbers.

The heat flux co-spectra $\varPhi _{v\theta }$, shown in figure 14, exhibit similar features to the power spectra. In figures 14(a,b), a sharp low-wavenumber peak is once again visible for RB and mixed RB systems, and for the heat flux, this peak increases towards the channel centre. A broad high-wavenumber peak is also observed in the spectrum close to the wall. This peak flattens and shifts to lower wavenumbers (i.e. larger scales) as distance from the wall increases, which can be interpreted as the emergence and coalescence of small-scale plume structures. A similar behaviour can be found in the VC configurations of figures 14(c,d), with the broad peak in both the $k_x$ and $k_z$ spectra shifting to lower wavenumbers for larger $y$. However, in the $k_x$ spectrum for natural VC (where $x$ is perpendicular to gravity), this peak also increases in amplitude towards the channel centre, highlighting that the majority of the wall-normal heat transport at the channel centre occurs due to structures of size $\lambda \approx H$. Comparing this to the mixed VC case in figure 14(d), the mid-scale peak appears significantly suppressed when the external shear is imposed, suggesting that the coalescence of plumes in the bulk is disrupted by the shear. This is consistent with the interpretation of Domaradzki & Metcalfe (Reference Domaradzki and Metcalfe1988) and Scagliarini et al. (Reference Scagliarini, Gylfason and Toschi2014), who argued that the drop in heat flux in mixed RB can be related phenomenologically to disruption of the organisation of small-scale convective plumes.

Figure 14. One-dimensional co-spectra $\varPhi _{v\theta }(k,y)$ of the wall-normal heat flux for the extended domain simulations described in table 2. See the caption of figure 13 for more details on the meanings of line styles and colours. The spectra are normalised by $\sqrt {Ra\,Pr}$ such that integration over $k$ recovers the dimensionless advective heat flux $Nu-1 = \langle v\theta \rangle /(\kappa \varDelta /H)$. As in figure 13, the configurations presented are (a) RB, (b) mixed RB, (c) VC, (d) mixed VC. This figure is also available as an interactive JFM notebook: https://www.cambridge.org/S0022112024005986/JFM-Notebooks/files/Figure-14.ipynb.

7. Heat flux statistics in the boundary layer

Whereas the previous section focused on the advective heat flux in the bulk and its modification due to shear, we now turn to the conductive heat flux that dominates the transport in the boundary layer. Specifically, we investigate the statistical distribution of the local conductive heat flux at the boundary plates. As defined earlier, in (3.1), we consider the local dimensionless heat flux $q_\theta$ whose time and plane average is equivalent to the Nusselt number. The pointwise data of $q_\theta (x,z,t)$, also shown in figure 2, are collected over time and space to construct histograms that are normalised to produce probability density functions (p.d.f.s) in figure 15. The data shown here are only for $Pr=1$ and $Gr=10^6, 10^7$ since these cases cover the widest range of Richardson numbers, highlighting the transition from natural convection to forced convection.

Figure 15. Probability density functions (p.d.f.s) of the dimensionless local heat flux $q_\theta$ as defined in (3.1) for (a,c) $Gr=10^6$ and (b,d) $Gr=10^7$. Line colours highlight variation in the Reynolds number, as shown by the colour bar, with the blue dashed line plotting the data for natural convection with $Re=0$. Plots (a,b) are presented with linear $y$ axes, and (c,d) present the same data on logarithmic $y$ axes to highlight the exponential tails. This figure is also available as an interactive JFM notebook: https://www.cambridge.org/S0022112024005986/JFM-Notebooks/files/Figure-15.ipynb.

In figures 15(a,b), where the data are presented on linear axes, we find the same effect of increasing shear at both Grashof numbers. As $Re$ increases from zero (shown by progressively darker lines, with $Re=0$ marked in blue for comparison), the peak of the p.d.f. increases in amplitude. In statistical terms, this means that the most common values of local heat flux become more common as shear is increased, up to $Re\approx 10^{3.5}$. Due to the skewed nature of the distributions, these increasingly common values of heat flux are below the mean, and although the peak shifts slightly to the right as $Re$ increases, these parameters are associated with the drop in heat flux observed in figure 6. This corresponds to the visualisation of figure 2, where the streaks of low heat flux span a larger proportion of the boundary at $Ri=1$ (figures 2e,i) than at high $Ri$ (figures 2b,c). Once $Re$ is sufficiently high, which in these cases is for $Re> 10^{3.5}$, the whole distribution shifts more significantly to higher $q_\theta$ as the transport becomes dominated by the pressure-driven shear flow.

At moderate $Re$, another key modification is to the tails of the distribution. These are visualised most clearly in figures 15(c,d), where the heat flux distributions are plotted on a logarithmic axis. In this representation, it is clear that the right tails of large heat flux decay exponentially for $Gr=10^7$, and close to exponentially for $Gr=10^6$. At both Grashof numbers, the tails are reduced as $Re$ increases relative to the natural convection case, showing that the probability of extreme local heat flux values is reduced due to the introduction of a mean horizontal flow. Again, this is consistent with what was observed in the visualisation of figure 2, where dark patches associated with large local heat flux became less prevalent as $Re$ increases. In natural VC, Pallares et al. (Reference Pallares, Vernet, Ferre and Grau2010) used conditional sampling to show that these patches are most often associated with instantaneous flow reversals at the walls, which are likely a result of impacting plumes originating from the opposing wall.