1. Introduction
The interplay between buoyancy and shear in mixed thermal convection can be studied by either adding Couette-type forcing to the Rayleigh–Bénard (RB) system (Ahlers, Grossmann & Lohse Reference Ahlers, Grossmann and Lohse2009; Lohse & Xia Reference Lohse and Xia2010; Chilla & Schumacher Reference Chilla and Schumacher2012; Xia Reference Xia2013; Shishkina Reference Shishkina2021; Ahlers et al. Reference Ahlers2022; Lohse & Shishkina Reference Lohse and Shishkina2023) to obtain the Couette–RB (CRB) system (Deardorff Reference Deardorff1965; Ingersoll Reference Ingersoll1966; Hathaway & Somerville Reference Hathaway and Somerville1986; Domaradzki & Metcalfe Reference Domaradzki and Metcalfe1988; Solomon & Gollub Reference Solomon and Gollub1990; Shevkar et al. Reference Shevkar, Gunasegarane, Mohanan and Puthenveetil2019; Blass et al. Reference Blass, Zhu, Verzicco, Lohse and Stevens2020, Reference Blass, Tabak, Verzicco, Stevens and Lohse2021), or by applying a Poiseuille-type forcing to obtain the Poiseuille–RB (PRB) system (Scagliarini, Gylfason & Toschi Reference Scagliarini, Gylfason and Toschi2014; Zonta & Soldati Reference Zonta and Soldati2014; Scagliarini et al. Reference Scagliarini, Einarsson, Gylfason and Toschi2015; Pirozzoli et al. Reference Pirozzoli, Bernardini, Verzicco and Orlandi2017). A schematic of the two systems is shown in figure 1.
Figure 1. Schematic of the (a) CRB and (b) PRB systems.
The CRB and PRB systems are described by the incompressible Navier–Stokes equations, the continuity equation and the temperature transport equation, within the Boussinesq approximation. In Cartesian coordinates, they read
$$\begin{gather} {\partial_t u_i} + u_j {\partial_j u_i} =- \rho^{-1} {\partial_i p} + \nu {\partial^2_j u_i} + \beta g \delta_{i3}\theta + \varPi \delta_{i1}, \quad {\partial_i u_i} = 0, \end{gather}$$
$$\begin{gather} {\partial_t \theta} + u_j {\partial_j \theta} = \kappa {\partial^2_j \theta}, \end{gather}$$
where 
$\boldsymbol {u} \equiv (u_x,u_y,u_z)$ is the velocity, 
$p$ the pressure, 
$\theta$ the reduced temperature, 
$\rho$ the density of the fluid, 
$g$ the acceleration due to gravity antiparallel to the 
$z$ direction, 
$\beta$ the isobaric thermal expansion coefficient, 
$\nu$ the kinematic viscosity, 
$\kappa$ the thermal diffusivity and 
$H$ is the distance between the horizontal walls. At the top wall (
$z=H$), the reduced temperature is set to 
$\theta = -\varDelta /2$ while at the bottom wall (
$z=0$), the reduced temperature is set to 
$\theta = \varDelta /2$. For the CRB system 
$\varPi = 0$, the bottom wall is at rest and a velocity of 
$2U_w$ is imposed on the top wall. For the PRB system, no-slip conditions are enforced at the walls and a volume forcing 
$\varPi$ is applied in the streamwise direction such that it induces a bulk velocity of 
$U_b$ averaged over the domain volume and time (for the details of the implementation of the shear forcing in the numerical simulations, we refer the reader to § 3). The streamwise direction is oriented along 
$x$ and the spanwise direction along 
$y$. The aspect ratios of the system are defined by 
$\varGamma _x = L_x/H$ and 
$\varGamma _y = L_y/H$, with 
$L_x$, 
$L_y$ being the dimensions of the system in the 
$x$ and 
$y$ directions, respectively.
The control parameters for the systems are the Rayleigh number, the Prandtl number and the Reynolds number associated with the shear forcing
\begin{equation} Ra \equiv \frac{\beta g H^3 \varDelta}{\nu \kappa}, \quad Pr \equiv \frac{\nu}{\kappa}, \quad Re_S \equiv \frac{U_S H}{\nu} . \end{equation}
The characteristic velocity scale associated with the shear forcing 
$U_S$ is given by 
$U_S \equiv U_w$ for the CRB system and 
$U_S \equiv U_b$ for the PRB system. Although 
$U_b$ is formally a response parameter, in our numerical simulations the volume forcing term 
$\varPi$ is computed at each time step to ensure a constant mass flow rate (Quadrio, Frohnapfel & Hasegawa Reference Quadrio, Frohnapfel and Hasegawa2016) dictated by 
$U_b$, making it a control parameter in our case. The shear forcing for the CRB system is given by the wall Reynolds number 
$Re_w \equiv {U_w H}/{\nu }$ whereas the shear forcing for the PRB system is given by the bulk Reynolds number 
$Re_b \equiv {U_b H}/{\nu }$. Henceforth, we use 
$Re_S$ to indicate shear forcing in equations that are applicable to both CRB and PRB systems.
Similarly, one can also define a Reynolds number associated with the thermal forcing in these systems. From the input parameters, we can construct a Reynolds number 
${Re_F = U_F H /\nu \equiv \sqrt {Ra/Pr}}$, using the free fall velocity scale 
$U_F = \sqrt {g \beta \Delta H}$. However, the free fall scale is often not a reliable estimate of the flow velocities that develop in natural convection flows. A more appropriate approach is to define the Reynolds number associated with the large-scale convection (LSC) roll given by 
$Re_L \equiv U_L H / \nu$, with 
$U_L$ indicating the mean velocity of the ‘wind of turbulence’ generated by the LSC roll. The parameter 
$Re_L$ is, however, a response parameter whose variation with 
$Ra$, 
$Pr$ and 
$Re_S$ is not known a priori.
In the limiting case of zero shear forcing, we can distinguish the Reynolds number associated with the LSC roll in pure RB flow as 
$Re_R(Ra,Pr) \equiv Re_L(Ra,Pr,Re_S=0)$. The dynamics of the sheared RB systems are governed by a combined effect of both shear and thermal forcing. Therefore, we can also introduce the Reynolds number 
$Re_T$, which is constructed using the total velocity 
$U_T$ comprising a vector sum of 
$U_S$ and 
$U_L$. Naturally, the time-averaged wall shear stress 
$\tau _T$ generated by the total velocity 
$U_T$ is also determined by the combined effect of 
$\tau _L$, which is the time-averaged shear stress locally generated on the walls by the LSC roll, and 
$\tau _S$, which is the mean streamwise shear stress generated on the walls due to the applied shear forcing. These shear stresses can be expressed in dimensionless form using the friction coefficients associated with the total shear, the LSC roll and the streamwise shear respectively, as
\begin{equation} C_T \equiv \frac{2 \tau_T}{\rho U_T^2}, \quad C_L \equiv \frac{2 \tau_L}{\rho U_L^2}, \quad C_S \equiv \frac{2 \tau_S}{\rho U_S^2}. \end{equation}
Once again, in the limiting case of zero shear forcing, we can distinguish the friction coefficient associated with the LSC roll in pure RB flow as 
$C_R(Ra,Pr) \equiv C_L(Ra,Pr,Re_S=0)$. The non-dimensional heat flux from the hot bottom wall to the cold top wall is given by the Nusselt number
\begin{equation} Nu \equiv \frac{\left\langle u_z \theta - \kappa \partial_z \theta \right\rangle_{A,t}} {\kappa \Delta H^{-1}}, \end{equation}
with 
$\langle \rangle _{A,t}$ indicating the averaging in time and over any horizontal plane spanned by 
$x$ and 
$y$. A summary of all response and control parameters discussed above is given in table 1 for reference. In this study, we are primarily interested in understanding the dependence of the response parameters on the control parameters, and the physics underlying the connections between the response parameters.
Table 1. A summary of all control and response parameters discussed in § 1.
Using the exact relations for the global kinetic and thermal dissipation rates, Grossmann & Lohse (Reference Grossmann and Lohse2000, Reference Grossmann and Lohse2001, Reference Grossmann and Lohse2002, Reference Grossmann and Lohse2004) offered a unifying theory for RB convection (hereafter referred to as the GL-theory). For cylinders of unit aspect ratio, Stevens et al. (Reference Stevens, van der Poel, Grossmann and Lohse2013) have demonstrated that fitting the GL-theory at four data points from Funfschilling et al. (Reference Funfschilling, Brown, Nikolaenko and Ahlers2005) at 
$Ra = 2.96 \times 10^7$ and 
$Ra = 1.92 \times 10^{10}$ with 
$Pr = 4.38$, from Xia, Lam & Zhou (Reference Xia, Lam and Zhou2002) at 
$Ra = 2.24 \times 10^8$ with 
$Pr = 554$ and from Kerr & Herring (Reference Kerr and Herring2000) at 
$Ra = 10^7$ with 
$Pr = 0.07$ using four free parameters can predict 
$Nu$ within 
$4\,\%$ of experimental and numerical results in most of the parameter space given by 
$10^{4} \leq Ra \leq 10^{14}$ and 
$10^{-3} \leq Pr \leq 10^{2}$, with only two small ranges that exhibit a greater than 
$10\,\%$ disagreement. In the same paper Stevens et al. (Reference Stevens, van der Poel, Grossmann and Lohse2013), a fit of similar quality is achieved also for an aspect ratio 
$\varGamma = 1/2$, with slightly different prefactors. This work has further been extended by Ahlers et al. (Reference Ahlers2022) to include the effects of aspect ratios between 
$1/32$ and 
$32$. The GL-theory has also been extended to the ultimate regime (Grossmann & Lohse Reference Grossmann and Lohse2011), where the heat transport is considerably enhanced, as the laminar-type boundary layers become turbulent due to a non-normal–nonlinear instability; for a detailed discussion see Roche (Reference Roche2020) and Lohse & Shishkina (Reference Lohse and Shishkina2023).
Presently, the GL-theory has been applied to RB convection without imposed shear. The objective of this work is to extend the theoretical approach to sheared thermal convection. Scagliarini et al. (Reference Scagliarini, Gylfason and Toschi2014, Reference Scagliarini, Einarsson, Gylfason and Toschi2015); 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) have made progress in understanding the variation of heat transfer in RB convection with imposed shear. Scagliarini et al. (Reference Scagliarini, Gylfason and Toschi2014, Reference Scagliarini, Einarsson, Gylfason and Toschi2015) proposed a model based on the concept of eddy viscosity and eddy diffusivity to explain the counter-intuitive initial decrease and subsequent increase in 
$Nu$ with increasing 
$Re_b$ for the PRB system. Blass et al. (Reference Blass, Zhu, Verzicco, Lohse and Stevens2020) observed a similar effect in the CRB system and attributed it to the initial destruction of the large-scale flow organisation and the subsequent formation of large meandering flow structures (Hutchins & Marusic Reference Hutchins and Marusic2007). They divided the flow into a buoyancy-dominated, a transitional and a shear-dominated regime, based on the Monin–Obukhov length scale. Blass et al. (Reference Blass, Tabak, Verzicco, Stevens and Lohse2021) further investigated the effect of 
$Pr$ on the variation of 
$Nu$ with 
$Re_w$ and concluded that the non-monotonic behaviour of 
$Nu(Re_w)$ is a consequence of flow layering, plume sweeping and bulk heat entrapment. Building on these findings, in this paper we propose a more general formulation applicable to sheared thermal convection, in the spirit of the GL-theory.
The objective of this paper is to extend the GL-theory to CRB and PRB systems by providing scaling relations for 
$Nu(Ra,Pr,Re_S)$ and 
$C_S(Ra,Pr,Re_S)$. Furthermore, we will show similarities in the dependence of 
$C_S$ on 
$Re_S$ in Couette and Poiseuille flows and 
$C_L$ on 
$Re_L$ in RB convection, suggesting some sort of universality in the behaviour of the friction coefficient in shear-driven flows and thermal convection. In § 2, we build the theoretical framework, using rigorous relations for globally averaged kinetic and thermal dissipation rates. In § 3, we validate the theoretical scaling relations against direct numerical simulations (DNS). Finally, the conclusions are presented in § 4.
2. Extending the Grossmann–Lohse theory to CRB and PRB
2.1. Kinetic and thermal dissipation rates
To extend the GL-theory to sheared thermal convection, we formulate exact global relations for the kinetic 
$(\epsilon _u)$ and thermal 
$(\epsilon _\theta )$ dissipation rates in the CRB and PRB systems. These arise from the time- and volume-averaged equations for the kinetic energy and temperature variance, respectively. The relation for the mean thermal dissipation rate is the same as in the classical RB convection
\begin{equation} \epsilon_\theta = \left\langle \kappa {(\partial_j \theta)}^2 \right\rangle_{V,t} = \kappa\frac{\varDelta^2}{H^2}Nu, \end{equation}see e.g. Shraiman & Siggia (Reference Shraiman and Siggia1990) and Siggia (Reference Siggia1994). The relation for the mean kinetic dissipation rate reads
\begin{equation} \epsilon_u = \left\langle \nu {(\partial_j u_i)}^2 \right\rangle_{V,t} = \frac{\nu^3}{H^4}\left( \underbrace{(Nu-1) Ra Pr^{-2}}_{\text{Buoyancy term}} + \underbrace{C_S Re_S^3}_{\text{Shear term}}\right), \end{equation}
with 
$\langle \dots \rangle _{V,t}$ indicating the average over time and volume. Note that the expression for thermal dissipation (2.1) is the same as that in classical RB convection but the expression for kinetic dissipation (2.2) includes contributions from both buoyancy and shear forcing.
2.2. Kinetic energy, large-scale convection rolls and boundary layer thickness
One of the central ideas of the GL-theory is the presence of persistent LSC rolls that churn the bulk of the system and generate boundary layers at the walls. As a result, the mean kinetic energy of the RB system is expected to scale as 
${\sim }U_R^2$, where 
$U_R$ is the velocity scale of the LSC roll (in the absence of any shear forcing). In pure RB flow, this mean kinetic energy is solely generated by the buoyancy forcing. However, in the case of sheared RB flow, where the LSC roll has a velocity scale 
$U_L$, the mean kinetic energy consists of contributions from both the LSC roll and the imposed shear flow. Therefore, we add the kinetic energy of the mean flow 
$U_S^2$ and the associated turbulent kinetic energy (TKE), which scales as the square of the friction velocity 
$u_\tau =\sqrt {\tau _S/\rho }$, to write
\begin{equation} U_T^2 \approx U_L^2 + U_S^2 + 2\gamma u_\tau^2, \end{equation}
with 
$\gamma$ being a prefactor. This can also be written in terms of the corresponding Reynolds numbers as
\begin{equation} Re_T^2 \approx Re_L^2 + Re_S^2 + \gamma C_S Re_S^2. \end{equation}
For a laminar Prandtl–Blasius-type (Prandtl Reference Prandtl1904; Blasius Reference Blasius1908) boundary layer, there is no TKE, so the contribution 
$C_S Re_S^2$ vanishes and we can approximate (2.4) as
\begin{equation} Re_T^2 \approx Re_L^2 + Re_S^2. \end{equation}
An interpretation of the above equation (2.5) is that the velocity associated with the LSC roll is preferentially oriented along the direction orthogonal to the shear forcing, consistent with the flow structures observed by Pirozzoli et al. (Reference Pirozzoli, Bernardini, Verzicco and Orlandi2017); Blass et al. (Reference Blass, Zhu, Verzicco, Lohse and Stevens2020). Note that the validity of (2.5) is limited to sufficiently low values of 
$Re_S$ wherein the contribution to spanwise shear stresses from shear forcing is negligible in comparison with the contribution from the LSC rolls. At high shear forcing, the spanwise shear stresses generated by velocity fluctuations arising purely from shear forcing may no longer be negligible, in which case (2.5) no longer holds.
In the buoyancy-dominated regime, relation (2.5) can be better understood by considering many LSC rolls each orientated at an angle 
$\alpha$ with the streamwise direction and studying the probability distribution of 
$\alpha$ given by 
$\phi (\alpha )$. Here, we mean LSC rolls in a broad sense, without addressing the exact details of the flow organisation at this stage. Empirical observations regarding flow organisation are reported in section § 3.3. Since the total velocity 
$U_T$ arises from a vector addition of the shear velocity and the LSC velocities, we can express it in terms of 
$\phi (\alpha )$ as
\begin{equation} U_T^2 = \int_{0}^{2{\rm \pi}} { \phi(\alpha) \left(U_S^2 + U_L^2 + 2 U_S U_L \cos(\alpha) \right) {\rm d}\alpha}. \end{equation}
Due to the symmetry of the system and the periodic boundary conditions in the horizontal directions, it is reasonable to consider that there are an equal number of clockwise and counter-clockwise LSC rolls within the sheared RB system. When averaged over the entire volume of the system, we postulate that the probability distribution 
$\phi (\alpha )$ is symmetric about the spanwise direction, i.e. about 
$\alpha = \pm {\rm \pi}/2$. Applying this symmetry condition to (2.6) gives us the relation (2.5). The velocity scale 
$U_T$ associated with the kinetic energy of the system can thus be considered as a vector sum of perpendicular velocity contributions from the shear forcing in the streamwise direction and the LSC roll in the spanwise direction. Following this approach, we can also decompose the total shear stress 
$\tau _T$ generated by 
$U_T$ into the streamwise component 
$\tau _s = \tau _T U_S/U_T$, generated by 
$U_S$, and the spanwise component 
$\tau _L = \tau _T U_L/U_T$, generated by 
$U_L$. This assumption gives us three equivalent definitions of the kinetic boundary layer thickness 
$\lambda _u$, namely
\begin{equation} \lambda_u \equiv \frac{2H}{C_L Re_L} = \frac{2H}{C_S Re_S} = \frac{2H}{C_T Re_T}. \end{equation}
Here, we used the slope criterion from Shishkina et al. (Reference Shishkina, Stevens, Grossmann and Lohse2010) for the definition of the kinetic boundary layer thickness. Similarly, we define the thermal boundary layer thickness 
$\lambda _\theta$ also with the slope criterion as
\begin{equation} \lambda_\theta \equiv \frac{H}{2Nu}. \end{equation}2.3. Estimating bulk and boundary layer contributions
Another key idea of the GL-theory is the splitting of 
$\epsilon _u$ and 
$\epsilon _\theta$ into their bulk and boundary layer contributions as follows:
\begin{equation} \epsilon_u = \epsilon_{u,BL} + \epsilon_{u,bulk}, \quad \epsilon_\theta = \epsilon_{\theta,BL} + \epsilon_{\theta,bulk}, \end{equation}
with 
$\epsilon _{u,BL}$ being the boundary layer contribution to the kinetic dissipation, 
$\epsilon _{u,bulk}$ being the bulk contribution to the kinetic dissipation, 
$\epsilon _{\theta,BL}$ being the boundary layer contribution to the thermal dissipation and 
$\epsilon _{\theta,bulk}$ being the bulk contribution to the thermal dissipation. Therefore, we focus on estimating 
$\epsilon _{u,BL}$, 
$\epsilon _{u,bulk}$, 
$\epsilon _{\theta,BL}$ and 
$\epsilon _{\theta,bulk}$ for sheared RB systems.
First, we estimate 
$\epsilon _{u,BL}$ using (2.7) as
\begin{equation} \epsilon_{u,BL} \sim \nu \frac{U_T^2}{\lambda_u^2} \frac{\lambda_u}{H} \sim \frac{\nu^3}{H^4} C_T Re_T^3 \equiv \frac{\nu^3}{H^4} \left(\underbrace{C_L Re_L^3}_{\text{LSC term}} + \underbrace{C_S Re_S^3}_{\text{Shear term}}\right), \end{equation}
which is, in turn, a sum of contributions from the LSC rolls and the applied shear forcing. Note that the shear contribution here exactly matches that in the global relation (2.2). Next, we estimate 
$\epsilon _{u,bulk}$. Here, it is important to note that the bulk dissipation rate is dominated by the contribution from LSC rolls, while the contribution from applied shear forcing is much smaller. Therefore, we only focus on estimating the contribution to 
$\epsilon _{u,bulk}$ from the LSC rolls by assume that the velocity scale 
$U_L$ associated with the LSC rolls is responsible for stirring the bulk with a kinetic energy that scales with 
$U_L^2$. However, LSC rolls are swept at the boundary layer height due to the applied shear forcing and the time scale of the stirring process is governed not by the velocity 
$U_L$ associated with the LSC rolls but by 
$U_T$ which is the total velocity scale. As we shall find later in § 3.2, this is a key assumption that explains the trend of 
$Nu$ with increasing 
$Re_S$ in the buoyancy-dominated regime. Following these assumptions, we write
\begin{equation} \epsilon_{u,bulk} \sim U_L^2 \frac{U_T}{H} \equiv Re_L^2 Re_T. \end{equation}
At high-shear forcing, the contribution of the boundary layer to the dissipation rate arising from shear forcing dominates the bulk contribution from LSC rolls. In the limiting case of passive transport in Couette/Poiseuille flow, one can rigorously derive that the total kinetic dissipation rate 
$\epsilon _u = (\nu ^3/H^4) C_S Re_S^3$. In this sense, the kinetic dissipation rate of sheared RB at high shear forcing will always be dominated by the boundary layer contribution and estimating 
$\epsilon _{u,bulk}$ in the shear-dominated regimes is redundant.
The analogous estimate for 
$\epsilon _{\theta,BL}$
\begin{equation} \epsilon_{\theta,BL} \sim \kappa \frac{\varDelta^2}{\lambda_\theta^2} \frac{\lambda_\theta}{H} \sim \kappa\frac{\varDelta^2}{H^2} Nu, \end{equation}is identical to the exact relation (2.1), on the one hand showing consistency of the approach, but on the other hand not giving new information. Therefore, following the GL-theory, we match the magnitude of the advective and diffusive terms of (1.2) at the thermal boundary layer height to obtain
\begin{equation} u_y \partial_y \sim \kappa \partial_{zz}. \end{equation}
As in the GL-theory, for regimes where the thermal boundary layer is thicker than the kinetic boundary layer (
$\lambda _\theta > \lambda _u$, associated with low 
$Pr$), we estimate 
$u_y \sim U_L$, 
$\partial _y \sim 1/H$, and 
$\partial _{z z} \sim \lambda _\theta ^{-2}$. Using these estimates in (2.13) with 
$\lambda _u$ from (2.7) and 
$\lambda _\theta$ from (2.8), we obtain
\begin{equation} Nu \sim Pr^{1/2} Re_L^{1/2} \equiv \frac{Pr^{1/2}}{C_L Re_L^{1/2}} C_T Re_T. \end{equation}
For high 
$Pr$ regimes, where 
$\lambda _\theta < \lambda _u$, we estimate 
$u_y \sim U_L \lambda _\theta / \lambda _u$ due to the fact that the relevant velocity scale at the thermal boundary height is smaller than the velocity 
$U_L$ by a factor 
$\lambda _\theta / \lambda _u$, exactly as in the GL-theory. Once again using the estimates 
$\partial _y \sim 1/H$, and 
$\partial _{z z} \sim \lambda _\theta ^{-2}$ in (2.13) with 
$\lambda _u$ from (2.7) and 
$\lambda _\theta$ from (2.8), we get
\begin{equation} Nu \sim Pr^{1/3} Re_L^{1/3} \left(C_T Re_T\right)^{1/3} \equiv \frac{Pr^{1/3}}{C_L^{2/3} Re_L^{1/3}} C_T Re_T. \end{equation}
Also the bulk contribution to the thermal dissipation is estimated in an identical manner to the corresponding equation in the GL-theory, only with minor changes to reflect the new dependence on 
$C_L$, 
$C_T$ and 
$Re_T$. For 
$\lambda _\theta > \lambda _u$
\begin{equation} \epsilon_{\theta,bulk} \sim \varDelta^2\frac{U_L}{H} \sim \kappa\frac{\varDelta^2}{H^2} Pr Re_L \equiv \kappa\frac{\varDelta^2}{H^2} \frac{Pr}{C_L} \left(C_T Re_T\right), \end{equation}
and for 
$\lambda _\theta < \lambda _u$, where the relevant velocity scale is 
$U_L \lambda _\theta /\lambda _u$, we get
\begin{equation} \epsilon_{\theta,bulk} \sim \varDelta^2\frac{U_L}{H} \frac{\lambda_\theta}{\lambda_u} \sim \kappa\frac{\varDelta^2}{H^2} Pr C_L Re_L^2 Nu^{-1} \equiv \kappa\frac{\varDelta^2}{H^2} \frac{Pr}{C_L} \frac{\left(C_T Re_T\right)^2}{Nu}. \end{equation}
We expect (2.14)–(2.17) to be valid in both buoyancy-dominated and shear-dominated regimes. This is accommodated by the fact that the dependencies 
$C_T(Re_T)$ and 
$C_L(Re_L)$ behave differently in the buoyancy-dominated and shear-dominated regimes. Compatibility of these relations in the limiting regimes is explained in the subsequent § 2.4.
2.4. Limiting regimes
It is important to note that there exist no pure power laws for 
$Nu$, 
$Re_L$ and 
$C_S$ as functions of 
$Ra$, 
$Pr$ and 
$Re_S$. Nonetheless, it is useful to study the pure scaling power laws that arise from limiting regimes, in the interest of understanding the physics of the system. Based on the dominance of boundary layer and bulk contributions to the kinetic and thermal dissipation rates, the GL-theory provides four regimes 
$I$, 
$II$, 
$III$ and 
$IV$ for the pure RB system. Furthermore, each regime can be divided into two subregimes based on whether the thermal boundary layer is nested into the kinetic boundary layer or vice versa. We find that this classification of regimes is also applicable to buoyancy-dominated sheared RB. The phase space of 
$Ra$, 
$Pr$ and 
$Re_S$ is divided by four different transitions – (i) transition from boundary layer dominated regimes to bulk-dominated regimes, (ii) transition from buoyancy-dominated regime to shear-dominated regime, (iii) transition between regimes where the thermal boundary layer is nested inside the kinetic boundary layer or vice versa and (iv) transition from a laminar Prandtl–Blasius-type boundary layer to a turbulent Prandtl–von Kármán-type boundary layer.
It should also be noted that some of these limiting regimes may not exist for sheared RB in the shear-dominated state. For example, at high shear forcing, the shear term in 
$\epsilon _{u,BL}$ always dominates 
$\epsilon _u$ because shear forcing primarily increases the boundary layer contribution of the kinetic dissipation 
$(\nu ^3/H^4)C_S Re_S^3$. Since the friction coefficient 
$C_S$ associated with the applied shear can become independent of 
$Re_S$ only asymptotically at infinite 
$Re_S$, shear-dominated systems with 
$\epsilon _u \sim \epsilon _{u,bulk}$ cannot exist.
We now proceed to first analyse the shear-dominated regime with 
$Re_S/Re_L \gg 1$. This can be realised in two ways – namely either 
$Re_L$ is small or 
$Re_L$ is not necessarily small but 
$Re_S \gg Re_L$. When 
$Re_L$ is small, we assume the existence of a laminar Prandtl–Blasius-type boundary layer with
\begin{equation} C_L \sim Re_L^{-1/2}, \end{equation}which, along with (2.7) and (2.18), gives
\begin{equation} C_S \sim Re_L^{1/2} Re_S^{-1}. \end{equation}
With the assumption (2.18), we see that (2.14) associated with small 
$Pr$ becomes
\begin{equation} Nu \sim Pr^{1/2} C_T Re_T, \end{equation}
and (2.15) associated with large 
$Pr$ becomes
\begin{equation} Nu \sim Pr^{1/3} C_T Re_T. \end{equation}
For the limiting case of 
$Re_L=0$ one can see that relations (2.20) and (2.21) recover the scaling laws for passive transport in Couette (Yerragolam et al. Reference Yerragolam, Stevens, Verzicco, Lohse and Shishkina2022a) or Poiseuille flow (Kays & Crawford Reference Kays and Crawford1993) where the relation between 
$C_T$ and 
$Re_T$ depends on whether the kinetic boundary layer is laminar or turbulent. In the presence of a laminar boundary layer, the trivial scaling, 
$C_T \sim Re_T^{-1}$, applies. When the boundary layer turns turbulent with increased shear forcing, the relation between 
$C_T$ and 
$Re_T$ is given by Prandtl (Reference Prandtl1932) friction law obtained from the log-law mean velocity profile which states
\begin{equation} \sqrt{\frac{2}{C_T}} = \frac{1}{k} \ln{ \left( Re_T \sqrt{\frac{C_T}{8}} \right) } + B, \end{equation}
with 
$k \approx 0.41$ (Pirozzoli, Bernardini & Orlandi Reference Pirozzoli, Bernardini and Orlandi2014) being the von Kármán (Reference von Kármán1934) constant and 
$B\approx 5$ (Pirozzoli et al. Reference Pirozzoli, Bernardini and Orlandi2014) indicating the log-law intercept.
When 
$Re_L$ is not necessarily small and 
$Re_S \gg Re_L$, we consider that the passive transport relations of (2.20) or (2.21) remain relevant, and that the dependence of 
$C_T(Re_T)$ is unchanged to that described above. For this case of larger 
$Re_L$, we assume that the thermal dissipation is dominated by contributions from the bulk, so that the global dissipation relation (2.1) can be estimated by (2.16) or (2.17). When we compare the passive transport relations with the dissipation estimates, we find that 
$C_L$ must become independent of 
$Re_L$, with 
$C_L\sim Pr^{1/2}$ for low 
$Pr$ and 
$C_L\sim Pr^{1/3}$ for high 
$Pr$. Furthermore, the behaviour of 
$Re_L$ in these cases can be revealed by combining the passive transport relations (2.20) or (2.21) with the boundary-layer estimate for the kinetic dissipation rate (2.10). For low 
$Pr$ this produces 
$Re_L\sim Ra^{1/2} Pr^{-3/4}$, and for high 
$Pr$ we get 
$Re_L\sim Ra^{1/2} Pr^{-5/6}$, which exactly match the Reynolds number scaling relations found in the boundary layer dominated regime 
$I$ of the GL-theory for classical RB convection.
Now, we focus on the buoyancy-dominated regimes where 
$Re_S/Re_L \lesssim 1$ and 
$C_S Re_S^3 \ll (Nu-1) Ra Pr^{-2}$ such that the shear contribution to the kinetic dissipation can be neglected in comparison with the buoyancy contribution. With this restriction, we approximate (2.2) as
\begin{equation} \epsilon_u \approx \frac{\nu^3}{H^4} (Nu-1) Ra Pr^{-2}, \end{equation}and (2.10) as
\begin{equation} \epsilon_{u,BL} \approx \frac{\nu^3}{H^4} C_L Re_L^3. \end{equation}
Using these approximations, we can provide scaling relations between 
$Nu/Nu_R$ and 
$Re/Re_R$ for buoyancy-dominated sheared RB system with 
$Nu_R$ and 
$Re_R$ being the Nusselt number and Reynolds number associated with the LSC rolls for the pure RB system. Following the GL-theory, there are various regimes that can be relevant depending on whether the dissipation rates are dominated by boundary layer or bulk contributions, and whether the thermal BL is thicker than the kinetic BL. For each of these regimes, we combine (2.7) and (2.23) with the relevant estimates for the dominant dissipation rate contribution to give expressions for 
$Nu(Ra,Pr,C_L, Re_L, Re_T)$ and 
$Re_L(Ra,Pr,C_L,Re_L, Re_T)$. Since our approach is consistent with the GL-theory, the 
$Ra$ and 
$Pr$ dependence simply recovers the scaling relations 
$Nu_R(Ra,Pr)$ and 
$Re_R(Ra,Pr)$ found for the various regimes of pure RB, providing us with scaling relations for 
$Nu/Nu_R$ and 
$Re_L/Re_R$ that only depend on 
$C_L$, 
$Re_L$ and 
$Re_T$. In table 2, we outline the relevant estimates for the dissipation rates and present the resulting scaling relations. Is it important to note that these scaling relations are only applicable to buoyancy-dominated sheared RB flows with scaling-wise laminar Prandtl–Blasius-type BLs.
Table 2. Scaling relations for the Nusselt number 
$Nu$ and LSC Reynolds number 
$Re_L$ in the buoyancy-dominated regime of sheared RB convection with a laminar Prandtl–Blasius-type boundary layer. The first column indicates the GL regime, the second column indicates the bulk or boundary layer (BL) dominance of the kinetic dissipation rate with the applicable scaling estimate in the parenthesis and the third column indicates the bulk or BL dominance of the thermal dissipation rate with the applicable scaling estimate in the parenthesis. The fourth column indicates whether the kinetic BL is thicker than the thermal boundary layer or vice versa. The fourth and fifth columns indicate the scaling relations for 
$Nu/Nu_R$ and 
$Re_L/Re_R$ using the values of 
$Nu$ and 
$Re_R$ estimated for the pure RB system from GL-theory.
In the buoyancy-dominated classical GL regimes, we can consider the Prandtl–Blasius scaling (2.18) to hold for small shear forcing (i.e. 
$Re_S \lesssim Re_L$). With this assumption, 
$C_L^2 Re_L \approx 1$, so the values of 
$Nu$ and 
$Re_L$ remain unchanged for buoyancy-dominated regimes 
$I$ and 
$III$, whereas for buoyancy-dominated regimes 
$II$ and 
$IV$, the non-monotonic behaviour of 
$Nu$ with increasing 
$Re_S$ becomes apparent. Although 
$Nu$ seems to decrease with increasing 
$Re_S$ in the buoyancy-dominant 
$II$ and 
$IV$ regimes, it is important to note that this behaviour is subject to the condition that the BL is a laminar one of the Prandtl–Blasius type. If the BL becomes turbulent, the expected decrease in 
$Nu$ in the buoyancy-dominated regime might disappear. In this study, we will explore the 
$Nu$ response in the buoyancy-dominated 
$II_u$ regime where we can still observe the decrease in 
$Nu$ with increasing 
$Re_S$ within reasonable computational cost.
3. Results from the direct numerical simulations
3.1. Scheme and procedure
In this section, we will compare the scaling relations derived in the previous section against the results from our DNS. Equations (1.1a,b) and (1.2) are solved numerically using the in-house open-source code ‘AFiD’, which is based on a second-order finite-difference scheme (van der Poel et al. Reference van der Poel, Ostilla-Mónico, Donners and Verzicco2015). The code has been extensively validated (Verzicco & Orlandi Reference Verzicco and Orlandi1996; Verzicco & Camussi Reference Verzicco and Camussi1997; Stevens, Verzicco & Lohse Reference Stevens, Verzicco and Lohse2010; Stevens, Lohse & Verzicco Reference Stevens, Lohse and Verzicco2011; Kooij et al. Reference Kooij, Botchev, Frederix, Geurts, Horn, Lohse, van der Poel, Shishkina, Stevens and Verzicco2018). We impose periodic boundary conditions in the horizontal directions and no-slip boundary conditions at the top and bottom walls. For most simulations, we use domains of aspect ratios 
$\varGamma _x=8$ and 
$\varGamma _y=4$. We also performed CRB simulations with 
$\varGamma _x=48$ and 
$\varGamma _y=24$ for 
$Ra=10^7$, 
$Pr=1$ to study large-scale flow structures. Due to the need for high resolution at large 
$Ra$, the RB simulations for 
$Ra=10^{10}$, 
$Pr=1$ and 
$Ra=10^{11}$, 
$Pr=1$ were performed in domains of aspect ratios 
$\varGamma _x=\varGamma _y=4$, while the RB simulation for 
$Ra=10^{12}$, 
$Pr=1$ was performed in domain of aspect ratios 
$\varGamma _x=\varGamma _y=2$. For the CRB simulations, the wall velocities 
$(-U_w)$ and 
$U_w$ were imposed as Dirichlet boundary conditions on the bottom and top walls, respectively. This is done for numerical reasons (Bernardini et al. Reference Bernardini, Pirozzoli, Quadrio and Orlandi2013) and does not affect the analysis of the results. The equivalent velocity fields of the CRB system with the bottom wall at rest and the top wall at 
$2U_w$ can be obtained by a simple Galilean transformation, i.e. by adding 
$U_w$ to the numerically obtained flow field. For the PRB simulations, the volume forcing term 
$\varPi$ is computed at each time step to ensure a constant mass flow rate (Quadrio et al. Reference Quadrio, Frohnapfel and Hasegawa2016).
We use a uniform discretisation in the horizontal, periodic directions and a non-uniform grid in the wall-normal direction, in which we employ higher grid resolution in the BLs next to the walls. The thermal BL was ensured to be sufficiently resolved according to the resolution requirements put forward by Shishkina et al. (Reference Shishkina, Stevens, Grossmann and Lohse2010). The near-wall resolution is comparable to that of Lozano-Durán & Jiménez (Reference Lozano-Durán and Jiménez2014); Pirozzoli et al. (Reference Pirozzoli, Bernardini and Orlandi2014); Lee & Moser (Reference Lee and Moser2018) to ensure that the kinetic BL is sufficiently resolved. The simulations were run for a long enough physical time for the standard deviation of 
$Nu$ to converge within approximately 1 % of its mean value. In our previous work (Yerragolam et al. Reference Yerragolam, Verzicco, Lohse and Stevens2022b), we verified that the 
$Nu$ and 
$Re_{\tau }$ obtained from a domain with 
$\varGamma _x=8$ and 
$\varGamma _y=4$ shows a difference of less than 1 % from the 
$Nu$ obtained from a domain with 
$\varGamma _x=48$ and 
$\varGamma _y=24$. This observation is also supported by the fact that 
$Nu$ for the RB system converges at an approximate aspect ratio 
$\varGamma _x=\varGamma _y=4$ (Stevens et al. Reference Stevens, Blass, Zhu, Verzicco and Lohse2018).
Since many of the scaling relations in § 2.4 rely on the value of the wind Reynolds number 
$Re_R$, obtaining an estimate for 
$Re_R$ from the numerical simulations of pure RB is necessary. The value of 
$Re_R$ can be estimated in two possible ways. The first estimate can be obtained by using 
${Re_{max} = U_{max}H/\nu }$ where 
$U_{max}$ is the maximum value of the root mean squared (r.m.s.) horizontal velocity profile 
${u_R(z) \equiv \sqrt {\langle u_x^2 + u_y^2\rangle _{A,t}}}$ at a height 
$z = \delta$ as shown in figure 2(a). The second estimate can be obtained by using the global r.m.s. velocity 
$Re_{rms} = U_{rms}H/\nu$ with 
$U_{rms} \equiv \sqrt {\langle u_x^2 + u_y^2 + u_z^2\rangle _{V,t}}$. In figure 2(b) and 2(c), we can see that these estimates are almost identical, providing strong evidence that the LSC rolls driving the wind at the wall also provide the dominant contribution to the mean kinetic energy in RB convection. In all the results discussed henceforth, we adopt 
$Re_{rms}$ as an estimate for 
$Re_R$. Unlike in Couette or Poiseuille flow, the mean shear stress at the wall is zero in RB convection. However, we can use the mean gradient of the r.m.s. horizontal velocity 
$\langle \partial _z u_R(z)\rangle _{W,t}$ to calculate the friction coefficient 
$C_R$ associated with the large-scale circulation. Here, 
$\langle \cdots \rangle _{W,t}$ indicates time averaging over the surface of the walls. The variation of 
$C_R$ with 
$Re_R$ is discussed separately in § 3.5.
Figure 2. (a) The r.m.s. horizontal velocity profile 
$u_R(z)$ normalised with the free-fall velocity for RB system with 
$Ra=10^7$, 
$Pr=1.0$ is shown as a function of wall-normal height 
$z$, indicating the maximum value 
$U_{max}$ occurring at a wall-normal height 
$\delta$. (b) Value of 
$Re_{max}$ obtained from 
$U_{max}$ is plotted against the globally averaged r.m.s. velocity 
$Re_R$ of the RB system. The dashed grey line indicates 
$Re_{max} = Re_{rms}$, showing that these two estimates are virtually identical. (c) 
$Re_{max}$ compensated with 
$Re_{rms}$ plotted against 
$Re_{rms}$, once again emphasising that the two estimates are equivalent within the 5 decades of 
$Ra$ and nearly 3 decades of 
$Re_{rms}$ investigated here. The values of 
$Ra$ are shown on the top for reference. This figure is also available at https://www.cambridge.org/S0022112024008723/JFM-Notebooks/files/figure2.
3.2. Global response parameters
We now validate the scaling relations for 
$Nu$ and 
$C_S$ derived in § 2.4. Within the parameter range of 
$Ra$, 
$Pr$ and 
$Re$ simulated, we already observe multiple transitions. As shear forcing is increased, we undergo transition from the buoyancy-dominated regime to the shear-dominated regime. For low shear forcing, 
$\epsilon _u \sim \epsilon _{u,bulk}$ with 
$\lambda _\theta < \lambda _u$ whereas for high shear forcing, 
$\epsilon _u \sim \epsilon _{u,BL}$ with 
$\lambda _\theta /\lambda _u \sim Pr^{1/2}$ (Yerragolam et al. Reference Yerragolam, Stevens, Verzicco, Lohse and Shishkina2022a). Additionally, at low shear forcing, we observe a laminar Prandtl–Blasius-type BL which undergoes a transition into a turbulent one at high shear.
In the buoyancy-dominated regime, we assume the presence of Prandtl–Blasius-type kinetic boundary layer with the friction coefficient 
$C_L$ given by (2.18). For the parameter range of our simulations, the relevant convection regime is 
$II_u$, so combining the relevant relation from table 2 with (2.18), we arrive at
\begin{equation} \frac{Nu}{Nu_R} \sim \left(\frac{Re_T}{Re_L}\right)^{-1/5}. \end{equation}
In the buoyancy-dominated regimes, we also take (2.19) for the friction coefficient 
$C_S$ associated with the imposed shear. In order to further simplify these equations, we approximate 
$Re_L \approx Re_R$ in the buoyancy-dominated regime, which is justified by the weak variation of 
$Re_L$ with increasing 
$Re_S$ in (table 2). With this assumption, we can rewrite (2.19) as
\begin{equation} C_S \sim \sqrt{Re_R}/Re_S, \end{equation}and use (2.5) to rewrite (3.1) as
\begin{equation} Nu/Nu_R \sim \left( \sqrt{1+\left(Re_S^2/Re_R^2\right)} \right)^{-1/5}. \end{equation} These equations show good agreement with the numerical data plotted in figures 3(a)–3(d) for the buoyancy-dominated regime. Note that (3.1) and (3.3) do not explicitly state the dependence of 
$Nu_R(Ra, Pr)$ or 
$Re_R(Ra, Pr)$. There is no pure scaling exponent for 
$Nu_R(Ra)$ in pure RB for regime 
$II_u$. The values of 
$Nu_R$ and 
$Re_R$ used for figures 3(b) and 3(d) are obtained from numerical simulations of pure RB. The present extension to the GL-theory assumes that the values of 
$Nu_R$ and 
$Re_R$ are known a priori, and only attempts to provide scaling relations for the normalised quantities 
$Nu/Nu_R (Re_S/Re_R)$ and 
$C_S/C_R (Re_S/Re_R)$ in the buoyancy-dominated regime.
Figure 3. (a) Value of 
$Nu/Pr^{1/2}$ plotted against 
$Re_S$. The black solid line indicates (3.5), the dashed lines indicate (2.20) and the dotted lines indicate (3.3). (b) Value of 
$Nu/Nu_R$ plotted against 
$Re_S/Re_R$. Dotted lines indicate (3.3). (c) Value of 
$C_S$ plotted against 
$Re_S$. The black solid line indicates (3.4) and coloured dashed lines indicate (2.19). (d) Value of 
$C_S$ normalised with 
$C_R$ plotted against 
$Re_S/Re_R$. Coloured dashed lines indicate (3.4) and black dashed line indicates (2.19). The data for PRB are indicated with plus markers and the data for CRB are indicated with dot markers. This figure is also available at https://www.cambridge.org/S0022112024008723/JFM-Notebooks/files/figure3.
For the shear-dominated regime, we observe that the BL becomes turbulent. In this case, (2.19) is no longer valid. Instead, the relation between 
$C_T$ and 
$Re_T$ is given by (2.22). Note that, in the limiting case of very high-shear forcing, 
$Re_T \approx Re_S$. Equation (2.22) can be rewritten as
\begin{equation} \sqrt{\frac{2}{C_S}} = \frac{1}{k} \ln{ \left( Re_S \sqrt{\frac{C_S}{8}} \right) } + B, \end{equation}and (2.20) can be rewritten using (2.7) as
\begin{equation} Nu \sim Pr^{1/2} C_S Re_S, \end{equation}
which agrees well with the numerical data points in figures 3(c) and 3(d) at very high shear forcing. However, it is more useful to substitute the value of 
$C_T$ obtained from (2.22) into (2.20) and approximate 
$Re_T$ from (2.5) as
\begin{equation} Re_T \approx \sqrt{Re_R^2 + Re_S^2}, \end{equation} to obtain the dashed lines plotted in figure 3(a), which show better agreement for a larger range of 
$Re_S/Re_R$ in the shear-dominated regime. The approximation given by (3.6) is then validated in figure 4(a). The additional energy term 
$\gamma C_S Re_S^2$ from (2.4) that corresponds to the turbulent fluctuations arising from shear forcing is shown in figure 4(b) with the value of the prefactor 
$\gamma \approx 10.24$. Note that the contribution from the fluctuations is much smaller than the contribution from the mean streamwise velocity, thereby making (3.6) a good approximation.
Figure 4. (a) Value of 
$Re_T/Re_R$ plotted against 
$Re_S/Re_R$. The black solid line indicates the relation (2.5). (b) Value of 
$(\sqrt {Re_T^2 - Re_S^2})/Re_R$ plotted against 
$\sqrt {C_S} Re_S/Re_R$. The black dashed line indicates the relation (2.4) with 
$\gamma \approx 10.24$. The data for PRB are indicated with plus markers and the data for CRB are indicated with dot markers. This figure is also available at https://www.cambridge.org/S0022112024008723/JFM-Notebooks/files/figure4.
3.3. Large-scale convection rolls
Next, we confirm the theoretical assumptions on the LSC rolls made in § 2.2. The three-dimensional volume visualisation in figure 5 shows the streamlines associated with these LSC rolls and it can be seen that they are predominantly oriented in the spanwise direction. However, the large-scale temperature structures in figure 5 are seen to be aligned neither fully along the streamwise direction, nor fully along the spanwise direction but along a diagonal. The thermal plumes that comprise these large-scale structures experience the advective effects of both 
$U_L$ in the spanwise direction and 
$U_S$ in the streamwise direction. Therefore, the orientation of these large scale temperature flow structures in the 
$x\unicode{x2013}y$ plane is tilted along a diagonal whose slope is approximately given by 
$U_L/U_S$. This is made clearer in figure 6 through the visualisation of the non-dimensional time-averaged local shear stress 
$\boldsymbol {\tau _w^{\prime }} \equiv (\tau _x^{\prime },\tau _y^{\prime },0)$ computed at the bottom wall in the large aspect ratio CRB system by subtracting the wall-averaged streamwise shear stress in the following way:
\begin{equation} \tau_x^{\prime} = H U_{F}^{-1}\left\langle \partial_z u_x - \left\langle \partial_z u_x \right\rangle_{x,y,t}\right\rangle_{t}, \quad \tau_y^{\prime} = H U_{F}^{-1}\left\langle \partial_z u_y \right\rangle_{t}. \end{equation}
Figure 5. Three-dimensional volume visualisation of the time-averaged reduced temperature field of a CRB system with 
$Ra=10^7$, 
$Pr=1.0$, 
$Re_w=1414$, 
$\varGamma _x=48$, 
$\varGamma _y=24$ with the averaging time being 100 free-fall time units. The green curves indicate the streamlines of the time-averaged velocity field 
$\langle \boldsymbol {u}\rangle _{t} - U_w \hat {x}$. The spanwise reorientation and streamwise sweeping of the plumes are evident. The zoomed inset shows the time-averaged two-dimensional visualisation of the temperature and velocity vectors in the 
$x\unicode{x2013}z$ plane. The plumes carry streamwise momentum along with temperature and are swept in the streamwise direction, causing a reduction in the heat transport.
Figure 6. Visualisation of 
$\boldsymbol {\tau _w^{\prime }}$ for the CRB system with 
$Ra=10^7$ and 
$Pr=1$. Panels (c), (e) and (g) are in the buoyancy-dominated regime and panel (i) is in the shear-dominated regime. The colour bars indicate the magnitude of 
$\boldsymbol {\tau _w^{\prime }}$ while the black arrows in the magnified panels (b), (d), (f), (h) and (j) indicate the direction of 
$\boldsymbol {\tau _w^{\prime }}$. The black lines in (c–h) indicate the slope of 
$(Re_w$/
$Re_R)^{-1}$. This figure is also available at https://www.cambridge.org/S0022112024008723/JFM-Notebooks/files/figure6.
Figure 6(a) shows the wall shear for the RB system i.e. for 
$Re_w=0$ and the zoomed inset in figure 6(b) shows the vectors of 
$\boldsymbol {\tau _w^{\prime }}$. As expected, the LSC rolls are randomly oriented and no global alignment of 
$\boldsymbol {\tau _w^{\prime }}$ is observed. A visual inspection of figures 6(c), 6(e) and 6(g) reveals that the large-scale flow structures seem to be oriented along the diagonal whose slope is approximately given by 
$Re_R$/
$Re_w$. However, figures 6(d), 6(f) and 6(h) show that 
$\boldsymbol {\tau _w^{\prime }}$ is primarily oriented along the spanwise direction in the transitional regime. Figure 6(i) shows the breakdown of the LSC rolls and the formation of large meandering flow structures (Hutchins & Marusic Reference Hutchins and Marusic2007; Blass et al. Reference Blass, Zhu, Verzicco, Lohse and Stevens2020, Reference Blass, Tabak, Verzicco, Stevens and Lohse2021) in the shear-dominated regime, while figure 6(j) shows that 
$\boldsymbol {\tau _w^{\prime }}$ is predominantly aligned in the streamwise direction in the shear-dominated regime.
For further confirmation of the changes in the LSC rolls, we study the probability distribution function 
$\phi (\alpha )$ of the angle 
$\alpha$ spanned by the horizontal velocity component fluctuations 
$\boldsymbol {u_h^\prime } \equiv (u_x^{\prime }, u_y^{\prime }, 0)$ with the streamwise direction 
$x$. In figures 7(a)–7(e) it can be seen that the behaviour of 
$\phi (\alpha )$ is qualitatively quite similar for the CRB and PRB systems. In the RB system, the LSC rolls are randomly oriented as shown by a uniform 
$\phi (\alpha )$ in figure 7(a). In the buoyancy-dominated regime, the LSC roll has a strong tendency to align in the spanwise directions in the transitional regime, as shown in figures 7(b) and 7(d). In the shear-dominated regime, the velocity fluctuations are predominantly aligned in the streamwise direction as in the case of turbulent Couette/Poiseuille flows, as seen in figures 7(c) and 7(e). At the thermal BL height, the symmetry of 
$\phi (\alpha )$ about 
$\alpha =\pm {\rm \pi}/2$ is strongly suggested for all three regimes by the data from the numerical simulations as shown in figures 7(f)–7(h), confirming the assumption made in (2.6). In figure 7(f), the small non-uniformity in 
$\phi (\alpha )$ is attributed to the numerical confinement experienced by the flow structures in domains of a smaller aspect ratio of 
$\varGamma _x=8$ and 
$\varGamma _y=4$. For the relatively unconfined case with 
$\varGamma _x=48$ and 
$\varGamma _y=24$, the probability distribution is nearly uniform for all values of 
$\alpha$.
Figure 7. The probability distribution 
$\phi (\alpha )$ of the flow orientation angle 
$\alpha$ for all heights 
$0 \leq z/\lambda _{\theta } \leq Nu$ for 
$Ra=10^7$, 
$Pr=1$, 
$\varGamma _x=48$ and 
$\varGamma _y=24$ in (a) RB flow (i.e. with 
$Re_S=0)$, (b,d) buoyancy-dominated regime with 
$Re_S/Re_R\approx 2.11$ and (c,e) shear-dominated regime with 
$Re_S/Re_R\geq 10$. Panels (b,c) are for the CRB system while panels (d,e) are for the PRB system. The black dashed lines indicate the height of the thermal BL and the colour bar indicates the magnitude of the probability. The probability distribution 
$\phi (\alpha )$ at the thermal BL height plotted against angle 
$\alpha$ for various 
$Ra$ and 
$Pr$ in (f) RB flow, (g) buoyancy-dominated regime and (h) shear-dominated regime. The solid lines are the CRB system with 
$\varGamma _x=48$ and 
$\varGamma _y=24$, the dashed lines are for the CRB system with 
$\varGamma _x=8$ and 
$\varGamma _y=4$ and the dotted lines are for the PRB system with 
$\varGamma _x=8$ and 
$\varGamma _y=4$. This figure is also available at https://www.cambridge.org/S0022112024008723/JFM-Notebooks/files/figure7.
3.4. Dissipation rates
We now investigate the bulk and BL contributions to the global kinetic and thermal dissipation rates described in § 2.1. We validate the rigorous relations given by (2.1) and (2.2) using the data obtained from numerical simulations as shown in figures 8(a) and 8(b), respectively. For the range of 
$Ra$ and 
$Pr$ studied in this work, as long as the flow is buoyancy dominated, the system is in the 
$II_u$ regime with 
$\epsilon _u \sim \epsilon _{u,bulk}$. With increasing shear forcing, the shear term of (2.2) increases the BL contribution of the kinetic dissipation due to the formation of streamwise velocity gradients close to the wall. For sufficiently strong shear, the kinetic dissipation will be dominated by the BL contribution with 
$\epsilon _u \sim \epsilon _{u,BL}$ as shown in 8(c). On the contrary, figure 8(d) shows that the thermal dissipation rate is dominated by the BL contribution for the entire range of 
$Re_S$ considered in this study but the contribution reduces noticeably towards higher 
$Re_S$. For extremely strong shear, a possibility of a transition towards bulk dominance in thermal dissipation cannot be ruled out.
Figure 8. (a) Ratio 
$(H^2 \epsilon _\theta \kappa ^{-1} \varDelta ^2 )/Nu$ showing good agreement with (2.1) indicated using the solid black line, and (b) ratio 
$( H^4 \epsilon _u/ \nu ^3 ) / ( (Nu - 1) Ra Pr^{-2} + C_S Re_S^3)$ plotted against 
$Re_S$ showing good agreement with (2.2) of the manuscript shown using the solid black line. (c) Boundary layer contribution to the global kinetic dissipation rate 
$\epsilon _{u,BL}$ and (d) the BL contribution to the global thermal dissipation rate 
$\epsilon _{\theta,BL}$. The circle markers are for the CRB system, while plus markers are for the PRB system. This figure is also available at https://www.cambridge.org/S0022112024008723/JFM-Notebooks/files/figure8.
3.5. Friction coefficient in Rayleigh–Bénard flow
So far, we have considered only the laminar Prandtl–Blasius-type kinetic BLs in the buoyancy-dominated regime with the BL only becoming turbulent in highly shear-dominated regime. However, we now consider the possibility of turbulent kinetic BLs even in pure RB flow which corresponds to the so-called ‘ultimate’ regime in RB flow (Kraichnan Reference Kraichnan1962; Grossmann & Lohse Reference Grossmann and Lohse2011; Roche Reference Roche2020; Lohse & Shishkina Reference Lohse and Shishkina2023). While the existence of logarithmic temperature profiles has already been observed at high 
$Ra$ (Grossmann & Lohse Reference Grossmann and Lohse2012; Ahlers et al. Reference Ahlers, Bodenschatz, Funfschilling, Grossmann, He, Lohse, Stevens and Verzicco2012; Ahlers, Bodenschatz & He Reference Ahlers, Bodenschatz and He2014), it is yet to be seen if logarithmic behaviour given by
\begin{equation} u_R^+ (z^+) = \frac{1}{k} \log (z^+) + B, \end{equation}with
\begin{equation} u_R^+ (z^+) = u_R(z) /u_\tau, \quad z^+= \frac{z u_\tau}{2 \nu}, \quad u_R(z) \equiv \sqrt{\left\langle u_x^2 + u_y^2\right\rangle_{A,t}}, \quad u_\tau = \sqrt{\nu \left\langle\partial_z u_R\right\rangle_{W,t}}, \end{equation}can be observed in the velocity profiles.
At the highest thermal forcing studied in this work with 
$Ra=10^{12}$, we start to also observe some hints of what could possibly be the onset of a log layer (see figure 9a), although it cannot be conclusively confirmed with the currently available data. Assuming that such a log layer could exist, an estimate of the modified von Kármán constant 
$k$ is obtained from the inflection point of the diagnostic function plotted in figure 9(b), giving 
$k \approx 1$. Correspondingly, the intercept 
$B \approx 5.2$ is found by fitting the data, as shown in figure 9(a).
Figure 9. (a) Value of 
$u_R^+$ plotted against 
$z^+$, showing the onset of a log layer with the grey dashed line. The velocity profile for 
$Ra=10^{12}$, 
$Pr=1.0$ is indicated in magenta. (b) Diagnostic function plotted against 
$z^+$ with the inflection point and the corresponding von Kármán constant indicated. (c) Value of 
$\delta ^+$ plotted against 
$Re_R$ shows good agreement with 
$\delta ^+ = (1/2) Re_R^{1/2}$, indicated with the dashed grey line. (d) Value of 
$C_R$ plotted against 
$Re_\delta$ showing a collapse without any prefactors. The solid grey line indicates 
$C_R \sim 4Re_\delta ^{-2/3}$ and the dashed grey line indicates the modified Prandtl (Reference Prandtl1932) friction law given by (3.10a,b). (e) Value of 
$C_R/a$ plotted against 
$Re_R$, showing the Prandtl–Blasius scaling (2.18) at low 
$Re_R$ with the solid grey line and modified Prandtl (Reference Prandtl1932) friction law given by (3.10a,b) at high 
$Re_R$ with the dashed grey line. (f) Value of 
$C_\epsilon /b$ given by (3.11b) plotted against 
$Re_R$. The solid grey line indicates 
$C_\epsilon \sim Re_R^{-1/2}$ for the regime with 
$\epsilon _u \sim \epsilon _{u,BL}$, the dashed grey line indicates the modified Prandtl (Reference Prandtl1932) friction law (3.10a,b). At higher 
$Re_R$, 
$\epsilon _u \sim \epsilon _{u,bulk}$, leading to 
$C_\epsilon$ becoming independent of 
$Re_R$, as indicated with the dash-dotted and dotted grey lines. This figure is also available at https://www.cambridge.org/S0022112024008723/JFM-Notebooks/files/figure9.
In the presence of such a logarithmic layer, we can now hypothesise about a relation for 
$C_R$ which is analogous to relation (3.4) for 
$C_S$ and (2.22) for 
$C_T$ by stating that in the limit of highly turbulent BL, 
$u_R^+=U_R$ at 
$z^+=\delta ^+$, where 
$\delta$ is the wall-normal distance to the peak velocity shown in figure 4. This gives us
\begin{equation} \sqrt{\frac{2}{C_R}} = \frac{1}{k} \ln{ \left( Re_\delta \sqrt{\frac{C_R}{8}} \right) } + B, \quad Re_\delta = U_R \delta / \nu .\end{equation}
If applicable to high 
$Ra$ RB flow, the general form of this equation could suggest universality in the behaviour of wall-bounded flows even though the values of 
$k$ and 
$B$ might be different from those observed for pipe or channel flows. By plotting 
$\delta ^+$ against 
$Re_R$, we find a very good fit with the scaling 
$\delta ^+ \sim Re_R^{1/2}$, as shown in figure 9(c). At present, we can only provide this scaling empirically because more data at extremely high 
$Ra$ are needed to understand the dynamics of the turbulent BL which is difficult due to the high computational expense of such numerical simulations. Plotting 
$C_R$ against 
$Re_\delta$ in 9(d), we find a nice collapse of the data for different aspect ratios and geometries. At low 
$Re_R$, 
$C_R$ seems to scale as 
$Re_\delta ^{-2/3}$, which can be obtained from the empirical scaling 
$\delta ^+ \sim Re_R^{1/2}$ shown in figure 9(c) and the assumption of a Prandtl–Blasius-type scaling of 
$C_R \sim Re_R^{-1/2}$. At high 
$Re_R$, we see good qualitative agreement between the data and the modified Prandtl (Reference Prandtl1932) friction law given by (3.10a,b).
We can further investigate the behaviour of 
$C_R$ by computing
$$\begin{gather} C_R \equiv a Re_R^{-1/2}, \end{gather}$$
$$\begin{gather} C_\epsilon \equiv b (Nu - 1) Ra Pr^{-2} Re_R^{-3}, \quad \begin{cases} C_\epsilon \sim C_R & \text{if} \epsilon_u \sim \epsilon_{u,BL} \\ C_\epsilon = \text{const.} & \text{if} \epsilon_u \sim \epsilon_{u,bulk}, \end{cases} \end{gather}$$
where 
$a$ and 
$b$ are prefactors that are obtained by fitting the data from figures 9(e) and 9(f), accounting for the effects of aspect ratio and geometry. Note that, when the kinetic dissipation of the RB system is dominated by the contribution from the BL, 
$C_\epsilon \sim C_R$ in (3.11b). The data points at lower 
$Re_R$ in figure 9(e) are observed to follow (2.18), which is consistent with the Prandtl–Blasius scaling. For higher 
$Re_R$, corresponding to 
$Ra>10^{10}$, 
$C_R$ shows better agreement with the Prandtl (Reference Prandtl1932) friction law given by (3.10a,b). Figure 9(f) shows a very similar behaviour to figure 9(e) but, in addition, we observe that 
$C_\epsilon$ becomes independent of 
$Re_R$ at very high values of 
$Re_R$. Although this does not reflect the true dependence of 
$C_R$ on 
$Re_R$, this apparent dependence is expected because the kinetic dissipation of the RB system undergoes a transition from being dominated by the BL to being dominated by the bulk (Lohse Reference Lohse1994). It should also be noted that this transition occurs at higher 
$Re_R$ for the more confined cylindrical RB simulations because more kinetic driving is required to overcome the viscous dissipation in the additional BLs on the side walls that are not present in the periodic box RB simulations (Ahlers et al. Reference Ahlers2022).
4. Conclusions
In summary, we have developed a framework by extending the GL-theory for RB turbulence to sheared RB turbulence. As in the case of RB flow, we observe that there are no pure scaling exponents for the Nusselt number 
$Nu$ and the friction coefficient 
$C_S$. This also holds for high thermal or shear driving where the BLs no longer obey scaling relations associated with the Prandtl–Blasius (Prandtl Reference Prandtl1904; Blasius Reference Blasius1908) BL theory but start to become more turbulent. In such cases, we observe that the relation for 
$C_S(Re_S)$ is well described by the friction law of Prandtl (Reference Prandtl1932). In addition, we find that a modified version of the Reference PrandtlPrandtl's (Reference Prandtl1932) friction law for the large-scale circulation rolls 
$C_L(Re_L)$ analogous to 
$C_S(Re_S)$ agrees well with the DNS data, suggesting some sort of universality in the relation between the shear stress and the flow velocity that generates that shear.
It is also interesting to note that the relations are identical for CRB or PRB systems once the appropriate velocity scale is chosen as a control parameter. This suggests that the flow physics is not strongly affected by the geometry of the system or by the way in which shear forcing is applied. The flow characteristics of these systems are essentially determined by the ratio of shear driving to thermal driving, given by 
$Re_S/Re_R$. Applying shear to the RB system causes increased coherence in the streamwise direction and leads to a re-orientation of the LSC rolls and causes them to align more in the spanwise direction, with the thermal plumes also transporting the momentum imparted by the shear forcing. In the buoyancy-dominated regime with 
$Re_S \simeq Re_R$, this may lead to enhanced streamwise mixing between hot and cold plumes at a time scale that is smaller than the time scale of heat diffusion at the wall. This leads to heat entrapment in the bulk and a reduction in 
$Nu$.
Taking into account the orientation of the LSC rolls and the bulk dominance of 
$\epsilon _u$ in the buoyancy-dominated regime, we show that the orientation of large-scale flow structures can also be predicted to a reasonable degree by the ratio 
$Re_S/Re_R$, and we provide scaling relations for the 
$Nu(Re_S/Re_R)$ and 
$C_S(Re_S/Re_R)$, which are shown to agree well with the numerical simulations for 
$10^6 \leq Ra \leq 10^8$, 
$0.5 \leq Pr \leq 5.0$ and 
$0 \leq Re_S \leq 10^4$. However, the evidence from the DNS is limited at the moment due to its high computational costs, thereby restricting the parameter range in which the proposed scaling laws can be validated. Simulations for very high or very low values of 
$Pr$ as well as for high 
$Ra$ or 
$Re_S$ can be very demanding, and it remains to be seen if the assumptions made in this work and the extended theory hold well in other control parameter ranges.
Supplementary material
Computational Notebook files are available as supplementary material at https://doi.org/10.1017/jfm.2024.872 and online at https://www.cambridge.org/S0022112024008723/JFM-Notebooks.
Acknowledgements
The authors gratefully acknowledge R. Hartmann, A. Blass and M. De Paoli for fruitful discussions.
Funding
This work was financially supported by the ERC starting grant (2018) for the project ‘UltimateRB’ and the Twente Max-Planck Center. This work was carried out on the Dutch national e-infrastructure with the support of SURF Cooperative. We acknowledge PRACE for awarding us access to MareNostrum at Barcelona Supercomputing Center (BSC), Spain (Project 2020225335 and 2020235589). The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project by providing computing time on the GCS Supercomputer SuperMUC at Leibniz Supercomputing Centre (www.lrz.de).
Declaration of interests
The authors report no conflict of interest.
Author contributions
The theoretical framework was developed by D.L., O.S., R.J.A.M.S., C.J.H. and G.S.Y. The code ‘AFiD’ for the DNS was developed by R.V. while the numerical simulations were performed by G.S.Y. All authors were involved in the preparation of the manuscript.
Appendix. Simulation parameters
In table 3 we provide the physical and numerical input parameters used for the new sheared RB simulations conducted for this study. In addition to the new simulations, data for large aspect ratio simulations with 
$\varGamma _x = 48$ and 
$\varGamma _y = 24$ are taken from Yerragolam et al. (Reference Yerragolam, Verzicco, Lohse and Stevens2022b), and data for high 
$Ra$ cylindrical simulations are taken from Hartmann et al. (Reference Hartmann, Yerragolam, Verzicco, Lohse and Stevens2023).
Table 3. Simulations considered in this work. The aspect ratios of the computational domain are given by 
$\varGamma _x$ in the streamwise direction and 
$\varGamma _y$ in the spanwise direction. The values of 
$N_x$, 
$N_y$ and 
$N_z$ indicate the number of grid points in the streamwise, spanwise and wall-normal directions. The grid spacing in wall units in the streamwise and spanwise directions is given by 
$\Delta x^+$ and 
$\Delta y^+$, respectively. The wall-normal grid spacing in wall units at the wall and the mid-height are given by 
$\Delta z_w^+$ and 
$\Delta z_c^+$, respectively. The number of grid points in the wall-normal direction within the thermal BL is given by 
$N_{BL}$. Here, 
$Nu_w$ indicates the value of Nusselt number computed using the average gradient of reduced temperature profile at the walls using (1.5) and 
$Nu_{\epsilon _\theta }$ and 
$Nu_{\epsilon _u}$ indicate the Nusselt numbers that satisfy the global balance of thermal and kinetic dissipation (2.1) and (2.2).
Open access
